Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
svd.hpp
Go to the documentation of this file.
1 /* -*- 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
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2005-03-27
7 
8  Copyright (C) 2005,2006 EPFL
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
29 #ifndef __SVD_HPP
30 #define __SVD_HPP 1
31 
32 #include <cmath>
33 #include <limits>
34 #include <algorithm>
35 
36 #include <boost/numeric/ublas/matrix.hpp>
37 #include <boost/numeric/ublas/vector.hpp>
38 #include <boost/numeric/ublas/io.hpp>
39 
40 #include <feel/feelcore/feel.hpp>
41 #include <feel/feelcore/traits.hpp>
42 
43 
44 namespace Feel
45 {
46 namespace ublas = boost::numeric::ublas;
72 template<typename MatrixA>
73 class SVD
74 {
75 public:
76 
77 
81 
82  typedef typename MatrixA::value_type value_type;
83  typedef MatrixA matrix_type;
84  typedef boost::numeric::ublas::vector<value_type> vector_type;
85 
87 
91 
92  SVD( matrix_type const& __A );
93  ~SVD() {}
94 
96 
100 
101 
103 
107 
109  matrix_type const& U() const
110  {
111  return M_U;
112  }
113 
115  matrix_type const& V() const
116  {
117  return M_V;
118  }
119 
121  vector_type const& S() const
122  {
123  return M_S;
124  }
125 
132  value_type conditionNumber() const
133  {
134  return max( M_S )/min( M_S );
135  }
136 
137  void conditionNumber( value_type& __max, value_type& __min ) const
138  {
139  __max = *std::max_element( M_S.begin(), M_S.end() );
140  __min = *std::min_element( M_S.begin(), M_S.end() );
141  }
142 
144 
148 
149 
151 
155 
156 
158 
159 protected:
160 
161 private:
162 
163  SVD( SVD const & );
164 
233  value_type leftHouseholder( MatrixA& A, const int i );
234 
235 
297  value_type rightHouseholder( MatrixA& A, const int i );
298 
309  value_type bidiagonalize( vector_type& super_diag, const matrix_type& _A );
310 
311  /*
312  \brief QR-diagonalization of a bidiagonal matrix
313 
314  After bidiagonalization we get a bidiagonal matrix J:
315  (1) J = U' * A * V
316  The present method turns J into a matrix JJ by applying a set of
317  orthogonal transforms
318  (2) JJ = S' * J * T
319  Orthogonal matrices S and T are chosen so that JJ were also a
320  bidiagonal matrix, but with superdiag elements smaller than those of J.
321  We repeat (2) until non-diag elements of JJ become smaller than EPS
322  and can be disregarded.
323  Matrices S and T are constructed as
324  (3) S = S1 * S2 * S3 ... Sn, and similarly T
325  where Sk and Tk are matrices of simple rotations
326  (4) Sk[i,j] = i==j ? 1 : 0 for all i>k or i<k-1
327  Sk[k-1,k-1] = cos(Phk), Sk[k-1,k] = -sin(Phk),
328  SK[k,k-1] = sin(Phk), Sk[k,k] = cos(Phk), k=2..N
329  Matrix Tk is constructed similarly, only with angle Thk rather than Phk.
330 
331  Thus left multiplication of J by SK' can be spelled out as
332  (5) (Sk' * J)[i,j] = J[i,j] when i>k or i<k-1,
333  [k-1,j] = cos(Phk)*J[k-1,j] + sin(Phk)*J[k,j]
334  [k,j] = -sin(Phk)*J[k-1,j] + cos(Phk)*J[k,j]
335  That is, k-1 and k rows of J are replaced by their linear combinations;
336  the rest of J is unaffected. Right multiplication of J by Tk similarly
337  changes only k-1 and k columns of J.
338  Matrix T2 is chosen the way that T2'J'JT2 were a QR-transform with a
339  shift. Note that multiplying J by T2 gives rise to a J[2,1] element of
340  the product J (which is below the main diagonal). Angle Ph2 is then
341  chosen so that multiplication by S2' (which combines 1 and 2 rows of J)
342  gets rid of that elemnent. But this will create a [1,3] non-zero element.
343  T3 is made to make it disappear, but this leads to [3,2], etc.
344  In the end, Sn removes a [N,N-1] element of J and matrix S'JT becomes
345  bidiagonal again. However, because of a special choice
346  of T2 (QR-algorithm), its non-diag elements are smaller than those of J.
347 
348  All this process in more detail is described in
349  J.H. Wilkinson, C. Reinsch. Linear algebra - Springer-Verlag, 1971
350 
351  If during transforms (1), JJ[N-1,N] turns 0, then JJ[N,N] is a singular
352  number (possibly with a wrong (that is, negative) sign). This is a
353  consequence of Frantsis' Theorem, see the book above. In that case, we can
354  eliminate the N-th row and column of JJ and carry out further transforms
355  with a smaller matrix. If any other superdiag element of JJ turns 0,
356  then JJ effectively falls into two independent matrices. We will process
357  them independently (the bottom one first).
358 
359  Since matrix J is a bidiagonal, it can be stored efficiently. As a matter
360  of fact, its N diagonal elements are in array Sig, and superdiag elements
361  are stored in array super_diag.
362 
363  Carry out U * S with a rotation matrix
364  S (which combines i-th and j-th columns
365  of U, i>j)
366 
367  */
368  void rotate( matrix_type& U, const int i, const int j, const value_type cos_ph, const value_type sin_ph );
369 
370  /*
371  A diagonal element J[l-1,l-1] turns out 0 at the k-th step of the
372  algorithm. That means that one of the original matrix' singular numbers
373  shall be zero. In that case, we multiply J by specially selected
374  matrices S' of simple rotations to eliminate a superdiag element J[l-1,l].
375  After that, matrix J falls into two pieces, which can be dealt with
376  in a regular way (the bottom piece first).
377 
378  These special S transformations are accumulated into matrix U: since J
379  is left-multiplied by S', U would be right-multiplied by S. Transform
380  formulas for doing these rotations are similar to (5) above. See the
381  book cited above for more details.
382  */
383  void ripThrough( vector_type& super_diag, const int k, const int l, const double eps );
384 
393  int getWorkSubmatrix( vector_type& super_diag, const int k, const double eps );
394 
396  void diagonalize( vector_type& super_diag, const double eps );
397 
398 private:
399 
400  uint M_rows;
401  uint M_cols;
402 
403  matrix_type M_U;
404  matrix_type M_V;
405  vector_type M_S;
406 
407 };
408 
409 template<typename MatrixA>
410 class SOrth
411 {
412 public:
413 
417 
418  typedef typename MatrixA::value_type value_type;
419  typedef MatrixA matrix_type;
420  typedef boost::numeric::ublas::vector<value_type> vector_type;
421 
423 
427 
428  SOrth( matrix_type const& __A );
429  ~SOrth() {}
430 
432 
436 
437 
439 
443 
445  matrix_type const& Q() const
446  {
447  return M_Q;
448  }
449 
451 
455 
456 
458 
462 
463 
465 private:
466 
467  uint M_rows;
468  uint M_cols;
469  matrix_type M_Q;
470 };
471 
475 template<typename MatrixA>
476 SOrth<MatrixA>::SOrth( MatrixA const& A )
477  :
478  M_rows( A.size1() ),
479  M_cols( A.size2() ),
480  M_Q()
481 {
482  using namespace boost::numeric::ublas;
483 
484  SVD<matrix_type> __svd( A );
485  vector_type __S;
486 
487  if ( M_rows > 1 )
488  __S = __svd.S();
489 
490  else if ( M_rows == 1 )
491  __S( 0 ) = __svd.S()( 0 );
492 
493  const value_type eps = std::numeric_limits<value_type>::epsilon() * *std::max_element( __S.begin(), __S.end() )* std::max( M_rows, M_cols );
494 
495  /*[m,n] = size(A);
496  if m > 1, s = diag(S);
497  elseif m == 1, s = S(1);
498  else s = 0;
499  end
500  tol = max(m,n) * max(s) * eps;
501  r = sum(s > tol);
502  Q = U(:,1:r);*/
503 
504 }
505 /*
506  * SVD
507  */
508 
509 template<typename MatrixA>
510 SVD<MatrixA>::SVD( MatrixA const& A )
511  :
512  M_rows( A.size1() ),
513  M_cols( A.size2() ),
514  M_U( M_rows, M_rows ),
515  M_V( M_cols, M_cols ),
516  M_S( M_cols )
517 {
518  //FEELPP_ASSERT( M_rows >= M_cols ).error("The Matrix A should have at least as many rows as it has columns");
519 
520  MatrixA B( A );
521 
522  bool do_transpose = M_rows < M_cols;
523 
524  if ( do_transpose )
525  {
526  // transpose
527  M_rows = A.size2();
528  M_cols = A.size1();
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 );
533  }
534 
535  M_U = boost::numeric::ublas::identity_matrix<value_type>( M_rows );
536  M_V = boost::numeric::ublas::identity_matrix<value_type>( M_cols );
537 
538  vector_type super_diag( M_cols );
539 
540  const value_type bidiag_norm = bidiagonalize( super_diag,B );
541 
542  // define Significance threshold
543  const value_type eps = std::numeric_limits<value_type>::epsilon() * bidiag_norm;
544 
545  diagonalize( super_diag,eps );
546 
547  if ( do_transpose )
548  {
549  matrix_type Temp( M_V );
550  M_V = ublas::trans( M_U );
551  M_U = ublas::trans( Temp );
552  }
553 }
554 
555 /*
556  Bidiagonalization
557 */
558 
559 
560 template<typename MatrixA>
561 typename SVD<MatrixA>::value_type
562 SVD<MatrixA>::leftHouseholder( matrix_type& A, const int i )
563 {
564  value_type scale = 0; // Compute the scaling factor
565 
566  for ( int k=i; k< M_rows; k++ )
567  scale += math::abs( A( k,i ) );
568 
569  if ( scale == 0 ) // If A[,i] is a null vector, no
570  return 0; // transform is required
571 
572  value_type Norm_sqr = 0; // Scale UPi (that is, A[,i])
573 
574  for ( int k=i; k< M_rows; k++ ) // and compute its norm, Norm^2
575  {
576  A( k,i ) /= scale;
577  Norm_sqr += A( k,i )*A( k,i );
578  }
579 
580  value_type new_Aii = math::sqrt( Norm_sqr ); // new_Aii = -Norm, Norm has the
581 
582  if ( A( i,i ) > 0 ) new_Aii = -new_Aii; // same sign as Aii (that is, UPi[i])
583 
584  const value_type beta = - A( i,i )*new_Aii + Norm_sqr;
585  A( i,i ) -= new_Aii; // UPi[i] = A[i,i] - (-Norm)
586 
587  for ( int j=i+1; j<M_cols; j++ ) // Transform i+1:N columns of A
588  {
589  value_type factor = 0;
590 
591  for ( int k=i; k< M_rows; k++ )
592  factor += A( k,i ) * A( k,j );
593 
594  factor /= beta;
595 
596  for ( int k=i; k< M_rows; k++ )
597  A( k,j ) -= A( k,i ) * factor;
598  }
599 
600  for ( size_type j=0; j< M_rows; j++ ) // Accumulate the transform in U
601  {
602  value_type factor = 0;
603 
604  for ( size_type k=i; k< M_rows; k++ )
605  factor += A( k,i ) * M_U( j,k ); // Compute U[j,] * UPi
606 
607  factor /= beta;
608 
609  for ( size_type k=i; k< M_rows; k++ )
610  M_U( j,k ) -= A( k,i ) * factor;
611  }
612 
613  return new_Aii * scale; // Scale new Aii back (our new Sig[i])
614 }
615 
616 template<typename MatrixA>
617 typename SVD<MatrixA>::value_type
618 SVD<MatrixA>::rightHouseholder( matrix_type& A, const int i )
619 {
620  value_type scale = 0; // Compute the scaling factor
621 
622  for ( uint16_type k=i+1; k<M_cols; k++ )
623  scale += math::abs( A( i,k ) );
624 
625  if ( scale == 0 ) // If A[i,] is a null vector, no
626  return 0; // transform is required
627 
628  value_type Norm_sqr = 0; // Scale VPi (that is, A[i,])
629 
630  for ( uint16_type k=i+1; k<M_cols; k++ ) // and compute its norm, Norm^2
631  {
632  A( i,k ) /= scale;
633  Norm_sqr += A( i,k )*A( i,k );
634  }
635 
636  value_type new_Aii1 = math::sqrt( Norm_sqr ); // new_Aii1 = -Norm, Norm has the
637 
638  if ( A( i,i+1 ) > 0 ) // same sign as
639  new_Aii1 = -new_Aii1; // Aii1 (that is, VPi[i+1])
640 
641  const value_type beta = - A( i,i+1 )*new_Aii1 + Norm_sqr;
642  A( i,i+1 ) -= new_Aii1; // VPi[i+1] = A[i,i+1] - (-Norm)
643 
644  for ( uint16_type j=i+1; j<M_rows; j++ ) // Transform i+1:M rows of A
645  {
646  value_type factor = 0;
647 
648  for ( uint16_type k=i+1; k<M_cols; k++ )
649  factor += A( i,k ) * A( j,k ); // Compute A[j,] * VPi
650 
651  factor /= beta;
652 
653  for ( uint16_type k=i+1; k<M_cols; k++ )
654  A( j,k ) -= A( i,k ) * factor;
655  }
656 
657  for ( uint16_type j=0; j<M_cols; j++ ) // Accumulate the transform in V
658  {
659  value_type factor = 0;
660 
661  for ( uint16_type k=i+1; k<M_cols; k++ )
662  factor += A( i,k ) * M_V( j,k ); // Compute V[j,] * VPi
663 
664  factor /= beta;
665 
666  for ( uint16_type k=i+1; k<M_cols; k++ )
667  M_V( j,k ) -= A( i,k ) * factor;
668  }
669 
670  return new_Aii1 * scale; // Scale new Aii1 back
671 }
672 
673 template<typename MatrixA>
674 typename SVD<MatrixA>::value_type
675 SVD<MatrixA>::bidiagonalize( vector_type& super_diag, const matrix_type& _A )
676 {
677  value_type __norm_acc = 0;
678 
679  // No superdiag elem above A(0,0)
680  super_diag( 0 ) = 0;
681 
682  // A being transformed
683  matrix_type A = _A;
684 
685  for ( uint16_type __i = 0; __i < M_cols; ++__i )
686  {
687  M_S( __i ) = leftHouseholder( A,__i );
688  //std::cout << "S=" << M_S << "\n";
689  const value_type& diagi = M_S( __i );
690 
691  if ( __i < M_cols-1 )
692  super_diag( __i+1 ) = rightHouseholder( A,__i );
693 
694  __norm_acc = std::max( __norm_acc, math::abs( diagi ) + math::abs( super_diag( __i ) ) );
695  }
696 
697  return __norm_acc;
698 }
699 
700 
701 template<typename MatrixA>
702 void
703 SVD<MatrixA>::rotate( matrix_type& U, const int i, const int j,
704  const value_type cos_ph, const value_type sin_ph )
705 {
706  using namespace boost::numeric::ublas;
707 
708  matrix_column<matrix_type> Ui ( U, i );
709  matrix_column<matrix_type> Uj ( U, j );
710 
711  for ( unsigned __k = 0; __k < Ui.size(); ++ __k )
712  {
713  const value_type Ujk_was = Uj( __k );
714 
715  Uj( __k ) = cos_ph * Ujk_was + sin_ph * Ui( __k );
716  Ui( __k ) = -sin_ph * Ujk_was + cos_ph * Ui( __k );
717  }
718 }
719 
720 
721 template<typename MatrixA>
722 void
723 SVD<MatrixA>::ripThrough( vector_type& super_diag, const int k, const int l, const double eps )
724 {
725  // Accumulate cos,sin of Ph
726  value_type cos_ph = 0, sin_ph = 1;
727 
728  // The first step of the loop below
729  // (when i==l) would eliminate J[l-1,l],
730  // which is stored in super_diag(l)
731  // However, it gives rise to J[l-1,l+1]
732  // and J[l,l+2]
733  // The following steps eliminate these
734  // until they fall below
735  // significance
736  for ( int i=l; i<=k; i++ )
737  {
738  const value_type f = sin_ph * super_diag( i );
739 
740  super_diag( i ) *= cos_ph;
741 
742  // Current J[l-1,l] may become unsignificant
743  if ( math::abs( f ) <= eps )
744  break;
745 
746  // unnormalized sin/cos
747  cos_ph = M_S( i );
748  sin_ph = -f;
749 
750  // sqrt(sin^2+cos^2)
751  const value_type norm = ( M_S( i ) = hypot( cos_ph,sin_ph ) );
752 
753  // Normalize sin/cos
754  cos_ph /= norm;
755  sin_ph /= norm;
756 
757  rotate( M_U, i, l-1, cos_ph, sin_ph );
758  }
759 }
760 template<typename MatrixA>
761 int
762 SVD<MatrixA>::getWorkSubmatrix( vector_type& super_diag, const int k, const double eps )
763 {
764  for ( int l=k; l > 0; --l )
765  {
766  if ( math::abs( super_diag( l ) ) <= eps )
767  {
768  // The breaking point: zero J[l-1,l]
769  return l;
770  }
771 
772  else if ( math::abs( M_S( l-1 ) ) <= eps )
773  {
774  // Diagonal J[l,l] turns out 0
775  // meaning J[l-1,l] _can_ be made
776  // zero after some rotations
777  ripThrough( super_diag,k,l,eps );
778  return l;
779  }
780  }
781 
782  // Deal with J[1:k,1:k] as a whole
783  return 0;
784 }
785 
786 template<typename MatrixA>
787 void
788 SVD<MatrixA>::diagonalize( vector_type& super_diag, const double eps )
789 {
790  for ( int k=M_cols-1; k >= 0; --k ) // QR-iterate upon J[l:k,l:k]
791  {
792  //std::cerr << "k = " << k << "\n";
793  int l;
794 
795  // until superdiag J[k-1,k] becomes 0
796  while ( math::abs( super_diag( k ) ) > eps )
797  {
798  l=getWorkSubmatrix( super_diag,k,eps );
799  //std::cerr << "l = " << l << "\n";
800  value_type shift; // Compute a QR-shift from a bottom
801  {
802  // corner minor of J[l:k,l:k] order 2
803  value_type Jk2k1 = super_diag( k-1 ), // J[k-2,k-1]
804  Jk1k = super_diag( k ),
805  Jk1k1 = M_S( k-1 ), // J[k-1,k-1]
806  Jkk = M_S( k ),
807  Jll = M_S( l ); // J[l,l]
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;
812  }
813  // Carry on multiplications by T2, S2, T3...
814  value_type cos_th = 1;
815  value_type sin_th = 1;
816  value_type Ji1i1 = M_S( l ); // J[i-1,i-1] at i=l+1...k
817 
818  for ( int i=l+1; i<=k; i++ )
819  {
820  value_type Ji1i = super_diag( i ); // J[i-1,i]
821  value_type Jii = M_S( i ); // J[i,i]
822 
823  sin_th *= Ji1i;
824  Ji1i *= cos_th;
825  cos_th = shift;
826 
827  value_type norm_f = ( super_diag( i-1 ) = hypot( cos_th,sin_th ) );
828  cos_th /= norm_f, sin_th /= norm_f;
829 
830  // Rotate J[i-1:i,i-1:i] by Ti
831  shift = cos_th*Ji1i1 + sin_th*Ji1i; // new J[i-1,i-1]
832  Ji1i = -sin_th*Ji1i1 + cos_th*Ji1i; // J[i-1,i] after rotation
833  const value_type Jii1 = Jii*sin_th; // Emerged J[i,i-1]
834  Jii *= cos_th; // new J[i,i]
835 
836  // Accumulate T rotations in V
837  rotate( M_V, i, i-1, cos_th, sin_th );
838 
839  value_type cos_ph = shift, sin_ph = Jii1;// Make Si to get rid of J[i,i-1]
840  M_S( i-1 ) = ( norm_f = hypot( cos_ph,sin_ph ) ); // New J[i-1,i-1]
841 
842  if ( norm_f == 0 ) // If norm =0, rotation angle
843  {
844  cos_ph = cos_th;
845  sin_ph = sin_th; // can be anything now
846  }
847 
848  else
849  {
850  cos_ph /= norm_f;
851  sin_ph /= norm_f;
852  }
853 
854  // Rotate J[i-1:i,i-1:i] by Si
855  shift = cos_ph * Ji1i + sin_ph*Jii; // New J[i-1,i]
856  Ji1i1 = -sin_ph*Ji1i + cos_ph*Jii; // New Jii, would carry over as J[i-1,i-1] for next i
857 
858  // Accumulate S rotations in U
859  rotate( M_U, i, i-1, cos_ph, sin_ph );
860 
861  // Jii1 disappears, sin_th would
862  cos_th = cos_ph, sin_th = sin_ph; // carry over a (scaled) J[i-1,i+1]
863  // to eliminate on the next i, cos_th
864  // would carry over a scaled J[i,i+1]
865  }
866 
867  super_diag( l ) = 0; // Supposed to be eliminated by now
868  super_diag( k ) = shift;
869  M_S( k ) = Ji1i1;
870  } // --- end-of-QR-iterations
871 
872  if ( M_S( k ) < 0 ) // Correct the sign of the sing number
873  {
874  M_S( k ) = -M_S( k );
875 
876  using namespace boost::numeric::ublas;
877  //matrix_column<matrix<value_type> > Vk ( M_V, k);
878  matrix_column<matrix_type> Vk ( M_V, k );
879 
880  for ( unsigned l = 0; l < Vk.size(); ++ l )
881  {
882  Vk( l ) = -Vk( l );
883  }
884  }
885  }
886 }
887 
888 
889 }
890 #endif /* __SVD_HPP */
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

Generated on Sun Dec 22 2013 13:11:13 for Feel++ by doxygen 1.8.5