diff options
author | Benoit Jacob <jacob.benoit.1@gmail.com> | 2009-09-02 06:36:55 -0400 |
---|---|---|
committer | Benoit Jacob <jacob.benoit.1@gmail.com> | 2009-09-02 06:36:55 -0400 |
commit | e6b77bcc6bc915ec38640ecf414726fa2ba56fba (patch) | |
tree | 35bf233dbe4df303af9c5883186e56b51f200346 /Eigen | |
parent | c16d65f01585b51f41b715a22c43983faab4299a (diff) |
JacobiSVD: implement general R-SVD using full-pivoting QR, so we now support any rectangular matrix size by reducing to the smaller of the two dimensions (which is also an optimization)
Diffstat (limited to 'Eigen')
-rw-r--r-- | Eigen/SVD | 2 | ||||
-rw-r--r-- | Eigen/src/QR/FullPivotingHouseholderQR.h | 2 | ||||
-rw-r--r-- | Eigen/src/SVD/JacobiSVD.h | 49 |
3 files changed, 40 insertions, 13 deletions
@@ -1,7 +1,7 @@ #ifndef EIGEN_SVD_MODULE_H #define EIGEN_SVD_MODULE_H -#include "Core" +#include "QR" #include "Householder" #include "Jacobi" diff --git a/Eigen/src/QR/FullPivotingHouseholderQR.h b/Eigen/src/QR/FullPivotingHouseholderQR.h index cee41b152..52405170b 100644 --- a/Eigen/src/QR/FullPivotingHouseholderQR.h +++ b/Eigen/src/QR/FullPivotingHouseholderQR.h @@ -286,7 +286,7 @@ FullPivotingHouseholderQR<MatrixType>& FullPivotingHouseholderQR<MatrixType>::co m_cols_permutation.resize(matrix.cols()); int number_of_transpositions = 0; - RealScalar biggest; + RealScalar biggest(0); for (int k = 0; k < size; ++k) { diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index b3d289706..654f8981a 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -188,15 +188,42 @@ void ei_real_2x2_jacobi_svd(const MatrixType& matrix, int p, int q, template<typename MatrixType, unsigned int Options> JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute(const MatrixType& matrix) { - MatrixType work_matrix(matrix); - int size = matrix.rows(); - if(ComputeU) m_matrixU = MatrixUType::Identity(size,size); - if(ComputeV) m_matrixV = MatrixUType::Identity(size,size); - m_singularValues.resize(size); + MatrixType work_matrix; + int rows = matrix.rows(); + int cols = matrix.cols(); + int diagSize = std::min(rows, cols); + if(ComputeU) m_matrixU = MatrixUType::Zero(rows,rows); + if(ComputeV) m_matrixV = MatrixVType::Zero(cols,cols); + m_singularValues.resize(diagSize); const RealScalar precision = 2 * epsilon<Scalar>(); + if(rows > cols) + { + FullPivotingHouseholderQR<MatrixType> qr(matrix); + work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView<UpperTriangular>(); + if(ComputeU) m_matrixU = qr.matrixQ(); + if(ComputeV) + for(int i = 0; i < cols; i++) + m_matrixV.coeffRef(qr.colsPermutation().coeff(i),i) = Scalar(1); + } + else if(rows < cols) + { + FullPivotingHouseholderQR<MatrixType> qr(MatrixType(matrix.adjoint())); + work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView<UpperTriangular>().adjoint(); + if(ComputeV) m_matrixV = qr.matrixQ(); + if(ComputeU) + for(int i = 0; i < rows; i++) + m_matrixU.coeffRef(qr.colsPermutation().coeff(i),i) = Scalar(1); + + } + else + { + work_matrix = matrix; + if(ComputeU) m_matrixU.diagonal().setOnes(); + if(ComputeV) m_matrixV.diagonal().setOnes(); + } sweep_again: - for(int p = 1; p < size; ++p) + for(int p = 1; p < diagSize; ++p) { for(int q = 0; q < p; ++q) { @@ -219,27 +246,27 @@ sweep_again: RealScalar biggestOnDiag = work_matrix.diagonal().cwise().abs().maxCoeff(); RealScalar maxAllowedOffDiag = biggestOnDiag * precision; - for(int p = 0; p < size; ++p) + for(int p = 0; p < diagSize; ++p) { for(int q = 0; q < p; ++q) if(ei_abs(work_matrix.coeff(p,q)) > maxAllowedOffDiag) goto sweep_again; - for(int q = p+1; q < size; ++q) + for(int q = p+1; q < diagSize; ++q) if(ei_abs(work_matrix.coeff(p,q)) > maxAllowedOffDiag) goto sweep_again; } - for(int i = 0; i < size; ++i) + for(int i = 0; i < diagSize; ++i) { RealScalar a = ei_abs(work_matrix.coeff(i,i)); m_singularValues.coeffRef(i) = a; if(ComputeU && (a!=RealScalar(0))) m_matrixU.col(i) *= work_matrix.coeff(i,i)/a; } - for(int i = 0; i < size; i++) + for(int i = 0; i < diagSize; i++) { int pos; - m_singularValues.end(size-i).maxCoeff(&pos); + m_singularValues.end(diagSize-i).maxCoeff(&pos); if(pos) { pos += i; |