Go to the documentation of this file.

/* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4 This file is part of the Feel library Author(s): Christophe Prud'homme <christophe.prudhomme@ujf-grenoble.fr> Date: 2005-03-27 Copyright (C) 2005,2006 EPFL This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 3.0 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ /** \file svd.hpp \author Christophe Prud'homme <christophe.prudhomme@ujf-grenoble.fr> \date 2005-03-27 */ #ifndef __SVD_HPP #define __SVD_HPP 1 #include <cmath> #include <limits> #include <algorithm> #include <boost/numeric/ublas/matrix.hpp> #include <boost/numeric/ublas/vector.hpp> #include <boost/numeric/ublas/io.hpp> #include <feel/feelcore/feel.hpp> #include <feel/feelcore/traits.hpp> namespace Feel { namespace ublas = boost::numeric::ublas; /*! \class SVD \brief Singular Value Decomposition of a rectangular matrix \f$ A = U * S * V^T \f$ where matrices \f$ U\f$ and \f$ V\f$ are orthogonal and Sig is a digonal matrix. The singular value decomposition is performed by constructing an SVD object from an M*N matrix A with M>=N (that is, at least as many rows as columns). Note, in case M > N, matrix Sig has to be a M*N diagonal matrix. However, it has only N diag elements, which we store in a 1:N Vector sig. Algorithm Bidiagonalization with Householder reflections followed by a modification of a QR-algorithm. For more details, see G.E. Forsythe, M.A. Malcolm, C.B. Moler Computer methods for mathematical computations. - Prentice-Hall, 1977 However, in the present implementation, matrices U and V are computed right away rather than delayed until after all Householder reflections. @author Christophe Prud'homme @see */ template<typename MatrixA> 00073 class SVD { public: /** @name Typedefs */ //@{ typedef typename MatrixA::value_type value_type; typedef MatrixA matrix_type; typedef boost::numeric::ublas::vector<value_type> vector_type; //@} /** @name Constructors, destructor */ //@{ SVD( matrix_type const& __A ); ~SVD() {} //@} /** @name Operator overloads */ //@{ //@} /** @name Accessors */ //@{ /** \return U of the Singular Value Decomposition */ 00109 matrix_type const& U() const { return _M_U; } /** \return V of the Singular Value Decomposition */ 00112 matrix_type const& V() const { return _M_V; } /** \return S of the Singular Value Decomposition */ 00115 vector_type const& S() const { return _M_S; } /** Computes \f$ C = \frac{\sigma^\mathrm{max}}{\sigma^\mathrm{min}} \f$ which represents the condition number of the matrix that was decomposed in singular values \return $C$ the condition number of the matrix */ 00123 value_type conditionNumber() const { return max(_M_S)/min(_M_S); } void conditionNumber( value_type& __max, value_type& __min ) const { __max = *std::max_element(_M_S.begin(), _M_S.end()); __min = *std::min_element(_M_S.begin(), _M_S.end()); } //@} /** @name Mutators */ //@{ //@} /** @name Methods */ //@{ //@} protected: private: SVD( SVD const & ); /** \brief Left Householder Transformations Zero out an entire subdiagonal of the i-th column of A and compute the modified A[i,i] by multiplication (UP' * A) with a matrix UP' (1) UP' = E - UPi * UPi' / beta where a column-vector UPi is as follows (2) UPi = [ (i-1) zeros, A[i,i] + Norm, vector A[i+1:M,i] ] where beta = UPi' A[,i] and Norm is the norm of a vector A[i:M,i] (sub-diag part of the i-th column of A). Note we assign the Norm the same sign as that of A[i,i]. By construction, (1) does not affect the first (i-1) rows of A. Since A[*,1:i-1] is bidiagonal (the result of the i-1 previous steps of the bidiag algorithm), transform (1) doesn't affect these i-1 columns either as one can easily verify. The i-th column of A is transformed as (3) UP' A[*,i] = A[*,i] - UPi (since UPi'*A[*,i]/beta = 1 by construction of UPi and beta) This means effectively zeroing out A[i+1:M,i] (the entire subdiagonal of the i-th column of A) and replacing A[i,i] with the -Norm. Thus modified A[i,i] is returned by the present function. The other (i+1:N) columns of A are transformed as (4) UP' A[,j] = A[,j] - UPi * ( UPi' * A[,j] / beta ) Note, due to (2), only elements of rows i+1:M actually participate in above transforms; the upper i-1 rows of A are not affected. As was mentioned earlier, (5) beta = UPi' * A[,i] = (A[i,i] + Norm)*A[i,i] + A[i+1:M,i]*A[i+1:M,i] = ||A[i:M,i]||^2 + Norm*A[i,i] = Norm^2 + Norm*A[i,i] (note the sign of the Norm is the same as A[i,i]) For extra precision, vector UPi (and so is Norm and beta) are scaled, which would not affect (4) as easy to see. To satisfy the definition (6) .SIG = U' * A * V the result of consecutive transformations (1) over matrix A is accumulated in matrix U' (which is initialized to be a unit matrix). At each step, U' is left-multiplied by UP' = UP (UP is symmetric by construction, see (1)). That is, U is right-multiplied by UP, that is, rows of U are transformed similarly to columns of A, see eq. (4). We also keep in mind that multiplication by UP at the i-th step does not affect the first i-1 columns of U. Note that the vector UPi doesn't have to be allocated explicitly: its first i-1 components are zeros (which we can always imply in computations), and the rest of the components (but the UPi[i]) are the same as those of A[i:M,i], the subdiagonal of A[,i]. This column, A[,i] is affected only trivially as explained above, that is, we don't need to carry this transformation explicitly (only A[i,i] is going to be non-trivially affected, that is, replaced by -Norm, but we will use sig[i] to store the result). */ value_type leftHouseholder( MatrixA& A, const int i ); /** \brief Right Householder Transformations Zero out i+2:N elements of a row A[i,] of matrix A by right multiplication (A VP) with a matrix VP (1) VP = E - VPi VPi' / beta where a vector-column .VPi is as follows (2) VPi = [ i zeros, A[i,i+1] + Norm, vector A[i,i+2:N] ] where beta = A[i,] VPi and Norm is the norm of a vector A[i,i+1:N] (right-diag part of the i-th row of A). Note we assign the Norm the same sign as that of A[i,i+1]. By construction, (1) does not affect the first i columns of A. Since A[1:i-1,] is bidiagonal (the result of the previous steps of the bidiag algorithm), transform (1) doesn't affect these i-1 rows either as one can easily verify. The i-th row of A is transformed as (3) A[i,*] VP = A[i,*] - VPi' (since A[i,*]*VPi/beta = 1 by construction of VPi and beta) This means effectively zeroing out A[i,i+2:N] (the entire right super- diagonal of the i-th row of A, but ONE superdiag element) and replacing A[i,i+1] with - Norm. Thus modified A[i,i+1] is returned as the result of the present function. The other (i+1:M) rows of A are transformed as (4) A[j,] VP = A[j,] - VPi' ( A[j,] VPi / beta ) Note, due to (2), only elements of columns i+1:N actually participate in above transforms; the left i columns of A are not affected. As was mentioned earlier, (5) beta = A[i,] VPi = (A[i,i+1] + Norm)*A[i,i+1]+ A[i,i+2:N]*A[i,i+2:N] = ||A[i,i+1:N]||^2 + Norm*A[i,i+1] = Norm^2 + Norm*A[i,i+1] (note the sign of the Norm is the same as A[i,i+1]) For extra precision, vector VPi (and so is Norm and beta) are scaled, which would not affect (4) as easy to see. The result of consecutive transformations (1) over matrix A is accumulated in matrix V (which is initialized to be a unit matrix). At each step, V is right-multiplied by VP. That is, rows of V are transformed similarly to rows of A, see eq. (4). We also keep in mind that multiplication by VP at the i-th step does not affect the first i rows of V. Note that vector VPi doesn't have to be allocated explicitly: its first i components are zeros (which we can always imply in computations), and the rest of the components (but the VPi[i+1]) are the same as those of A[i,i+1:N], the superdiagonal of A[i,]. This row, A[i,] is affected only trivially as explained above, that is, we don't need to carry this transformation explicitly (only A[i,i+1] is going to be non-trivially affected, that is, replaced by -Norm, but we will use super_diag[i+1] to store the result). */ value_type rightHouseholder( MatrixA& A, const int i ); /** \brief Bidiagonalization This nethod turns matrix A into a bidiagonal one. Its N diagonal elements are stored in a vector sig, while N-1 superdiagonal elements are stored in a vector super_diag(2:N) (with super_diag(1) being always 0). Matrices U and V store the record of orthogonal Householder reflections that were used to convert A to this form. The method returns the norm of the resulting bidiagonal matrix, that is, the maximal column sum. */ value_type bidiagonalize( vector_type& super_diag, const matrix_type& _A); /* \brief QR-diagonalization of a bidiagonal matrix After bidiagonalization we get a bidiagonal matrix J: (1) J = U' * A * V The present method turns J into a matrix JJ by applying a set of orthogonal transforms (2) JJ = S' * J * T Orthogonal matrices S and T are chosen so that JJ were also a bidiagonal matrix, but with superdiag elements smaller than those of J. We repeat (2) until non-diag elements of JJ become smaller than EPS and can be disregarded. Matrices S and T are constructed as (3) S = S1 * S2 * S3 ... Sn, and similarly T where Sk and Tk are matrices of simple rotations (4) Sk[i,j] = i==j ? 1 : 0 for all i>k or i<k-1 Sk[k-1,k-1] = cos(Phk), Sk[k-1,k] = -sin(Phk), SK[k,k-1] = sin(Phk), Sk[k,k] = cos(Phk), k=2..N Matrix Tk is constructed similarly, only with angle Thk rather than Phk. Thus left multiplication of J by SK' can be spelled out as (5) (Sk' * J)[i,j] = J[i,j] when i>k or i<k-1, [k-1,j] = cos(Phk)*J[k-1,j] + sin(Phk)*J[k,j] [k,j] = -sin(Phk)*J[k-1,j] + cos(Phk)*J[k,j] That is, k-1 and k rows of J are replaced by their linear combinations; the rest of J is unaffected. Right multiplication of J by Tk similarly changes only k-1 and k columns of J. Matrix T2 is chosen the way that T2'J'JT2 were a QR-transform with a shift. Note that multiplying J by T2 gives rise to a J[2,1] element of the product J (which is below the main diagonal). Angle Ph2 is then chosen so that multiplication by S2' (which combines 1 and 2 rows of J) gets rid of that elemnent. But this will create a [1,3] non-zero element. T3 is made to make it disappear, but this leads to [3,2], etc. In the end, Sn removes a [N,N-1] element of J and matrix S'JT becomes bidiagonal again. However, because of a special choice of T2 (QR-algorithm), its non-diag elements are smaller than those of J. All this process in more detail is described in J.H. Wilkinson, C. Reinsch. Linear algebra - Springer-Verlag, 1971 If during transforms (1), JJ[N-1,N] turns 0, then JJ[N,N] is a singular number (possibly with a wrong (that is, negative) sign). This is a 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) */ void rotate( matrix_type& U, const int i, const int j, const value_type cos_ph, const value_type sin_ph ); /* 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. */ void ripThrough( vector_type& super_diag, const int k, const int l, const double eps ); /** 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] */ int getWorkSubmatrix( vector_type& super_diag, const int k, const double eps ); /** main algorithm */ void diagonalize( vector_type& super_diag, const double eps ); private: uint _M_rows; uint _M_cols; matrix_type _M_U; matrix_type _M_V; vector_type _M_S; }; template<typename MatrixA> 00398 class SOrth { public: /** @name Typedefs */ //@{ typedef typename MatrixA::value_type value_type; typedef MatrixA matrix_type; typedef boost::numeric::ublas::vector<value_type> vector_type; //@} /** @name Constructors, destructor */ //@{ SOrth( matrix_type const& __A ); ~SOrth() {} //@} /** @name Operator overloads */ //@{ //@} /** @name Accessors */ //@{ /** \return U of the Singular Value Decomposition */ 00433 matrix_type const& Q() const { return _M_Q; } //@} /** @name Mutators */ //@{ //@} /** @name Methods */ //@{ //@} private: uint _M_rows; uint _M_cols; matrix_type _M_Q; }; /** \todo unfinished implementation: need to extract */ template<typename MatrixA> SOrth<MatrixA>::SOrth( MatrixA const& A ) : _M_rows( A.size1() ), _M_cols( A.size2() ), _M_Q() { using namespace boost::numeric::ublas; SVD<matrix_type> __svd( A ); vector_type __S; if ( _M_rows > 1 ) __S = __svd.S(); else if ( _M_rows == 1 ) __S(0) = __svd.S()( 0 ); const value_type eps = std::numeric_limits<value_type>::epsilon() * *std::max_element( __S.begin(), __S.end() )* std::max( _M_rows, _M_cols ); /*[m,n] = size(A); if m > 1, s = diag(S); elseif m == 1, s = S(1); else s = 0; end tol = max(m,n) * max(s) * eps; r = sum(s > tol); Q = U(:,1:r);*/ } /* * SVD */ template<typename MatrixA> SVD<MatrixA>::SVD( MatrixA const& A ) : _M_rows( A.size1() ), _M_cols( A.size2() ), _M_U( _M_rows, _M_rows ), _M_V( _M_cols, _M_cols ), _M_S( _M_cols ) { //FEEL_ASSERT( _M_rows >= _M_cols ).error("The Matrix A should have at least as many rows as it has columns"); MatrixA B( A ); bool do_transpose = _M_rows < _M_cols; if ( do_transpose ) { // transpose _M_rows = A.size2(); _M_cols = A.size1(); _M_U.resize( _M_rows, _M_rows ); _M_V.resize( _M_cols, _M_cols ); _M_S.resize( _M_cols ); B = ublas::trans( A ); } _M_U = boost::numeric::ublas::identity_matrix<value_type>( _M_rows ); _M_V = boost::numeric::ublas::identity_matrix<value_type>( _M_cols ); vector_type super_diag( _M_cols ); const value_type bidiag_norm = bidiagonalize(super_diag,B); // define Significance threshold const value_type eps = std::numeric_limits<value_type>::epsilon() * bidiag_norm; diagonalize(super_diag,eps); if ( do_transpose ) { matrix_type Temp( _M_V ); _M_V = ublas::trans( _M_U ); _M_U = ublas::trans( Temp ); } } /* Bidiagonalization */ template<typename MatrixA> typename SVD<MatrixA>::value_type SVD<MatrixA>::leftHouseholder( matrix_type& A, const int i ) { value_type scale = 0; // Compute the scaling factor for(int k=i; k< _M_rows; k++) scale += math::abs(A(k,i)); if( scale == 0 ) // If A[,i] is a null vector, no return 0; // transform is required value_type Norm_sqr = 0; // Scale UPi (that is, A[,i]) for(int k=i; k< _M_rows; k++) // and compute its norm, Norm^2 { A(k,i) /= scale; Norm_sqr += A(k,i)*A(k,i); } value_type new_Aii = math::sqrt(Norm_sqr); // new_Aii = -Norm, Norm has the if( A(i,i) > 0 ) new_Aii = -new_Aii; // same sign as Aii (that is, UPi[i]) const value_type beta = - A(i,i)*new_Aii + Norm_sqr; A(i,i) -= new_Aii; // UPi[i] = A[i,i] - (-Norm) for(int j=i+1; j<_M_cols; j++) // Transform i+1:N columns of A { value_type factor = 0; for(int k=i; k< _M_rows; k++) factor += A(k,i) * A(k,j); factor /= beta; for(int k=i; k< _M_rows; k++) A(k,j) -= A(k,i) * factor; } for(size_type j=0; j< _M_rows; j++) // Accumulate the transform in U { value_type factor = 0; for(size_type k=i; k< _M_rows; k++) factor += A(k,i) * _M_U(j,k); // Compute U[j,] * UPi factor /= beta; for(size_type k=i; k< _M_rows; k++) _M_U(j,k) -= A(k,i) * factor; } return new_Aii * scale; // Scale new Aii back (our new Sig[i]) } template<typename MatrixA> typename SVD<MatrixA>::value_type SVD<MatrixA>::rightHouseholder( matrix_type& A, const int i ) { value_type scale = 0; // Compute the scaling factor for(uint16_type k=i+1; k<_M_cols; k++) scale += math::abs(A(i,k)); if( scale == 0 ) // If A[i,] is a null vector, no return 0; // transform is required value_type Norm_sqr = 0; // Scale VPi (that is, A[i,]) for(uint16_type k=i+1; k<_M_cols; k++) // and compute its norm, Norm^2 { A(i,k) /= scale; Norm_sqr += A(i,k)*A(i,k); } value_type new_Aii1 = math::sqrt(Norm_sqr); // new_Aii1 = -Norm, Norm has the if( A(i,i+1) > 0 ) // same sign as new_Aii1 = -new_Aii1; // Aii1 (that is, VPi[i+1]) const value_type beta = - A(i,i+1)*new_Aii1 + Norm_sqr; A(i,i+1) -= new_Aii1; // VPi[i+1] = A[i,i+1] - (-Norm) for(uint16_type j=i+1; j<_M_rows; j++) // Transform i+1:M rows of A { value_type factor = 0; for(uint16_type k=i+1; k<_M_cols; k++) factor += A(i,k) * A(j,k); // Compute A[j,] * VPi factor /= beta; for(uint16_type k=i+1; k<_M_cols; k++) A(j,k) -= A(i,k) * factor; } for(uint16_type j=0; j<_M_cols; j++) // Accumulate the transform in V { value_type factor = 0; for(uint16_type k=i+1; k<_M_cols; k++) factor += A(i,k) * _M_V(j,k); // Compute V[j,] * VPi factor /= beta; for(uint16_type k=i+1; k<_M_cols; k++) _M_V(j,k) -= A(i,k) * factor; } return new_Aii1 * scale; // Scale new Aii1 back } template<typename MatrixA> typename SVD<MatrixA>::value_type 00629 SVD<MatrixA>::bidiagonalize( vector_type& super_diag, const matrix_type& _A) { value_type __norm_acc = 0; // No superdiag elem above A(0,0) super_diag(0) = 0; // A being transformed matrix_type A = _A; for(uint16_type __i = 0; __i < _M_cols; ++__i) { _M_S(__i) = leftHouseholder(A,__i); //std::cout << "S=" << _M_S << "\n"; const value_type& diagi = _M_S( __i ); if( __i < _M_cols-1 ) super_diag(__i+1) = rightHouseholder(A,__i); __norm_acc = std::max(__norm_acc, math::abs(diagi) + math::abs(super_diag(__i))); } return __norm_acc; } template<typename MatrixA> void SVD<MatrixA>::rotate( matrix_type& U, const int i, const int j, const value_type cos_ph, const value_type sin_ph ) { using namespace boost::numeric::ublas; matrix_column<matrix_type> Ui ( U, i); matrix_column<matrix_type> Uj ( U, j); for ( unsigned __k = 0; __k < Ui.size(); ++ __k ) { const value_type Ujk_was = Uj(__k); Uj( __k ) = cos_ph * Ujk_was + sin_ph * Ui(__k); Ui( __k ) = -sin_ph * Ujk_was + cos_ph * Ui(__k); } } template<typename MatrixA> void SVD<MatrixA>::ripThrough( vector_type& super_diag, const int k, const int l, const double eps ) { // Accumulate cos,sin of Ph value_type cos_ph = 0, sin_ph = 1; // 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( int i=l; i<=k; i++) { const value_type f = sin_ph * super_diag(i); super_diag(i) *= cos_ph; // Current J[l-1,l] may become unsignificant if( math::abs(f) <= eps ) break; // unnormalized sin/cos cos_ph = _M_S(i); sin_ph = -f; // sqrt(sin^2+cos^2) const value_type norm = ( _M_S(i) = hypot(cos_ph,sin_ph)); // Normalize sin/cos cos_ph /= norm; sin_ph /= norm; rotate( _M_U, i, l-1, cos_ph, sin_ph ); } } template<typename MatrixA> int 00714 SVD<MatrixA>::getWorkSubmatrix( vector_type& super_diag, const int k, const double eps ) { for( int l=k; l > 0; --l ) { if( math::abs(super_diag(l)) <= eps ) { // The breaking point: zero J[l-1,l] return l; } else if( math::abs(_M_S(l-1)) <= eps ) { // Diagonal J[l,l] turns out 0 // meaning J[l-1,l] _can_ be made // zero after some rotations ripThrough(super_diag,k,l,eps); return l; } } // Deal with J[1:k,1:k] as a whole return 0; } template<typename MatrixA> void 00738 SVD<MatrixA>::diagonalize( vector_type& super_diag, const double eps ) { for(int k=_M_cols-1; k >= 0; --k) // QR-iterate upon J[l:k,l:k] { //std::cerr << "k = " << k << "\n"; int l; // until superdiag J[k-1,k] becomes 0 while(math::abs(super_diag(k)) > eps) { l=getWorkSubmatrix(super_diag,k,eps); //std::cerr << "l = " << l << "\n"; value_type shift; // Compute a QR-shift from a bottom { // corner minor of J[l:k,l:k] order 2 value_type Jk2k1 = super_diag(k-1), // J[k-2,k-1] Jk1k = super_diag(k), Jk1k1 = _M_S(k-1), // J[k-1,k-1] Jkk = _M_S(k), Jll = _M_S(l); // J[l,l] shift = (Jk1k1-Jkk)*(Jk1k1+Jkk) + (Jk2k1-Jk1k)*(Jk2k1+Jk1k); shift /= 2*Jk1k*Jk1k1; shift += (shift < 0 ? -1 : 1) * math::sqrt(shift*shift+1); shift = ( (Jll-Jkk)*(Jll+Jkk) + Jk1k*(Jk1k1/shift-Jk1k) )/Jll; } // Carry on multiplications by T2, S2, T3... value_type cos_th = 1; value_type sin_th = 1; value_type Ji1i1 = _M_S(l); // J[i-1,i-1] at i=l+1...k for(int i=l+1; i<=k; i++) { value_type Ji1i = super_diag(i); // J[i-1,i] value_type Jii = _M_S(i); // J[i,i] sin_th *= Ji1i; Ji1i *= cos_th; cos_th = shift; value_type 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 value_type Jii1 = Jii*sin_th; // Emerged J[i,i-1] Jii *= cos_th; // new J[i,i] // Accumulate T rotations in V rotate( _M_V, i, i-1, cos_th, sin_th ); value_type cos_ph = shift, sin_ph = Jii1;// Make Si to get rid of J[i,i-1] _M_S(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 // Accumulate S rotations in U rotate( _M_U, i, i-1, cos_ph, sin_ph ); // 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; _M_S(k) = Ji1i1; } // --- end-of-QR-iterations if( _M_S(k) < 0 ) // Correct the sign of the sing number { _M_S(k) = -_M_S(k); using namespace boost::numeric::ublas; //matrix_column<matrix<value_type> > Vk ( _M_V, k); matrix_column<matrix_type> Vk ( _M_V, k); for ( unsigned l = 0; l < Vk.size(); ++ l ) { Vk( l ) = -Vk( l ); } } } } } #endif /* __SVD_HPP */

Generated by Doxygen 1.6.0 Back to index