36 #include <boost/numeric/ublas/matrix.hpp>
37 #include <boost/numeric/ublas/vector.hpp>
38 #include <boost/numeric/ublas/io.hpp>
40 #include <feel/feelcore/feel.hpp>
41 #include <feel/feelcore/traits.hpp>
46 namespace ublas = boost::numeric::ublas;
72 template<
typename MatrixA>
82 typedef typename MatrixA::value_type value_type;
83 typedef MatrixA matrix_type;
84 typedef boost::numeric::ublas::vector<value_type> vector_type;
92 SVD( matrix_type
const& __A );
109 matrix_type
const&
U()
const
115 matrix_type
const&
V()
const
121 vector_type
const&
S()
const
134 return max( M_S )/min( M_S );
139 __max = *std::max_element( M_S.begin(), M_S.end() );
140 __min = *std::min_element( M_S.begin(), M_S.end() );
233 value_type leftHouseholder( MatrixA& A,
const int i );
297 value_type rightHouseholder( MatrixA& A,
const int i );
309 value_type bidiagonalize( vector_type& super_diag,
const matrix_type& _A );
368 void rotate( matrix_type&
U,
const int i,
const int j,
const value_type cos_ph,
const value_type sin_ph );
383 void ripThrough( vector_type& super_diag,
const int k,
const int l,
const double eps );
393 int getWorkSubmatrix( vector_type& super_diag,
const int k,
const double eps );
396 void diagonalize( vector_type& super_diag,
const double eps );
409 template<
typename MatrixA>
418 typedef typename MatrixA::value_type value_type;
419 typedef MatrixA matrix_type;
420 typedef boost::numeric::ublas::vector<value_type> vector_type;
428 SOrth( matrix_type
const& __A );
445 matrix_type
const& Q()
const
475 template<
typename MatrixA>
476 SOrth<MatrixA>::SOrth( MatrixA
const& A )
482 using namespace boost::numeric::ublas;
484 SVD<matrix_type> __svd( A );
490 else if ( M_rows == 1 )
491 __S( 0 ) = __svd.S()( 0 );
493 const value_type eps = std::numeric_limits<value_type>::epsilon() * *std::max_element( __S.begin(), __S.end() )* std::max( M_rows, M_cols );
509 template<
typename MatrixA>
510 SVD<MatrixA>::SVD( MatrixA
const& A )
514 M_U( M_rows, M_rows ),
515 M_V( M_cols, M_cols ),
522 bool do_transpose = M_rows < M_cols;
529 M_U.resize( M_rows, M_rows );
530 M_V.resize( M_cols, M_cols );
531 M_S.resize( M_cols );
532 B = ublas::trans( A );
535 M_U = boost::numeric::ublas::identity_matrix<value_type>( M_rows );
536 M_V = boost::numeric::ublas::identity_matrix<value_type>( M_cols );
538 vector_type super_diag( M_cols );
540 const value_type bidiag_norm = bidiagonalize( super_diag,B );
543 const value_type eps = std::numeric_limits<value_type>::epsilon() * bidiag_norm;
545 diagonalize( super_diag,eps );
549 matrix_type Temp( M_V );
550 M_V = ublas::trans( M_U );
551 M_U = ublas::trans( Temp );
560 template<
typename MatrixA>
561 typename SVD<MatrixA>::value_type
562 SVD<MatrixA>::leftHouseholder( matrix_type& A,
const int i )
564 value_type scale = 0;
566 for (
int k=i; k< M_rows; k++ )
567 scale += math::abs( A( k,i ) );
572 value_type Norm_sqr = 0;
574 for (
int k=i; k< M_rows; k++ )
577 Norm_sqr += A( k,i )*A( k,i );
580 value_type new_Aii = math::sqrt( Norm_sqr );
582 if ( A( i,i ) > 0 ) new_Aii = -new_Aii;
584 const value_type beta = - A( i,i )*new_Aii + Norm_sqr;
587 for (
int j=i+1; j<M_cols; j++ )
589 value_type factor = 0;
591 for (
int k=i; k< M_rows; k++ )
592 factor += A( k,i ) * A( k,j );
596 for (
int k=i; k< M_rows; k++ )
597 A( k,j ) -= A( k,i ) * factor;
602 value_type factor = 0;
605 factor += A( k,i ) * M_U( j,k );
610 M_U( j,k ) -= A( k,i ) * factor;
613 return new_Aii * scale;
616 template<
typename MatrixA>
617 typename SVD<MatrixA>::value_type
618 SVD<MatrixA>::rightHouseholder( matrix_type& A,
const int i )
620 value_type scale = 0;
622 for ( uint16_type k=i+1; k<M_cols; k++ )
623 scale += math::abs( A( i,k ) );
628 value_type Norm_sqr = 0;
630 for ( uint16_type k=i+1; k<M_cols; k++ )
633 Norm_sqr += A( i,k )*A( i,k );
636 value_type new_Aii1 = math::sqrt( Norm_sqr );
638 if ( A( i,i+1 ) > 0 )
639 new_Aii1 = -new_Aii1;
641 const value_type beta = - A( i,i+1 )*new_Aii1 + Norm_sqr;
642 A( i,i+1 ) -= new_Aii1;
644 for ( uint16_type j=i+1; j<M_rows; j++ )
646 value_type factor = 0;
648 for ( uint16_type k=i+1; k<M_cols; k++ )
649 factor += A( i,k ) * A( j,k );
653 for ( uint16_type k=i+1; k<M_cols; k++ )
654 A( j,k ) -= A( i,k ) * factor;
657 for ( uint16_type j=0; j<M_cols; j++ )
659 value_type factor = 0;
661 for ( uint16_type k=i+1; k<M_cols; k++ )
662 factor += A( i,k ) * M_V( j,k );
666 for ( uint16_type k=i+1; k<M_cols; k++ )
667 M_V( j,k ) -= A( i,k ) * factor;
670 return new_Aii1 * scale;
673 template<
typename MatrixA>
674 typename SVD<MatrixA>::value_type
675 SVD<MatrixA>::bidiagonalize( vector_type& super_diag,
const matrix_type& _A )
677 value_type __norm_acc = 0;
685 for ( uint16_type __i = 0; __i < M_cols; ++__i )
687 M_S( __i ) = leftHouseholder( A,__i );
689 const value_type& diagi = M_S( __i );
691 if ( __i < M_cols-1 )
692 super_diag( __i+1 ) = rightHouseholder( A,__i );
694 __norm_acc = std::max( __norm_acc, math::abs( diagi ) + math::abs( super_diag( __i ) ) );
701 template<
typename MatrixA>
703 SVD<MatrixA>::rotate( matrix_type& U,
const int i,
const int j,
704 const value_type cos_ph,
const value_type sin_ph )
706 using namespace boost::numeric::ublas;
708 matrix_column<matrix_type> Ui ( U, i );
709 matrix_column<matrix_type> Uj ( U, j );
711 for (
unsigned __k = 0; __k < Ui.size(); ++ __k )
713 const value_type Ujk_was = Uj( __k );
715 Uj( __k ) = cos_ph * Ujk_was + sin_ph * Ui( __k );
716 Ui( __k ) = -sin_ph * Ujk_was + cos_ph * Ui( __k );
721 template<
typename MatrixA>
723 SVD<MatrixA>::ripThrough( vector_type& super_diag,
const int k,
const int l,
const double eps )
726 value_type cos_ph = 0, sin_ph = 1;
736 for (
int i=l; i<=k; i++ )
738 const value_type f = sin_ph * super_diag( i );
740 super_diag( i ) *= cos_ph;
743 if ( math::abs( f ) <= eps )
751 const value_type norm = ( M_S( i ) = hypot( cos_ph,sin_ph ) );
757 rotate( M_U, i, l-1, cos_ph, sin_ph );
760 template<
typename MatrixA>
762 SVD<MatrixA>::getWorkSubmatrix( vector_type& super_diag,
const int k,
const double eps )
764 for (
int l=k; l > 0; --l )
766 if ( math::abs( super_diag( l ) ) <= eps )
772 else if ( math::abs( M_S( l-1 ) ) <= eps )
777 ripThrough( super_diag,k,l,eps );
786 template<
typename MatrixA>
788 SVD<MatrixA>::diagonalize( vector_type& super_diag,
const double eps )
790 for (
int k=M_cols-1; k >= 0; --k )
796 while ( math::abs( super_diag( k ) ) > eps )
798 l=getWorkSubmatrix( super_diag,k,eps );
803 value_type Jk2k1 = super_diag( k-1 ),
804 Jk1k = super_diag( k ),
808 shift = ( Jk1k1-Jkk )*( Jk1k1+Jkk ) + ( Jk2k1-Jk1k )*( Jk2k1+Jk1k );
809 shift /= 2*Jk1k*Jk1k1;
810 shift += ( shift < 0 ? -1 : 1 ) * math::sqrt( shift*shift+1 );
811 shift = ( ( Jll-Jkk )*( Jll+Jkk ) + Jk1k*( Jk1k1/shift-Jk1k ) )/Jll;
814 value_type cos_th = 1;
815 value_type sin_th = 1;
816 value_type Ji1i1 = M_S( l );
818 for (
int i=l+1; i<=k; i++ )
820 value_type Ji1i = super_diag( i );
821 value_type Jii = M_S( i );
827 value_type norm_f = ( super_diag( i-1 ) = hypot( cos_th,sin_th ) );
828 cos_th /= norm_f, sin_th /= norm_f;
831 shift = cos_th*Ji1i1 + sin_th*Ji1i;
832 Ji1i = -sin_th*Ji1i1 + cos_th*Ji1i;
833 const value_type Jii1 = Jii*sin_th;
837 rotate( M_V, i, i-1, cos_th, sin_th );
839 value_type cos_ph = shift, sin_ph = Jii1;
840 M_S( i-1 ) = ( norm_f = hypot( cos_ph,sin_ph ) );
855 shift = cos_ph * Ji1i + sin_ph*Jii;
856 Ji1i1 = -sin_ph*Ji1i + cos_ph*Jii;
859 rotate( M_U, i, i-1, cos_ph, sin_ph );
862 cos_th = cos_ph, sin_th = sin_ph;
868 super_diag( k ) = shift;
874 M_S( k ) = -M_S( k );
876 using namespace boost::numeric::ublas;
878 matrix_column<matrix_type> Vk ( M_V, k );
880 for (
unsigned l = 0; l < Vk.size(); ++ l )
matrix_type const & U() const
Definition: svd.hpp:109
vector_type const & S() const
Definition: svd.hpp:121
Singular Value Decomposition of a rectangular matrix.
Definition: svd.hpp:73
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
matrix_type const & V() const
Definition: svd.hpp:115
value_type conditionNumber() const
Definition: svd.hpp:132