svd.cc

basic linear algebra classes and applications (SVD,interpolation, multivariate optimization)

 * consequence of Frantsis' Theorem, see the book above. In that case, we can
 * eliminate the N-th row and column of JJ and carry out further transforms
 * with a smaller matrix. If any other superdiag element of JJ turns 0,
 * then JJ effectively falls into two independent matrices. We will process
 * them independently (the bottom one first).
 *
 * Since matrix J is a bidiagonal, it can be stored efficiently. As a matter
 * of fact, its N diagonal elements are in array Sig, and superdiag elements
 * are stored in array super_diag.
 */
 
				// Carry out U * S with a rotation matrix
				// S (which combines i-th and j-th columns
				// of U, i>j)
inline void SVD::rotate(Matrix& U, const int i, const int j,
		   const double cos_ph, const double sin_ph)
{
  LAStreamOut Ui(MatrixColumn (U,i));
  LAStreamOut Uj(MatrixColumn (U,j));
  while( !Ui.eof() )
  {
    REAL& Uil = Ui.get(); REAL& Ujl = Uj.get();
    const REAL Ujl_was = Ujl;
    Ujl =  cos_ph*Ujl_was + sin_ph*Uil;
    Uil = -sin_ph*Ujl_was + cos_ph*Uil;
  }
}

/*
 * A diagonal element J[l-1,l-1] turns out 0 at the k-th step of the
 * algorithm. That means that one of the original matrix' singular numbers
 * shall be zero. In that case, we multiply J by specially selected
 * matrices S' of simple rotations to eliminate a superdiag element J[l-1,l].
 * After that, matrix J falls into two pieces, which can be dealt with
 * in a regular way (the bottom piece first).
 * 
 * These special S transformations are accumulated into matrix U: since J
 * is left-multiplied by S', U would be right-multiplied by S. Transform
 * formulas for doing these rotations are similar to (5) above. See the
 * book cited above for more details.
 */
inline void SVD::rip_through(
	Vector& super_diag, const int k, const int l, const double eps)
{
  double cos_ph = 0, sin_ph = 1;	// Accumulate cos,sin of Ph
  				// The first step of the loop below
  				// (when i==l) would eliminate J[l-1,l],
  				// which is stored in super_diag(l)
  				// However, it gives rise to J[l-1,l+1]
  				// and J[l,l+2]
  				// The following steps eliminate these
  				// until they fall below
  				// significance
  for(register int i=l; i<=k; i++)
  {
    const double f = sin_ph * super_diag(i);
    super_diag(i) *= cos_ph;
    if( abs(f) <= eps )
      break;			// Current J[l-1,l] became unsignificant
    cos_ph = sig(i); sin_ph = -f;	// unnormalized sin/cos
    const double norm = (sig(i) = hypot(cos_ph,sin_ph)); // sqrt(sin^2+cos^2)
    cos_ph /= norm, sin_ph /= norm;	// Normalize sin/cos
    rotate(U,i,l-1,cos_ph,sin_ph);
  }
}

			// We're about to proceed doing QR-transforms
			// on a (bidiag) matrix J[1:k,1:k]. It may happen
			// though that the matrix splits (or can be
			// split) into two independent pieces. This function
			// checks for splitting and returns the lowerbound
			// index l of the bottom piece, J[l:k,l:k]
inline int SVD::get_submatrix_to_work_on(
	Vector& super_diag, const int k, const double eps)
{
  for(register int l=k; l>1; l--)
    if( abs(super_diag(l)) <= eps )
      return l;				// The breaking point: zero J[l-1,l]
    else if( abs(sig(l-1)) <= eps )	// Diagonal J[l,l] turns out 0
    {					// meaning J[l-1,l] _can_ be made
      rip_through(super_diag,k,l,eps);	// zero after some rotations
      return l;
    }
  return 1;			// Deal with J[1:k,1:k] as a whole
}

		// Diagonalize root module
void SVD::diagonalize(Vector& super_diag, const double eps)
{
  for(register int k=N; k>=1; k--)	// QR-iterate upon J[l:k,l:k]
  {
    register int l;
    while(l=get_submatrix_to_work_on(super_diag,k,eps),
    	  abs(super_diag(k)) > eps)	// until superdiag J[k-1,k] becomes 0
    {
      double shift;			// Compute a QR-shift from a bottom
      {					// corner minor of J[l:k,l:k] order 2
      	REAL Jk2k1 = super_diag(k-1),	// J[k-2,k-1]
      	     Jk1k  = super_diag(k),
      	     Jk1k1 = sig(k-1),		// J[k-1,k-1]
      	     Jkk   = sig(k),
      	     Jll   = sig(l);		// J[l,l]
      	shift = (Jk1k1-Jkk)*(Jk1k1+Jkk) + (Jk2k1-Jk1k)*(Jk2k1+Jk1k);
      	shift /= 2*Jk1k*Jk1k1;
      	shift += (shift < 0 ? -1 : 1) * sqrt(shift*shift+1);
      	shift = ( (Jll-Jkk)*(Jll+Jkk) + Jk1k*(Jk1k1/shift-Jk1k) )/Jll;
      }
      				// Carry on multiplications by T2, S2, T3...
      double cos_th = 1, sin_th = 1;
      REAL Ji1i1 = sig(l);	// J[i-1,i-1] at i=l+1...k
      for(register int i=l+1; i<=k; i++)
      {
      	REAL Ji1i = super_diag(i), Jii = sig(i);  // J[i-1,i] and J[i,i]
      	sin_th *= Ji1i; Ji1i *= cos_th; cos_th = shift;
      	double norm_f = (super_diag(i-1) = hypot(cos_th,sin_th));
      	cos_th /= norm_f, sin_th /= norm_f;
      					// Rotate J[i-1:i,i-1:i] by Ti
      	shift = cos_th*Ji1i1 + sin_th*Ji1i;	// new J[i-1,i-1]
      	Ji1i = -sin_th*Ji1i1 + cos_th*Ji1i;	// J[i-1,i] after rotation
      	const double Jii1 = Jii*sin_th;		// Emerged J[i,i-1]
      	Jii *= cos_th;				// new J[i,i]
        rotate(V,i,i-1,cos_th,sin_th); // Accumulate T rotations in V
        
        double cos_ph = shift, sin_ph = Jii1;// Make Si to get rid of J[i,i-1]
        sig(i-1) = (norm_f = hypot(cos_ph,sin_ph));	// New J[i-1,i-1]
        if( norm_f == 0 )		// If norm =0, rotation angle
          cos_ph = cos_th, sin_ph = sin_th; // can be anything now
        else
          cos_ph /= norm_f, sin_ph /= norm_f;
      					// Rotate J[i-1:i,i-1:i] by Si
        shift = cos_ph * Ji1i + sin_ph*Jii;	// New J[i-1,i]
        Ji1i1 = -sin_ph*Ji1i + cos_ph*Jii;	// New Jii, would carry over
        					// as J[i-1,i-1] for next i
        rotate(U,i,i-1,cos_ph,sin_ph);  // Accumulate S rotations in U
        				// Jii1 disappears, sin_th would
        cos_th = cos_ph, sin_th = sin_ph; // carry over a (scaled) J[i-1,i+1]
        				// to eliminate on the next i, cos_th
        				// would carry over a scaled J[i,i+1]
      }
      super_diag(l) = 0;		// Supposed to be eliminated by now
      super_diag(k) = shift;
      sig(k) = Ji1i1;
    }		// --- end-of-QR-iterations
    if( sig(k) < 0 )		// Correct the sign of the sing number
    {
      sig(k) = -sig(k);
      for(LAStreamOut Vk(MatrixColumn(V,k)); !Vk.eof(); )
        { REAL& vk = Vk.get(); vk = -vk; }
    }
  }
    
} 


/*
 *------------------------------------------------------------------------
 *				The root Module
 */

SVD::SVD(const Matrix& A)
   : M(A.q_nrows()), N(A.q_ncols()),
     U(A.q_nrows(),A.q_nrows()),
     V(A.q_ncols(),A.q_ncols()),
     sig(A.q_ncols())
{
  if( M < N )
    A.info(),
    _error("Matrix A should have at least as many rows as it has columns");
     
  U.unit_matrix(); V.unit_matrix();

  Vector super_diag(N);
  const double bidiag_norm = bidiagonalize(super_diag,A);
  const double eps = FLT_EPSILON * bidiag_norm;	// Significance threshold
  diagonalize(super_diag,eps);
}

/*
 *------------------------------------------------------------------------
 *		Print some info about the SVD that has been built
 */

				// Print the info about the SVD
void SVD::info(void) const
{
  U.is_valid();
  message("\nSVD of an %dx%d matrix",M,N);
}

				// Return min and max singular values
SVD::operator MinMax(void) const
{
  LAStreamIn sigs(sig);
  MinMax mm(sigs.get());
  while( !sigs.eof() )
    mm << sigs.get();
  return mm;
}

					// sig_max/sig_min
double SVD::q_cond_number(void) const
{
  return ((MinMax)(*this)).ratio();
}

/*
 *------------------------------------------------------------------------
 * 			class SVD_inv_mult
 * 	Solving an (overspecified) set of linear equations A*X=B
 * 		using a least squares method via SVD
 *
 * Matrix A is a rectangular M*N matrix with M>=N, B is a M*K matrix
 * with K either =1 (that is, B is merely a column-vector) or K>1.
 * If B is a unit matrix of the size that of A, the present LazyMatrix is
 * a regularized (pseudo)inverse of A.
 *
 * Algorithm
 *   Matrix A is decomposed first using SVD:
 *	(1) A = U*Sig*V'
 * where matrices U and V are orthogonal and Sig is diagonal.
 * The set of simultaneous linear equations AX=B can be written then as
 *	(2) Sig*V'*X = U'*B
 * (where we have used the fact that U'=inv(U)), or
 *	(3) Sig*Vx = Ub
 * where we introduced
 *	(4) Vx = V'*X and Ub = U'*b
 * Since Sig is a diag matrix, eq. (3) is solved trivially:
 * 	(5) Vx[i] = Ub[i]/sig[i], if sig[i] > tau
 *		  = 0	otherwise
 * In the latter case matrix Sig has an incomplete rank, or close to that.
 * That is, one or more singular values are _too_ small. This means that
 * there is linear dependence among the equations of the set. In that case,
 * we will print the corresponding "null coefficients", the corresponding
 * columns of V. Adding them to the solution X won't change A*X.
 *
 * Having computed Vx, X is simply recovered as
 *	(6) X = V*Vx
 * since V*V' = E.
 *
 * Threshold tau in (5) is chosen as N*FLT_EPSILON*max(Sig[i,i]) unless
 * specified otherwise by the user.
 */

SVD_inv_mult::SVD_inv_mult
	(const SVD& _svd, const Matrix& _B,const double _tau)
	: LazyMatrix(_svd.q_V().q_ncols(),_B.q_ncols()),
	  svd(_svd), B(_B), tau(_tau)
{
  if( svd.q_U().q_nrows() != B.q_nrows() )
    svd.info(), B.info(),
    _error("Unsuitable matrices for SVD*X=B set");
  MinMax sig_minmax = (MinMax)svd;
  if( tau == 0 )
   tau = svd.q_V().q_nrows() * sig_minmax.max() * FLT_EPSILON;
  are_zero_coeff = sig_minmax.min() < tau;
}

#if 0
				// Computing X in a special case where
				// B (and X) is a vector
inline void SVD_inv_mult::fill_in_vector(Matrix& X) const
{
  X.clear();
  bool x_is_vector = x.q_ncols() == 1;
  if( are_zero_coeff )
    message("\nSVD solver of AX=B detected a linear dependency among X"
    	    "\n  #  \tsingular value\tnull coefficients\n");
  for(register int i=1; i<=svd.q_V().q_nrows(); i++)
  {
     MatrixCol Ui(U,i);
     const double sigi = svd.q_sig()(i);
     if( sigi > tau )
       if( x_is_vector )
         X(i) = (Ui * B)/sigi;
       else
         X(i) =0, print null coeffs...
  }
  X *= V;			// ???
}
#endif
				// Computing X in a general case where
				// B (and X) is a rectangular matrix
void SVD_inv_mult::fill_in(Matrix& X) const
{
  if( are_zero_coeff )
    message("\nSVD solver of AX=B detected a linear dependency among X"
    	    "\n  #  singular value\tnull coefficients\n");
  Matrix Vx = zero(X);
  Vector Ui(svd.q_U().q_nrows());			// i-th col of U
  Vector Bj(B.q_nrows());				// j-th col of B
  const Matrix& V = svd.q_V();
  const Matrix& U = svd.q_U();
  LAStreamIn sigs(svd.q_sig());
  for(register int i=1; i<=V.q_ncols(); i++)
  {
    const double sigi = sigs.get();
    LAStrideStreamOut Vxi( MatrixRow(Vx,i) );
    if( sigi > tau )
    {
      to_every(Ui) = of_every(ConstMatrixColumn(U,i));
      for(register int j=1; j<=X.q_ncols(); j++)
        to_every(Bj) = of_every(ConstMatrixColumn(B,j)), Vxi.get() = Ui * Bj/sigi;
    }
    else
    {
      while( !Vxi.eof() ) Vxi.get() = 0;
      message(" %d %12.2g \t(",i,sigi);
      LAStreamIn Vi(ConstMatrixColumn(V,i));
      while( !Vi.eof() )
	message("%10g ",Vi.get());
      message(")\n");
    }
  }
  X.mult(V,Vx);				// that is, set X to be V*Vx
}