Jaesung Oh / Eigen

Eigne Matrix Class Library

Dependents:   MPC_current_control HydraulicControlBoard_SW AHRS Test_ekf ... more

src/SVD/JacobiSVD.h@0:13a5d365ba16, 2016-10-13 (annotated)

Committer:
ykuroda
Date:
Thu Oct 13 04:07:23 2016 +0000
Revision:
0:13a5d365ba16
First commint, Eigne Matrix Class Library

Who changed what in which revision?

UserRevisionLine numberNew contents of line
ykuroda 0:13a5d365ba16 1 // This file is part of Eigen, a lightweight C++ template library
ykuroda 0:13a5d365ba16 2 // for linear algebra.
ykuroda 0:13a5d365ba16 3 //
ykuroda 0:13a5d365ba16 4 // Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
ykuroda 0:13a5d365ba16 5 //
ykuroda 0:13a5d365ba16 6 // This Source Code Form is subject to the terms of the Mozilla
ykuroda 0:13a5d365ba16 7 // Public License v. 2.0. If a copy of the MPL was not distributed
ykuroda 0:13a5d365ba16 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
ykuroda 0:13a5d365ba16 9
ykuroda 0:13a5d365ba16 10 #ifndef EIGEN_JACOBISVD_H
ykuroda 0:13a5d365ba16 11 #define EIGEN_JACOBISVD_H
ykuroda 0:13a5d365ba16 12
ykuroda 0:13a5d365ba16 13 namespace Eigen {
ykuroda 0:13a5d365ba16 14
ykuroda 0:13a5d365ba16 15 namespace internal {
ykuroda 0:13a5d365ba16 16 // forward declaration (needed by ICC)
ykuroda 0:13a5d365ba16 17 // the empty body is required by MSVC
ykuroda 0:13a5d365ba16 18 template<typename MatrixType, int QRPreconditioner,
ykuroda 0:13a5d365ba16 19 bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
ykuroda 0:13a5d365ba16 20 struct svd_precondition_2x2_block_to_be_real {};
ykuroda 0:13a5d365ba16 21
ykuroda 0:13a5d365ba16 22 /*** QR preconditioners (R-SVD)
ykuroda 0:13a5d365ba16 23 ***
ykuroda 0:13a5d365ba16 24 *** Their role is to reduce the problem of computing the SVD to the case of a square matrix.
ykuroda 0:13a5d365ba16 25 *** This approach, known as R-SVD, is an optimization for rectangular-enough matrices, and is a requirement for
ykuroda 0:13a5d365ba16 26 *** JacobiSVD which by itself is only able to work on square matrices.
ykuroda 0:13a5d365ba16 27 ***/
ykuroda 0:13a5d365ba16 28
ykuroda 0:13a5d365ba16 29 enum { PreconditionIfMoreColsThanRows, PreconditionIfMoreRowsThanCols };
ykuroda 0:13a5d365ba16 30
ykuroda 0:13a5d365ba16 31 template<typename MatrixType, int QRPreconditioner, int Case>
ykuroda 0:13a5d365ba16 32 struct qr_preconditioner_should_do_anything
ykuroda 0:13a5d365ba16 33 {
ykuroda 0:13a5d365ba16 34 enum { a = MatrixType::RowsAtCompileTime != Dynamic &&
ykuroda 0:13a5d365ba16 35 MatrixType::ColsAtCompileTime != Dynamic &&
ykuroda 0:13a5d365ba16 36 MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime,
ykuroda 0:13a5d365ba16 37 b = MatrixType::RowsAtCompileTime != Dynamic &&
ykuroda 0:13a5d365ba16 38 MatrixType::ColsAtCompileTime != Dynamic &&
ykuroda 0:13a5d365ba16 39 MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime,
ykuroda 0:13a5d365ba16 40 ret = !( (QRPreconditioner == NoQRPreconditioner) ||
ykuroda 0:13a5d365ba16 41 (Case == PreconditionIfMoreColsThanRows && bool(a)) ||
ykuroda 0:13a5d365ba16 42 (Case == PreconditionIfMoreRowsThanCols && bool(b)) )
ykuroda 0:13a5d365ba16 43 };
ykuroda 0:13a5d365ba16 44 };
ykuroda 0:13a5d365ba16 45
ykuroda 0:13a5d365ba16 46 template<typename MatrixType, int QRPreconditioner, int Case,
ykuroda 0:13a5d365ba16 47 bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret
ykuroda 0:13a5d365ba16 48 > struct qr_preconditioner_impl {};
ykuroda 0:13a5d365ba16 49
ykuroda 0:13a5d365ba16 50 template<typename MatrixType, int QRPreconditioner, int Case>
ykuroda 0:13a5d365ba16 51 class qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false>
ykuroda 0:13a5d365ba16 52 {
ykuroda 0:13a5d365ba16 53 public:
ykuroda 0:13a5d365ba16 54 typedef typename MatrixType::Index Index;
ykuroda 0:13a5d365ba16 55 void allocate(const JacobiSVD<MatrixType, QRPreconditioner>&) {}
ykuroda 0:13a5d365ba16 56 bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&)
ykuroda 0:13a5d365ba16 57 {
ykuroda 0:13a5d365ba16 58 return false;
ykuroda 0:13a5d365ba16 59 }
ykuroda 0:13a5d365ba16 60 };
ykuroda 0:13a5d365ba16 61
ykuroda 0:13a5d365ba16 62 /*** preconditioner using FullPivHouseholderQR ***/
ykuroda 0:13a5d365ba16 63
ykuroda 0:13a5d365ba16 64 template<typename MatrixType>
ykuroda 0:13a5d365ba16 65 class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
ykuroda 0:13a5d365ba16 66 {
ykuroda 0:13a5d365ba16 67 public:
ykuroda 0:13a5d365ba16 68 typedef typename MatrixType::Index Index;
ykuroda 0:13a5d365ba16 69 typedef typename MatrixType::Scalar Scalar;
ykuroda 0:13a5d365ba16 70 enum
ykuroda 0:13a5d365ba16 71 {
ykuroda 0:13a5d365ba16 72 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ykuroda 0:13a5d365ba16 73 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
ykuroda 0:13a5d365ba16 74 };
ykuroda 0:13a5d365ba16 75 typedef Matrix<Scalar, 1, RowsAtCompileTime, RowMajor, 1, MaxRowsAtCompileTime> WorkspaceType;
ykuroda 0:13a5d365ba16 76
ykuroda 0:13a5d365ba16 77 void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
ykuroda 0:13a5d365ba16 78 {
ykuroda 0:13a5d365ba16 79 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
ykuroda 0:13a5d365ba16 80 {
ykuroda 0:13a5d365ba16 81 m_qr.~QRType();
ykuroda 0:13a5d365ba16 82 ::new (&m_qr) QRType(svd.rows(), svd.cols());
ykuroda 0:13a5d365ba16 83 }
ykuroda 0:13a5d365ba16 84 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
ykuroda 0:13a5d365ba16 85 }
ykuroda 0:13a5d365ba16 86
ykuroda 0:13a5d365ba16 87 bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
ykuroda 0:13a5d365ba16 88 {
ykuroda 0:13a5d365ba16 89 if(matrix.rows() > matrix.cols())
ykuroda 0:13a5d365ba16 90 {
ykuroda 0:13a5d365ba16 91 m_qr.compute(matrix);
ykuroda 0:13a5d365ba16 92 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
ykuroda 0:13a5d365ba16 93 if(svd.m_computeFullU) m_qr.matrixQ().evalTo(svd.m_matrixU, m_workspace);
ykuroda 0:13a5d365ba16 94 if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
ykuroda 0:13a5d365ba16 95 return true;
ykuroda 0:13a5d365ba16 96 }
ykuroda 0:13a5d365ba16 97 return false;
ykuroda 0:13a5d365ba16 98 }
ykuroda 0:13a5d365ba16 99 private:
ykuroda 0:13a5d365ba16 100 typedef FullPivHouseholderQR<MatrixType> QRType;
ykuroda 0:13a5d365ba16 101 QRType m_qr;
ykuroda 0:13a5d365ba16 102 WorkspaceType m_workspace;
ykuroda 0:13a5d365ba16 103 };
ykuroda 0:13a5d365ba16 104
ykuroda 0:13a5d365ba16 105 template<typename MatrixType>
ykuroda 0:13a5d365ba16 106 class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
ykuroda 0:13a5d365ba16 107 {
ykuroda 0:13a5d365ba16 108 public:
ykuroda 0:13a5d365ba16 109 typedef typename MatrixType::Index Index;
ykuroda 0:13a5d365ba16 110 typedef typename MatrixType::Scalar Scalar;
ykuroda 0:13a5d365ba16 111 enum
ykuroda 0:13a5d365ba16 112 {
ykuroda 0:13a5d365ba16 113 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ykuroda 0:13a5d365ba16 114 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
ykuroda 0:13a5d365ba16 115 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
ykuroda 0:13a5d365ba16 116 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
ykuroda 0:13a5d365ba16 117 Options = MatrixType::Options
ykuroda 0:13a5d365ba16 118 };
ykuroda 0:13a5d365ba16 119 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
ykuroda 0:13a5d365ba16 120 TransposeTypeWithSameStorageOrder;
ykuroda 0:13a5d365ba16 121
ykuroda 0:13a5d365ba16 122 void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
ykuroda 0:13a5d365ba16 123 {
ykuroda 0:13a5d365ba16 124 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
ykuroda 0:13a5d365ba16 125 {
ykuroda 0:13a5d365ba16 126 m_qr.~QRType();
ykuroda 0:13a5d365ba16 127 ::new (&m_qr) QRType(svd.cols(), svd.rows());
ykuroda 0:13a5d365ba16 128 }
ykuroda 0:13a5d365ba16 129 m_adjoint.resize(svd.cols(), svd.rows());
ykuroda 0:13a5d365ba16 130 if (svd.m_computeFullV) m_workspace.resize(svd.cols());
ykuroda 0:13a5d365ba16 131 }
ykuroda 0:13a5d365ba16 132
ykuroda 0:13a5d365ba16 133 bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
ykuroda 0:13a5d365ba16 134 {
ykuroda 0:13a5d365ba16 135 if(matrix.cols() > matrix.rows())
ykuroda 0:13a5d365ba16 136 {
ykuroda 0:13a5d365ba16 137 m_adjoint = matrix.adjoint();
ykuroda 0:13a5d365ba16 138 m_qr.compute(m_adjoint);
ykuroda 0:13a5d365ba16 139 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
ykuroda 0:13a5d365ba16 140 if(svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace);
ykuroda 0:13a5d365ba16 141 if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
ykuroda 0:13a5d365ba16 142 return true;
ykuroda 0:13a5d365ba16 143 }
ykuroda 0:13a5d365ba16 144 else return false;
ykuroda 0:13a5d365ba16 145 }
ykuroda 0:13a5d365ba16 146 private:
ykuroda 0:13a5d365ba16 147 typedef FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
ykuroda 0:13a5d365ba16 148 QRType m_qr;
ykuroda 0:13a5d365ba16 149 TransposeTypeWithSameStorageOrder m_adjoint;
ykuroda 0:13a5d365ba16 150 typename internal::plain_row_type<MatrixType>::type m_workspace;
ykuroda 0:13a5d365ba16 151 };
ykuroda 0:13a5d365ba16 152
ykuroda 0:13a5d365ba16 153 /*** preconditioner using ColPivHouseholderQR ***/
ykuroda 0:13a5d365ba16 154
ykuroda 0:13a5d365ba16 155 template<typename MatrixType>
ykuroda 0:13a5d365ba16 156 class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
ykuroda 0:13a5d365ba16 157 {
ykuroda 0:13a5d365ba16 158 public:
ykuroda 0:13a5d365ba16 159 typedef typename MatrixType::Index Index;
ykuroda 0:13a5d365ba16 160
ykuroda 0:13a5d365ba16 161 void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
ykuroda 0:13a5d365ba16 162 {
ykuroda 0:13a5d365ba16 163 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
ykuroda 0:13a5d365ba16 164 {
ykuroda 0:13a5d365ba16 165 m_qr.~QRType();
ykuroda 0:13a5d365ba16 166 ::new (&m_qr) QRType(svd.rows(), svd.cols());
ykuroda 0:13a5d365ba16 167 }
ykuroda 0:13a5d365ba16 168 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
ykuroda 0:13a5d365ba16 169 else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
ykuroda 0:13a5d365ba16 170 }
ykuroda 0:13a5d365ba16 171
ykuroda 0:13a5d365ba16 172 bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
ykuroda 0:13a5d365ba16 173 {
ykuroda 0:13a5d365ba16 174 if(matrix.rows() > matrix.cols())
ykuroda 0:13a5d365ba16 175 {
ykuroda 0:13a5d365ba16 176 m_qr.compute(matrix);
ykuroda 0:13a5d365ba16 177 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
ykuroda 0:13a5d365ba16 178 if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
ykuroda 0:13a5d365ba16 179 else if(svd.m_computeThinU)
ykuroda 0:13a5d365ba16 180 {
ykuroda 0:13a5d365ba16 181 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
ykuroda 0:13a5d365ba16 182 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
ykuroda 0:13a5d365ba16 183 }
ykuroda 0:13a5d365ba16 184 if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
ykuroda 0:13a5d365ba16 185 return true;
ykuroda 0:13a5d365ba16 186 }
ykuroda 0:13a5d365ba16 187 return false;
ykuroda 0:13a5d365ba16 188 }
ykuroda 0:13a5d365ba16 189
ykuroda 0:13a5d365ba16 190 private:
ykuroda 0:13a5d365ba16 191 typedef ColPivHouseholderQR<MatrixType> QRType;
ykuroda 0:13a5d365ba16 192 QRType m_qr;
ykuroda 0:13a5d365ba16 193 typename internal::plain_col_type<MatrixType>::type m_workspace;
ykuroda 0:13a5d365ba16 194 };
ykuroda 0:13a5d365ba16 195
ykuroda 0:13a5d365ba16 196 template<typename MatrixType>
ykuroda 0:13a5d365ba16 197 class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
ykuroda 0:13a5d365ba16 198 {
ykuroda 0:13a5d365ba16 199 public:
ykuroda 0:13a5d365ba16 200 typedef typename MatrixType::Index Index;
ykuroda 0:13a5d365ba16 201 typedef typename MatrixType::Scalar Scalar;
ykuroda 0:13a5d365ba16 202 enum
ykuroda 0:13a5d365ba16 203 {
ykuroda 0:13a5d365ba16 204 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ykuroda 0:13a5d365ba16 205 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
ykuroda 0:13a5d365ba16 206 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
ykuroda 0:13a5d365ba16 207 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
ykuroda 0:13a5d365ba16 208 Options = MatrixType::Options
ykuroda 0:13a5d365ba16 209 };
ykuroda 0:13a5d365ba16 210
ykuroda 0:13a5d365ba16 211 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
ykuroda 0:13a5d365ba16 212 TransposeTypeWithSameStorageOrder;
ykuroda 0:13a5d365ba16 213
ykuroda 0:13a5d365ba16 214 void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
ykuroda 0:13a5d365ba16 215 {
ykuroda 0:13a5d365ba16 216 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
ykuroda 0:13a5d365ba16 217 {
ykuroda 0:13a5d365ba16 218 m_qr.~QRType();
ykuroda 0:13a5d365ba16 219 ::new (&m_qr) QRType(svd.cols(), svd.rows());
ykuroda 0:13a5d365ba16 220 }
ykuroda 0:13a5d365ba16 221 if (svd.m_computeFullV) m_workspace.resize(svd.cols());
ykuroda 0:13a5d365ba16 222 else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
ykuroda 0:13a5d365ba16 223 m_adjoint.resize(svd.cols(), svd.rows());
ykuroda 0:13a5d365ba16 224 }
ykuroda 0:13a5d365ba16 225
ykuroda 0:13a5d365ba16 226 bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
ykuroda 0:13a5d365ba16 227 {
ykuroda 0:13a5d365ba16 228 if(matrix.cols() > matrix.rows())
ykuroda 0:13a5d365ba16 229 {
ykuroda 0:13a5d365ba16 230 m_adjoint = matrix.adjoint();
ykuroda 0:13a5d365ba16 231 m_qr.compute(m_adjoint);
ykuroda 0:13a5d365ba16 232
ykuroda 0:13a5d365ba16 233 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
ykuroda 0:13a5d365ba16 234 if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
ykuroda 0:13a5d365ba16 235 else if(svd.m_computeThinV)
ykuroda 0:13a5d365ba16 236 {
ykuroda 0:13a5d365ba16 237 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
ykuroda 0:13a5d365ba16 238 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
ykuroda 0:13a5d365ba16 239 }
ykuroda 0:13a5d365ba16 240 if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
ykuroda 0:13a5d365ba16 241 return true;
ykuroda 0:13a5d365ba16 242 }
ykuroda 0:13a5d365ba16 243 else return false;
ykuroda 0:13a5d365ba16 244 }
ykuroda 0:13a5d365ba16 245
ykuroda 0:13a5d365ba16 246 private:
ykuroda 0:13a5d365ba16 247 typedef ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
ykuroda 0:13a5d365ba16 248 QRType m_qr;
ykuroda 0:13a5d365ba16 249 TransposeTypeWithSameStorageOrder m_adjoint;
ykuroda 0:13a5d365ba16 250 typename internal::plain_row_type<MatrixType>::type m_workspace;
ykuroda 0:13a5d365ba16 251 };
ykuroda 0:13a5d365ba16 252
ykuroda 0:13a5d365ba16 253 /*** preconditioner using HouseholderQR ***/
ykuroda 0:13a5d365ba16 254
ykuroda 0:13a5d365ba16 255 template<typename MatrixType>
ykuroda 0:13a5d365ba16 256 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
ykuroda 0:13a5d365ba16 257 {
ykuroda 0:13a5d365ba16 258 public:
ykuroda 0:13a5d365ba16 259 typedef typename MatrixType::Index Index;
ykuroda 0:13a5d365ba16 260
ykuroda 0:13a5d365ba16 261 void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
ykuroda 0:13a5d365ba16 262 {
ykuroda 0:13a5d365ba16 263 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
ykuroda 0:13a5d365ba16 264 {
ykuroda 0:13a5d365ba16 265 m_qr.~QRType();
ykuroda 0:13a5d365ba16 266 ::new (&m_qr) QRType(svd.rows(), svd.cols());
ykuroda 0:13a5d365ba16 267 }
ykuroda 0:13a5d365ba16 268 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
ykuroda 0:13a5d365ba16 269 else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
ykuroda 0:13a5d365ba16 270 }
ykuroda 0:13a5d365ba16 271
ykuroda 0:13a5d365ba16 272 bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
ykuroda 0:13a5d365ba16 273 {
ykuroda 0:13a5d365ba16 274 if(matrix.rows() > matrix.cols())
ykuroda 0:13a5d365ba16 275 {
ykuroda 0:13a5d365ba16 276 m_qr.compute(matrix);
ykuroda 0:13a5d365ba16 277 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
ykuroda 0:13a5d365ba16 278 if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
ykuroda 0:13a5d365ba16 279 else if(svd.m_computeThinU)
ykuroda 0:13a5d365ba16 280 {
ykuroda 0:13a5d365ba16 281 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
ykuroda 0:13a5d365ba16 282 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
ykuroda 0:13a5d365ba16 283 }
ykuroda 0:13a5d365ba16 284 if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
ykuroda 0:13a5d365ba16 285 return true;
ykuroda 0:13a5d365ba16 286 }
ykuroda 0:13a5d365ba16 287 return false;
ykuroda 0:13a5d365ba16 288 }
ykuroda 0:13a5d365ba16 289 private:
ykuroda 0:13a5d365ba16 290 typedef HouseholderQR<MatrixType> QRType;
ykuroda 0:13a5d365ba16 291 QRType m_qr;
ykuroda 0:13a5d365ba16 292 typename internal::plain_col_type<MatrixType>::type m_workspace;
ykuroda 0:13a5d365ba16 293 };
ykuroda 0:13a5d365ba16 294
ykuroda 0:13a5d365ba16 295 template<typename MatrixType>
ykuroda 0:13a5d365ba16 296 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
ykuroda 0:13a5d365ba16 297 {
ykuroda 0:13a5d365ba16 298 public:
ykuroda 0:13a5d365ba16 299 typedef typename MatrixType::Index Index;
ykuroda 0:13a5d365ba16 300 typedef typename MatrixType::Scalar Scalar;
ykuroda 0:13a5d365ba16 301 enum
ykuroda 0:13a5d365ba16 302 {
ykuroda 0:13a5d365ba16 303 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ykuroda 0:13a5d365ba16 304 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
ykuroda 0:13a5d365ba16 305 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
ykuroda 0:13a5d365ba16 306 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
ykuroda 0:13a5d365ba16 307 Options = MatrixType::Options
ykuroda 0:13a5d365ba16 308 };
ykuroda 0:13a5d365ba16 309
ykuroda 0:13a5d365ba16 310 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
ykuroda 0:13a5d365ba16 311 TransposeTypeWithSameStorageOrder;
ykuroda 0:13a5d365ba16 312
ykuroda 0:13a5d365ba16 313 void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
ykuroda 0:13a5d365ba16 314 {
ykuroda 0:13a5d365ba16 315 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
ykuroda 0:13a5d365ba16 316 {
ykuroda 0:13a5d365ba16 317 m_qr.~QRType();
ykuroda 0:13a5d365ba16 318 ::new (&m_qr) QRType(svd.cols(), svd.rows());
ykuroda 0:13a5d365ba16 319 }
ykuroda 0:13a5d365ba16 320 if (svd.m_computeFullV) m_workspace.resize(svd.cols());
ykuroda 0:13a5d365ba16 321 else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
ykuroda 0:13a5d365ba16 322 m_adjoint.resize(svd.cols(), svd.rows());
ykuroda 0:13a5d365ba16 323 }
ykuroda 0:13a5d365ba16 324
ykuroda 0:13a5d365ba16 325 bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
ykuroda 0:13a5d365ba16 326 {
ykuroda 0:13a5d365ba16 327 if(matrix.cols() > matrix.rows())
ykuroda 0:13a5d365ba16 328 {
ykuroda 0:13a5d365ba16 329 m_adjoint = matrix.adjoint();
ykuroda 0:13a5d365ba16 330 m_qr.compute(m_adjoint);
ykuroda 0:13a5d365ba16 331
ykuroda 0:13a5d365ba16 332 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
ykuroda 0:13a5d365ba16 333 if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
ykuroda 0:13a5d365ba16 334 else if(svd.m_computeThinV)
ykuroda 0:13a5d365ba16 335 {
ykuroda 0:13a5d365ba16 336 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
ykuroda 0:13a5d365ba16 337 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
ykuroda 0:13a5d365ba16 338 }
ykuroda 0:13a5d365ba16 339 if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
ykuroda 0:13a5d365ba16 340 return true;
ykuroda 0:13a5d365ba16 341 }
ykuroda 0:13a5d365ba16 342 else return false;
ykuroda 0:13a5d365ba16 343 }
ykuroda 0:13a5d365ba16 344
ykuroda 0:13a5d365ba16 345 private:
ykuroda 0:13a5d365ba16 346 typedef HouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
ykuroda 0:13a5d365ba16 347 QRType m_qr;
ykuroda 0:13a5d365ba16 348 TransposeTypeWithSameStorageOrder m_adjoint;
ykuroda 0:13a5d365ba16 349 typename internal::plain_row_type<MatrixType>::type m_workspace;
ykuroda 0:13a5d365ba16 350 };
ykuroda 0:13a5d365ba16 351
ykuroda 0:13a5d365ba16 352 /*** 2x2 SVD implementation
ykuroda 0:13a5d365ba16 353 ***
ykuroda 0:13a5d365ba16 354 *** JacobiSVD consists in performing a series of 2x2 SVD subproblems
ykuroda 0:13a5d365ba16 355 ***/
ykuroda 0:13a5d365ba16 356
ykuroda 0:13a5d365ba16 357 template<typename MatrixType, int QRPreconditioner>
ykuroda 0:13a5d365ba16 358 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false>
ykuroda 0:13a5d365ba16 359 {
ykuroda 0:13a5d365ba16 360 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
ykuroda 0:13a5d365ba16 361 typedef typename SVD::Index Index;
ykuroda 0:13a5d365ba16 362 static void run(typename SVD::WorkMatrixType&, SVD&, Index, Index) {}
ykuroda 0:13a5d365ba16 363 };
ykuroda 0:13a5d365ba16 364
ykuroda 0:13a5d365ba16 365 template<typename MatrixType, int QRPreconditioner>
ykuroda 0:13a5d365ba16 366 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
ykuroda 0:13a5d365ba16 367 {
ykuroda 0:13a5d365ba16 368 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
ykuroda 0:13a5d365ba16 369 typedef typename MatrixType::Scalar Scalar;
ykuroda 0:13a5d365ba16 370 typedef typename MatrixType::RealScalar RealScalar;
ykuroda 0:13a5d365ba16 371 typedef typename SVD::Index Index;
ykuroda 0:13a5d365ba16 372 static void run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q)
ykuroda 0:13a5d365ba16 373 {
ykuroda 0:13a5d365ba16 374 using std::sqrt;
ykuroda 0:13a5d365ba16 375 Scalar z;
ykuroda 0:13a5d365ba16 376 JacobiRotation<Scalar> rot;
ykuroda 0:13a5d365ba16 377 RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p)));
ykuroda 0:13a5d365ba16 378
ykuroda 0:13a5d365ba16 379 if(n==0)
ykuroda 0:13a5d365ba16 380 {
ykuroda 0:13a5d365ba16 381 z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
ykuroda 0:13a5d365ba16 382 work_matrix.row(p) *= z;
ykuroda 0:13a5d365ba16 383 if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z);
ykuroda 0:13a5d365ba16 384 if(work_matrix.coeff(q,q)!=Scalar(0))
ykuroda 0:13a5d365ba16 385 {
ykuroda 0:13a5d365ba16 386 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
ykuroda 0:13a5d365ba16 387 work_matrix.row(q) *= z;
ykuroda 0:13a5d365ba16 388 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
ykuroda 0:13a5d365ba16 389 }
ykuroda 0:13a5d365ba16 390 // otherwise the second row is already zero, so we have nothing to do.
ykuroda 0:13a5d365ba16 391 }
ykuroda 0:13a5d365ba16 392 else
ykuroda 0:13a5d365ba16 393 {
ykuroda 0:13a5d365ba16 394 rot.c() = conj(work_matrix.coeff(p,p)) / n;
ykuroda 0:13a5d365ba16 395 rot.s() = work_matrix.coeff(q,p) / n;
ykuroda 0:13a5d365ba16 396 work_matrix.applyOnTheLeft(p,q,rot);
ykuroda 0:13a5d365ba16 397 if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint());
ykuroda 0:13a5d365ba16 398 if(work_matrix.coeff(p,q) != Scalar(0))
ykuroda 0:13a5d365ba16 399 {
ykuroda 0:13a5d365ba16 400 Scalar z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
ykuroda 0:13a5d365ba16 401 work_matrix.col(q) *= z;
ykuroda 0:13a5d365ba16 402 if(svd.computeV()) svd.m_matrixV.col(q) *= z;
ykuroda 0:13a5d365ba16 403 }
ykuroda 0:13a5d365ba16 404 if(work_matrix.coeff(q,q) != Scalar(0))
ykuroda 0:13a5d365ba16 405 {
ykuroda 0:13a5d365ba16 406 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
ykuroda 0:13a5d365ba16 407 work_matrix.row(q) *= z;
ykuroda 0:13a5d365ba16 408 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
ykuroda 0:13a5d365ba16 409 }
ykuroda 0:13a5d365ba16 410 }
ykuroda 0:13a5d365ba16 411 }
ykuroda 0:13a5d365ba16 412 };
ykuroda 0:13a5d365ba16 413
ykuroda 0:13a5d365ba16 414 template<typename MatrixType, typename RealScalar, typename Index>
ykuroda 0:13a5d365ba16 415 void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
ykuroda 0:13a5d365ba16 416 JacobiRotation<RealScalar> *j_left,
ykuroda 0:13a5d365ba16 417 JacobiRotation<RealScalar> *j_right)
ykuroda 0:13a5d365ba16 418 {
ykuroda 0:13a5d365ba16 419 using std::sqrt;
ykuroda 0:13a5d365ba16 420 using std::abs;
ykuroda 0:13a5d365ba16 421 Matrix<RealScalar,2,2> m;
ykuroda 0:13a5d365ba16 422 m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)),
ykuroda 0:13a5d365ba16 423 numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q));
ykuroda 0:13a5d365ba16 424 JacobiRotation<RealScalar> rot1;
ykuroda 0:13a5d365ba16 425 RealScalar t = m.coeff(0,0) + m.coeff(1,1);
ykuroda 0:13a5d365ba16 426 RealScalar d = m.coeff(1,0) - m.coeff(0,1);
ykuroda 0:13a5d365ba16 427 if(t == RealScalar(0))
ykuroda 0:13a5d365ba16 428 {
ykuroda 0:13a5d365ba16 429 rot1.c() = RealScalar(0);
ykuroda 0:13a5d365ba16 430 rot1.s() = d > RealScalar(0) ? RealScalar(1) : RealScalar(-1);
ykuroda 0:13a5d365ba16 431 }
ykuroda 0:13a5d365ba16 432 else
ykuroda 0:13a5d365ba16 433 {
ykuroda 0:13a5d365ba16 434 RealScalar t2d2 = numext::hypot(t,d);
ykuroda 0:13a5d365ba16 435 rot1.c() = abs(t)/t2d2;
ykuroda 0:13a5d365ba16 436 rot1.s() = d/t2d2;
ykuroda 0:13a5d365ba16 437 if(t<RealScalar(0))
ykuroda 0:13a5d365ba16 438 rot1.s() = -rot1.s();
ykuroda 0:13a5d365ba16 439 }
ykuroda 0:13a5d365ba16 440 m.applyOnTheLeft(0,1,rot1);
ykuroda 0:13a5d365ba16 441 j_right->makeJacobi(m,0,1);
ykuroda 0:13a5d365ba16 442 *j_left = rot1 * j_right->transpose();
ykuroda 0:13a5d365ba16 443 }
ykuroda 0:13a5d365ba16 444
ykuroda 0:13a5d365ba16 445 } // end namespace internal
ykuroda 0:13a5d365ba16 446
ykuroda 0:13a5d365ba16 447 /** \ingroup SVD_Module
ykuroda 0:13a5d365ba16 448 *
ykuroda 0:13a5d365ba16 449 *
ykuroda 0:13a5d365ba16 450 * \class JacobiSVD
ykuroda 0:13a5d365ba16 451 *
ykuroda 0:13a5d365ba16 452 * \brief Two-sided Jacobi SVD decomposition of a rectangular matrix
ykuroda 0:13a5d365ba16 453 *
ykuroda 0:13a5d365ba16 454 * \param MatrixType the type of the matrix of which we are computing the SVD decomposition
ykuroda 0:13a5d365ba16 455 * \param QRPreconditioner this optional parameter allows to specify the type of QR decomposition that will be used internally
ykuroda 0:13a5d365ba16 456 * for the R-SVD step for non-square matrices. See discussion of possible values below.
ykuroda 0:13a5d365ba16 457 *
ykuroda 0:13a5d365ba16 458 * SVD decomposition consists in decomposing any n-by-p matrix \a A as a product
ykuroda 0:13a5d365ba16 459 * \f[ A = U S V^* \f]
ykuroda 0:13a5d365ba16 460 * where \a U is a n-by-n unitary, \a V is a p-by-p unitary, and \a S is a n-by-p real positive matrix which is zero outside of its main diagonal;
ykuroda 0:13a5d365ba16 461 * the diagonal entries of S are known as the \em singular \em values of \a A and the columns of \a U and \a V are known as the left
ykuroda 0:13a5d365ba16 462 * and right \em singular \em vectors of \a A respectively.
ykuroda 0:13a5d365ba16 463 *
ykuroda 0:13a5d365ba16 464 * Singular values are always sorted in decreasing order.
ykuroda 0:13a5d365ba16 465 *
ykuroda 0:13a5d365ba16 466 * This JacobiSVD decomposition computes only the singular values by default. If you want \a U or \a V, you need to ask for them explicitly.
ykuroda 0:13a5d365ba16 467 *
ykuroda 0:13a5d365ba16 468 * You can ask for only \em thin \a U or \a V to be computed, meaning the following. In case of a rectangular n-by-p matrix, letting \a m be the
ykuroda 0:13a5d365ba16 469 * smaller value among \a n and \a p, there are only \a m singular vectors; the remaining columns of \a U and \a V do not correspond to actual
ykuroda 0:13a5d365ba16 470 * singular vectors. Asking for \em thin \a U or \a V means asking for only their \a m first columns to be formed. So \a U is then a n-by-m matrix,
ykuroda 0:13a5d365ba16 471 * and \a V is then a p-by-m matrix. Notice that thin \a U and \a V are all you need for (least squares) solving.
ykuroda 0:13a5d365ba16 472 *
ykuroda 0:13a5d365ba16 473 * Here's an example demonstrating basic usage:
ykuroda 0:13a5d365ba16 474 * \include JacobiSVD_basic.cpp
ykuroda 0:13a5d365ba16 475 * Output: \verbinclude JacobiSVD_basic.out
ykuroda 0:13a5d365ba16 476 *
ykuroda 0:13a5d365ba16 477 * This JacobiSVD class is a two-sided Jacobi R-SVD decomposition, ensuring optimal reliability and accuracy. The downside is that it's slower than
ykuroda 0:13a5d365ba16 478 * bidiagonalizing SVD algorithms for large square matrices; however its complexity is still \f$ O(n^2p) \f$ where \a n is the smaller dimension and
ykuroda 0:13a5d365ba16 479 * \a p is the greater dimension, meaning that it is still of the same order of complexity as the faster bidiagonalizing R-SVD algorithms.
ykuroda 0:13a5d365ba16 480 * In particular, like any R-SVD, it takes advantage of non-squareness in that its complexity is only linear in the greater dimension.
ykuroda 0:13a5d365ba16 481 *
ykuroda 0:13a5d365ba16 482 * If the input matrix has inf or nan coefficients, the result of the computation is undefined, but the computation is guaranteed to
ykuroda 0:13a5d365ba16 483 * terminate in finite (and reasonable) time.
ykuroda 0:13a5d365ba16 484 *
ykuroda 0:13a5d365ba16 485 * The possible values for QRPreconditioner are:
ykuroda 0:13a5d365ba16 486 * \li ColPivHouseholderQRPreconditioner is the default. In practice it's very safe. It uses column-pivoting QR.
ykuroda 0:13a5d365ba16 487 * \li FullPivHouseholderQRPreconditioner, is the safest and slowest. It uses full-pivoting QR.
ykuroda 0:13a5d365ba16 488 * Contrary to other QRs, it doesn't allow computing thin unitaries.
ykuroda 0:13a5d365ba16 489 * \li HouseholderQRPreconditioner is the fastest, and less safe and accurate than the pivoting variants. It uses non-pivoting QR.
ykuroda 0:13a5d365ba16 490 * This is very similar in safety and accuracy to the bidiagonalization process used by bidiagonalizing SVD algorithms (since bidiagonalization
ykuroda 0:13a5d365ba16 491 * is inherently non-pivoting). However the resulting SVD is still more reliable than bidiagonalizing SVDs because the Jacobi-based iterarive
ykuroda 0:13a5d365ba16 492 * process is more reliable than the optimized bidiagonal SVD iterations.
ykuroda 0:13a5d365ba16 493 * \li NoQRPreconditioner allows not to use a QR preconditioner at all. This is useful if you know that you will only be computing
ykuroda 0:13a5d365ba16 494 * JacobiSVD decompositions of square matrices. Non-square matrices require a QR preconditioner. Using this option will result in
ykuroda 0:13a5d365ba16 495 * faster compilation and smaller executable code. It won't significantly speed up computation, since JacobiSVD is always checking
ykuroda 0:13a5d365ba16 496 * if QR preconditioning is needed before applying it anyway.
ykuroda 0:13a5d365ba16 497 *
ykuroda 0:13a5d365ba16 498 * \sa MatrixBase::jacobiSvd()
ykuroda 0:13a5d365ba16 499 */
ykuroda 0:13a5d365ba16 500 template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
ykuroda 0:13a5d365ba16 501 {
ykuroda 0:13a5d365ba16 502 public:
ykuroda 0:13a5d365ba16 503
ykuroda 0:13a5d365ba16 504 typedef _MatrixType MatrixType;
ykuroda 0:13a5d365ba16 505 typedef typename MatrixType::Scalar Scalar;
ykuroda 0:13a5d365ba16 506 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
ykuroda 0:13a5d365ba16 507 typedef typename MatrixType::Index Index;
ykuroda 0:13a5d365ba16 508 enum {
ykuroda 0:13a5d365ba16 509 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ykuroda 0:13a5d365ba16 510 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
ykuroda 0:13a5d365ba16 511 DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime),
ykuroda 0:13a5d365ba16 512 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
ykuroda 0:13a5d365ba16 513 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
ykuroda 0:13a5d365ba16 514 MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime),
ykuroda 0:13a5d365ba16 515 MatrixOptions = MatrixType::Options
ykuroda 0:13a5d365ba16 516 };
ykuroda 0:13a5d365ba16 517
ykuroda 0:13a5d365ba16 518 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime,
ykuroda 0:13a5d365ba16 519 MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime>
ykuroda 0:13a5d365ba16 520 MatrixUType;
ykuroda 0:13a5d365ba16 521 typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime,
ykuroda 0:13a5d365ba16 522 MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime>
ykuroda 0:13a5d365ba16 523 MatrixVType;
ykuroda 0:13a5d365ba16 524 typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType;
ykuroda 0:13a5d365ba16 525 typedef typename internal::plain_row_type<MatrixType>::type RowType;
ykuroda 0:13a5d365ba16 526 typedef typename internal::plain_col_type<MatrixType>::type ColType;
ykuroda 0:13a5d365ba16 527 typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime,
ykuroda 0:13a5d365ba16 528 MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime>
ykuroda 0:13a5d365ba16 529 WorkMatrixType;
ykuroda 0:13a5d365ba16 530
ykuroda 0:13a5d365ba16 531 /** \brief Default Constructor.
ykuroda 0:13a5d365ba16 532 *
ykuroda 0:13a5d365ba16 533 * The default constructor is useful in cases in which the user intends to
ykuroda 0:13a5d365ba16 534 * perform decompositions via JacobiSVD::compute(const MatrixType&).
ykuroda 0:13a5d365ba16 535 */
ykuroda 0:13a5d365ba16 536 JacobiSVD()
ykuroda 0:13a5d365ba16 537 : m_isInitialized(false),
ykuroda 0:13a5d365ba16 538 m_isAllocated(false),
ykuroda 0:13a5d365ba16 539 m_usePrescribedThreshold(false),
ykuroda 0:13a5d365ba16 540 m_computationOptions(0),
ykuroda 0:13a5d365ba16 541 m_rows(-1), m_cols(-1), m_diagSize(0)
ykuroda 0:13a5d365ba16 542 {}
ykuroda 0:13a5d365ba16 543
ykuroda 0:13a5d365ba16 544
ykuroda 0:13a5d365ba16 545 /** \brief Default Constructor with memory preallocation
ykuroda 0:13a5d365ba16 546 *
ykuroda 0:13a5d365ba16 547 * Like the default constructor but with preallocation of the internal data
ykuroda 0:13a5d365ba16 548 * according to the specified problem size.
ykuroda 0:13a5d365ba16 549 * \sa JacobiSVD()
ykuroda 0:13a5d365ba16 550 */
ykuroda 0:13a5d365ba16 551 JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0)
ykuroda 0:13a5d365ba16 552 : m_isInitialized(false),
ykuroda 0:13a5d365ba16 553 m_isAllocated(false),
ykuroda 0:13a5d365ba16 554 m_usePrescribedThreshold(false),
ykuroda 0:13a5d365ba16 555 m_computationOptions(0),
ykuroda 0:13a5d365ba16 556 m_rows(-1), m_cols(-1)
ykuroda 0:13a5d365ba16 557 {
ykuroda 0:13a5d365ba16 558 allocate(rows, cols, computationOptions);
ykuroda 0:13a5d365ba16 559 }
ykuroda 0:13a5d365ba16 560
ykuroda 0:13a5d365ba16 561 /** \brief Constructor performing the decomposition of given matrix.
ykuroda 0:13a5d365ba16 562 *
ykuroda 0:13a5d365ba16 563 * \param matrix the matrix to decompose
ykuroda 0:13a5d365ba16 564 * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed.
ykuroda 0:13a5d365ba16 565 * By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU,
ykuroda 0:13a5d365ba16 566 * #ComputeFullV, #ComputeThinV.
ykuroda 0:13a5d365ba16 567 *
ykuroda 0:13a5d365ba16 568 * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not
ykuroda 0:13a5d365ba16 569 * available with the (non-default) FullPivHouseholderQR preconditioner.
ykuroda 0:13a5d365ba16 570 */
ykuroda 0:13a5d365ba16 571 JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0)
ykuroda 0:13a5d365ba16 572 : m_isInitialized(false),
ykuroda 0:13a5d365ba16 573 m_isAllocated(false),
ykuroda 0:13a5d365ba16 574 m_usePrescribedThreshold(false),
ykuroda 0:13a5d365ba16 575 m_computationOptions(0),
ykuroda 0:13a5d365ba16 576 m_rows(-1), m_cols(-1)
ykuroda 0:13a5d365ba16 577 {
ykuroda 0:13a5d365ba16 578 compute(matrix, computationOptions);
ykuroda 0:13a5d365ba16 579 }
ykuroda 0:13a5d365ba16 580
ykuroda 0:13a5d365ba16 581 /** \brief Method performing the decomposition of given matrix using custom options.
ykuroda 0:13a5d365ba16 582 *
ykuroda 0:13a5d365ba16 583 * \param matrix the matrix to decompose
ykuroda 0:13a5d365ba16 584 * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed.
ykuroda 0:13a5d365ba16 585 * By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU,
ykuroda 0:13a5d365ba16 586 * #ComputeFullV, #ComputeThinV.
ykuroda 0:13a5d365ba16 587 *
ykuroda 0:13a5d365ba16 588 * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not
ykuroda 0:13a5d365ba16 589 * available with the (non-default) FullPivHouseholderQR preconditioner.
ykuroda 0:13a5d365ba16 590 */
ykuroda 0:13a5d365ba16 591 JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions);
ykuroda 0:13a5d365ba16 592
ykuroda 0:13a5d365ba16 593 /** \brief Method performing the decomposition of given matrix using current options.
ykuroda 0:13a5d365ba16 594 *
ykuroda 0:13a5d365ba16 595 * \param matrix the matrix to decompose
ykuroda 0:13a5d365ba16 596 *
ykuroda 0:13a5d365ba16 597 * This method uses the current \a computationOptions, as already passed to the constructor or to compute(const MatrixType&, unsigned int).
ykuroda 0:13a5d365ba16 598 */
ykuroda 0:13a5d365ba16 599 JacobiSVD& compute(const MatrixType& matrix)
ykuroda 0:13a5d365ba16 600 {
ykuroda 0:13a5d365ba16 601 return compute(matrix, m_computationOptions);
ykuroda 0:13a5d365ba16 602 }
ykuroda 0:13a5d365ba16 603
ykuroda 0:13a5d365ba16 604 /** \returns the \a U matrix.
ykuroda 0:13a5d365ba16 605 *
ykuroda 0:13a5d365ba16 606 * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p,
ykuroda 0:13a5d365ba16 607 * the U matrix is n-by-n if you asked for #ComputeFullU, and is n-by-m if you asked for #ComputeThinU.
ykuroda 0:13a5d365ba16 608 *
ykuroda 0:13a5d365ba16 609 * The \a m first columns of \a U are the left singular vectors of the matrix being decomposed.
ykuroda 0:13a5d365ba16 610 *
ykuroda 0:13a5d365ba16 611 * This method asserts that you asked for \a U to be computed.
ykuroda 0:13a5d365ba16 612 */
ykuroda 0:13a5d365ba16 613 const MatrixUType& matrixU() const
ykuroda 0:13a5d365ba16 614 {
ykuroda 0:13a5d365ba16 615 eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
ykuroda 0:13a5d365ba16 616 eigen_assert(computeU() && "This JacobiSVD decomposition didn't compute U. Did you ask for it?");
ykuroda 0:13a5d365ba16 617 return m_matrixU;
ykuroda 0:13a5d365ba16 618 }
ykuroda 0:13a5d365ba16 619
ykuroda 0:13a5d365ba16 620 /** \returns the \a V matrix.
ykuroda 0:13a5d365ba16 621 *
ykuroda 0:13a5d365ba16 622 * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p,
ykuroda 0:13a5d365ba16 623 * the V matrix is p-by-p if you asked for #ComputeFullV, and is p-by-m if you asked for ComputeThinV.
ykuroda 0:13a5d365ba16 624 *
ykuroda 0:13a5d365ba16 625 * The \a m first columns of \a V are the right singular vectors of the matrix being decomposed.
ykuroda 0:13a5d365ba16 626 *
ykuroda 0:13a5d365ba16 627 * This method asserts that you asked for \a V to be computed.
ykuroda 0:13a5d365ba16 628 */
ykuroda 0:13a5d365ba16 629 const MatrixVType& matrixV() const
ykuroda 0:13a5d365ba16 630 {
ykuroda 0:13a5d365ba16 631 eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
ykuroda 0:13a5d365ba16 632 eigen_assert(computeV() && "This JacobiSVD decomposition didn't compute V. Did you ask for it?");
ykuroda 0:13a5d365ba16 633 return m_matrixV;
ykuroda 0:13a5d365ba16 634 }
ykuroda 0:13a5d365ba16 635
ykuroda 0:13a5d365ba16 636 /** \returns the vector of singular values.
ykuroda 0:13a5d365ba16 637 *
ykuroda 0:13a5d365ba16 638 * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p, the
ykuroda 0:13a5d365ba16 639 * returned vector has size \a m. Singular values are always sorted in decreasing order.
ykuroda 0:13a5d365ba16 640 */
ykuroda 0:13a5d365ba16 641 const SingularValuesType& singularValues() const
ykuroda 0:13a5d365ba16 642 {
ykuroda 0:13a5d365ba16 643 eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
ykuroda 0:13a5d365ba16 644 return m_singularValues;
ykuroda 0:13a5d365ba16 645 }
ykuroda 0:13a5d365ba16 646
ykuroda 0:13a5d365ba16 647 /** \returns true if \a U (full or thin) is asked for in this SVD decomposition */
ykuroda 0:13a5d365ba16 648 inline bool computeU() const { return m_computeFullU || m_computeThinU; }
ykuroda 0:13a5d365ba16 649 /** \returns true if \a V (full or thin) is asked for in this SVD decomposition */
ykuroda 0:13a5d365ba16 650 inline bool computeV() const { return m_computeFullV || m_computeThinV; }
ykuroda 0:13a5d365ba16 651
ykuroda 0:13a5d365ba16 652 /** \returns a (least squares) solution of \f$ A x = b \f$ using the current SVD decomposition of A.
ykuroda 0:13a5d365ba16 653 *
ykuroda 0:13a5d365ba16 654 * \param b the right-hand-side of the equation to solve.
ykuroda 0:13a5d365ba16 655 *
ykuroda 0:13a5d365ba16 656 * \note Solving requires both U and V to be computed. Thin U and V are enough, there is no need for full U or V.
ykuroda 0:13a5d365ba16 657 *
ykuroda 0:13a5d365ba16 658 * \note SVD solving is implicitly least-squares. Thus, this method serves both purposes of exact solving and least-squares solving.
ykuroda 0:13a5d365ba16 659 * In other words, the returned solution is guaranteed to minimize the Euclidean norm \f$ \Vert A x - b \Vert \f$.
ykuroda 0:13a5d365ba16 660 */
ykuroda 0:13a5d365ba16 661 template<typename Rhs>
ykuroda 0:13a5d365ba16 662 inline const internal::solve_retval<JacobiSVD, Rhs>
ykuroda 0:13a5d365ba16 663 solve(const MatrixBase<Rhs>& b) const
ykuroda 0:13a5d365ba16 664 {
ykuroda 0:13a5d365ba16 665 eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
ykuroda 0:13a5d365ba16 666 eigen_assert(computeU() && computeV() && "JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice).");
ykuroda 0:13a5d365ba16 667 return internal::solve_retval<JacobiSVD, Rhs>(*this, b.derived());
ykuroda 0:13a5d365ba16 668 }
ykuroda 0:13a5d365ba16 669
ykuroda 0:13a5d365ba16 670 /** \returns the number of singular values that are not exactly 0 */
ykuroda 0:13a5d365ba16 671 Index nonzeroSingularValues() const
ykuroda 0:13a5d365ba16 672 {
ykuroda 0:13a5d365ba16 673 eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
ykuroda 0:13a5d365ba16 674 return m_nonzeroSingularValues;
ykuroda 0:13a5d365ba16 675 }
ykuroda 0:13a5d365ba16 676
ykuroda 0:13a5d365ba16 677 /** \returns the rank of the matrix of which \c *this is the SVD.
ykuroda 0:13a5d365ba16 678 *
ykuroda 0:13a5d365ba16 679 * \note This method has to determine which singular values should be considered nonzero.
ykuroda 0:13a5d365ba16 680 * For that, it uses the threshold value that you can control by calling
ykuroda 0:13a5d365ba16 681 * setThreshold(const RealScalar&).
ykuroda 0:13a5d365ba16 682 */
ykuroda 0:13a5d365ba16 683 inline Index rank() const
ykuroda 0:13a5d365ba16 684 {
ykuroda 0:13a5d365ba16 685 using std::abs;
ykuroda 0:13a5d365ba16 686 eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
ykuroda 0:13a5d365ba16 687 if(m_singularValues.size()==0) return 0;
ykuroda 0:13a5d365ba16 688 RealScalar premultiplied_threshold = m_singularValues.coeff(0) * threshold();
ykuroda 0:13a5d365ba16 689 Index i = m_nonzeroSingularValues-1;
ykuroda 0:13a5d365ba16 690 while(i>=0 && m_singularValues.coeff(i) < premultiplied_threshold) --i;
ykuroda 0:13a5d365ba16 691 return i+1;
ykuroda 0:13a5d365ba16 692 }
ykuroda 0:13a5d365ba16 693
ykuroda 0:13a5d365ba16 694 /** Allows to prescribe a threshold to be used by certain methods, such as rank() and solve(),
ykuroda 0:13a5d365ba16 695 * which need to determine when singular values are to be considered nonzero.
ykuroda 0:13a5d365ba16 696 * This is not used for the SVD decomposition itself.
ykuroda 0:13a5d365ba16 697 *
ykuroda 0:13a5d365ba16 698 * When it needs to get the threshold value, Eigen calls threshold().
ykuroda 0:13a5d365ba16 699 * The default is \c NumTraits<Scalar>::epsilon()
ykuroda 0:13a5d365ba16 700 *
ykuroda 0:13a5d365ba16 701 * \param threshold The new value to use as the threshold.
ykuroda 0:13a5d365ba16 702 *
ykuroda 0:13a5d365ba16 703 * A singular value will be considered nonzero if its value is strictly greater than
ykuroda 0:13a5d365ba16 704 * \f$ \vert singular value \vert \leqslant threshold \times \vert max singular value \vert \f$.
ykuroda 0:13a5d365ba16 705 *
ykuroda 0:13a5d365ba16 706 * If you want to come back to the default behavior, call setThreshold(Default_t)
ykuroda 0:13a5d365ba16 707 */
ykuroda 0:13a5d365ba16 708 JacobiSVD& setThreshold(const RealScalar& threshold)
ykuroda 0:13a5d365ba16 709 {
ykuroda 0:13a5d365ba16 710 m_usePrescribedThreshold = true;
ykuroda 0:13a5d365ba16 711 m_prescribedThreshold = threshold;
ykuroda 0:13a5d365ba16 712 return *this;
ykuroda 0:13a5d365ba16 713 }
ykuroda 0:13a5d365ba16 714
ykuroda 0:13a5d365ba16 715 /** Allows to come back to the default behavior, letting Eigen use its default formula for
ykuroda 0:13a5d365ba16 716 * determining the threshold.
ykuroda 0:13a5d365ba16 717 *
ykuroda 0:13a5d365ba16 718 * You should pass the special object Eigen::Default as parameter here.
ykuroda 0:13a5d365ba16 719 * \code svd.setThreshold(Eigen::Default); \endcode
ykuroda 0:13a5d365ba16 720 *
ykuroda 0:13a5d365ba16 721 * See the documentation of setThreshold(const RealScalar&).
ykuroda 0:13a5d365ba16 722 */
ykuroda 0:13a5d365ba16 723 JacobiSVD& setThreshold(Default_t)
ykuroda 0:13a5d365ba16 724 {
ykuroda 0:13a5d365ba16 725 m_usePrescribedThreshold = false;
ykuroda 0:13a5d365ba16 726 return *this;
ykuroda 0:13a5d365ba16 727 }
ykuroda 0:13a5d365ba16 728
ykuroda 0:13a5d365ba16 729 /** Returns the threshold that will be used by certain methods such as rank().
ykuroda 0:13a5d365ba16 730 *
ykuroda 0:13a5d365ba16 731 * See the documentation of setThreshold(const RealScalar&).
ykuroda 0:13a5d365ba16 732 */
ykuroda 0:13a5d365ba16 733 RealScalar threshold() const
ykuroda 0:13a5d365ba16 734 {
ykuroda 0:13a5d365ba16 735 eigen_assert(m_isInitialized || m_usePrescribedThreshold);
ykuroda 0:13a5d365ba16 736 return m_usePrescribedThreshold ? m_prescribedThreshold
ykuroda 0:13a5d365ba16 737 : (std::max<Index>)(1,m_diagSize)*NumTraits<Scalar>::epsilon();
ykuroda 0:13a5d365ba16 738 }
ykuroda 0:13a5d365ba16 739
ykuroda 0:13a5d365ba16 740 inline Index rows() const { return m_rows; }
ykuroda 0:13a5d365ba16 741 inline Index cols() const { return m_cols; }
ykuroda 0:13a5d365ba16 742
ykuroda 0:13a5d365ba16 743 private:
ykuroda 0:13a5d365ba16 744 void allocate(Index rows, Index cols, unsigned int computationOptions);
ykuroda 0:13a5d365ba16 745
ykuroda 0:13a5d365ba16 746 static void check_template_parameters()
ykuroda 0:13a5d365ba16 747 {
ykuroda 0:13a5d365ba16 748 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
ykuroda 0:13a5d365ba16 749 }
ykuroda 0:13a5d365ba16 750
ykuroda 0:13a5d365ba16 751 protected:
ykuroda 0:13a5d365ba16 752 MatrixUType m_matrixU;
ykuroda 0:13a5d365ba16 753 MatrixVType m_matrixV;
ykuroda 0:13a5d365ba16 754 SingularValuesType m_singularValues;
ykuroda 0:13a5d365ba16 755 WorkMatrixType m_workMatrix;
ykuroda 0:13a5d365ba16 756 bool m_isInitialized, m_isAllocated, m_usePrescribedThreshold;
ykuroda 0:13a5d365ba16 757 bool m_computeFullU, m_computeThinU;
ykuroda 0:13a5d365ba16 758 bool m_computeFullV, m_computeThinV;
ykuroda 0:13a5d365ba16 759 unsigned int m_computationOptions;
ykuroda 0:13a5d365ba16 760 Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
ykuroda 0:13a5d365ba16 761 RealScalar m_prescribedThreshold;
ykuroda 0:13a5d365ba16 762
ykuroda 0:13a5d365ba16 763 template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex>
ykuroda 0:13a5d365ba16 764 friend struct internal::svd_precondition_2x2_block_to_be_real;
ykuroda 0:13a5d365ba16 765 template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything>
ykuroda 0:13a5d365ba16 766 friend struct internal::qr_preconditioner_impl;
ykuroda 0:13a5d365ba16 767
ykuroda 0:13a5d365ba16 768 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols;
ykuroda 0:13a5d365ba16 769 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows;
ykuroda 0:13a5d365ba16 770 MatrixType m_scaledMatrix;
ykuroda 0:13a5d365ba16 771 };
ykuroda 0:13a5d365ba16 772
ykuroda 0:13a5d365ba16 773 template<typename MatrixType, int QRPreconditioner>
ykuroda 0:13a5d365ba16 774 void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, unsigned int computationOptions)
ykuroda 0:13a5d365ba16 775 {
ykuroda 0:13a5d365ba16 776 eigen_assert(rows >= 0 && cols >= 0);
ykuroda 0:13a5d365ba16 777
ykuroda 0:13a5d365ba16 778 if (m_isAllocated &&
ykuroda 0:13a5d365ba16 779 rows == m_rows &&
ykuroda 0:13a5d365ba16 780 cols == m_cols &&
ykuroda 0:13a5d365ba16 781 computationOptions == m_computationOptions)
ykuroda 0:13a5d365ba16 782 {
ykuroda 0:13a5d365ba16 783 return;
ykuroda 0:13a5d365ba16 784 }
ykuroda 0:13a5d365ba16 785
ykuroda 0:13a5d365ba16 786 m_rows = rows;
ykuroda 0:13a5d365ba16 787 m_cols = cols;
ykuroda 0:13a5d365ba16 788 m_isInitialized = false;
ykuroda 0:13a5d365ba16 789 m_isAllocated = true;
ykuroda 0:13a5d365ba16 790 m_computationOptions = computationOptions;
ykuroda 0:13a5d365ba16 791 m_computeFullU = (computationOptions & ComputeFullU) != 0;
ykuroda 0:13a5d365ba16 792 m_computeThinU = (computationOptions & ComputeThinU) != 0;
ykuroda 0:13a5d365ba16 793 m_computeFullV = (computationOptions & ComputeFullV) != 0;
ykuroda 0:13a5d365ba16 794 m_computeThinV = (computationOptions & ComputeThinV) != 0;
ykuroda 0:13a5d365ba16 795 eigen_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U");
ykuroda 0:13a5d365ba16 796 eigen_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V");
ykuroda 0:13a5d365ba16 797 eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
ykuroda 0:13a5d365ba16 798 "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.");
ykuroda 0:13a5d365ba16 799 if (QRPreconditioner == FullPivHouseholderQRPreconditioner)
ykuroda 0:13a5d365ba16 800 {
ykuroda 0:13a5d365ba16 801 eigen_assert(!(m_computeThinU || m_computeThinV) &&
ykuroda 0:13a5d365ba16 802 "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
ykuroda 0:13a5d365ba16 803 "Use the ColPivHouseholderQR preconditioner instead.");
ykuroda 0:13a5d365ba16 804 }
ykuroda 0:13a5d365ba16 805 m_diagSize = (std::min)(m_rows, m_cols);
ykuroda 0:13a5d365ba16 806 m_singularValues.resize(m_diagSize);
ykuroda 0:13a5d365ba16 807 if(RowsAtCompileTime==Dynamic)
ykuroda 0:13a5d365ba16 808 m_matrixU.resize(m_rows, m_computeFullU ? m_rows
ykuroda 0:13a5d365ba16 809 : m_computeThinU ? m_diagSize
ykuroda 0:13a5d365ba16 810 : 0);
ykuroda 0:13a5d365ba16 811 if(ColsAtCompileTime==Dynamic)
ykuroda 0:13a5d365ba16 812 m_matrixV.resize(m_cols, m_computeFullV ? m_cols
ykuroda 0:13a5d365ba16 813 : m_computeThinV ? m_diagSize
ykuroda 0:13a5d365ba16 814 : 0);
ykuroda 0:13a5d365ba16 815 m_workMatrix.resize(m_diagSize, m_diagSize);
ykuroda 0:13a5d365ba16 816
ykuroda 0:13a5d365ba16 817 if(m_cols>m_rows) m_qr_precond_morecols.allocate(*this);
ykuroda 0:13a5d365ba16 818 if(m_rows>m_cols) m_qr_precond_morerows.allocate(*this);
ykuroda 0:13a5d365ba16 819 if(m_rows!=m_cols) m_scaledMatrix.resize(rows,cols);
ykuroda 0:13a5d365ba16 820 }
ykuroda 0:13a5d365ba16 821
ykuroda 0:13a5d365ba16 822 template<typename MatrixType, int QRPreconditioner>
ykuroda 0:13a5d365ba16 823 JacobiSVD<MatrixType, QRPreconditioner>&
ykuroda 0:13a5d365ba16 824 JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsigned int computationOptions)
ykuroda 0:13a5d365ba16 825 {
ykuroda 0:13a5d365ba16 826 check_template_parameters();
ykuroda 0:13a5d365ba16 827
ykuroda 0:13a5d365ba16 828 using std::abs;
ykuroda 0:13a5d365ba16 829 allocate(matrix.rows(), matrix.cols(), computationOptions);
ykuroda 0:13a5d365ba16 830
ykuroda 0:13a5d365ba16 831 // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number of iterations,
ykuroda 0:13a5d365ba16 832 // only worsening the precision of U and V as we accumulate more rotations
ykuroda 0:13a5d365ba16 833 const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon();
ykuroda 0:13a5d365ba16 834
ykuroda 0:13a5d365ba16 835 // limit for very small denormal numbers to be considered zero in order to avoid infinite loops (see bug 286)
ykuroda 0:13a5d365ba16 836 const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min();
ykuroda 0:13a5d365ba16 837
ykuroda 0:13a5d365ba16 838 // Scaling factor to reduce over/under-flows
ykuroda 0:13a5d365ba16 839 RealScalar scale = matrix.cwiseAbs().maxCoeff();
ykuroda 0:13a5d365ba16 840 if(scale==RealScalar(0)) scale = RealScalar(1);
ykuroda 0:13a5d365ba16 841
ykuroda 0:13a5d365ba16 842 /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
ykuroda 0:13a5d365ba16 843
ykuroda 0:13a5d365ba16 844 if(m_rows!=m_cols)
ykuroda 0:13a5d365ba16 845 {
ykuroda 0:13a5d365ba16 846 m_scaledMatrix = matrix / scale;
ykuroda 0:13a5d365ba16 847 m_qr_precond_morecols.run(*this, m_scaledMatrix);
ykuroda 0:13a5d365ba16 848 m_qr_precond_morerows.run(*this, m_scaledMatrix);
ykuroda 0:13a5d365ba16 849 }
ykuroda 0:13a5d365ba16 850 else
ykuroda 0:13a5d365ba16 851 {
ykuroda 0:13a5d365ba16 852 m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize) / scale;
ykuroda 0:13a5d365ba16 853 if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
ykuroda 0:13a5d365ba16 854 if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
ykuroda 0:13a5d365ba16 855 if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
ykuroda 0:13a5d365ba16 856 if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
ykuroda 0:13a5d365ba16 857 }
ykuroda 0:13a5d365ba16 858
ykuroda 0:13a5d365ba16 859 /*** step 2. The main Jacobi SVD iteration. ***/
ykuroda 0:13a5d365ba16 860
ykuroda 0:13a5d365ba16 861 bool finished = false;
ykuroda 0:13a5d365ba16 862 while(!finished)
ykuroda 0:13a5d365ba16 863 {
ykuroda 0:13a5d365ba16 864 finished = true;
ykuroda 0:13a5d365ba16 865
ykuroda 0:13a5d365ba16 866 // do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix
ykuroda 0:13a5d365ba16 867
ykuroda 0:13a5d365ba16 868 for(Index p = 1; p < m_diagSize; ++p)
ykuroda 0:13a5d365ba16 869 {
ykuroda 0:13a5d365ba16 870 for(Index q = 0; q < p; ++q)
ykuroda 0:13a5d365ba16 871 {
ykuroda 0:13a5d365ba16 872 // if this 2x2 sub-matrix is not diagonal already...
ykuroda 0:13a5d365ba16 873 // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't
ykuroda 0:13a5d365ba16 874 // keep us iterating forever. Similarly, small denormal numbers are considered zero.
ykuroda 0:13a5d365ba16 875 using std::max;
ykuroda 0:13a5d365ba16 876 RealScalar threshold = (max)(considerAsZero, precision * (max)(abs(m_workMatrix.coeff(p,p)),
ykuroda 0:13a5d365ba16 877 abs(m_workMatrix.coeff(q,q))));
ykuroda 0:13a5d365ba16 878 // We compare both values to threshold instead of calling max to be robust to NaN (See bug 791)
ykuroda 0:13a5d365ba16 879 if(abs(m_workMatrix.coeff(p,q))>threshold || abs(m_workMatrix.coeff(q,p)) > threshold)
ykuroda 0:13a5d365ba16 880 {
ykuroda 0:13a5d365ba16 881 finished = false;
ykuroda 0:13a5d365ba16 882
ykuroda 0:13a5d365ba16 883 // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal
ykuroda 0:13a5d365ba16 884 internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q);
ykuroda 0:13a5d365ba16 885 JacobiRotation<RealScalar> j_left, j_right;
ykuroda 0:13a5d365ba16 886 internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
ykuroda 0:13a5d365ba16 887
ykuroda 0:13a5d365ba16 888 // accumulate resulting Jacobi rotations
ykuroda 0:13a5d365ba16 889 m_workMatrix.applyOnTheLeft(p,q,j_left);
ykuroda 0:13a5d365ba16 890 if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.transpose());
ykuroda 0:13a5d365ba16 891
ykuroda 0:13a5d365ba16 892 m_workMatrix.applyOnTheRight(p,q,j_right);
ykuroda 0:13a5d365ba16 893 if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right);
ykuroda 0:13a5d365ba16 894 }
ykuroda 0:13a5d365ba16 895 }
ykuroda 0:13a5d365ba16 896 }
ykuroda 0:13a5d365ba16 897 }
ykuroda 0:13a5d365ba16 898
ykuroda 0:13a5d365ba16 899 /*** step 3. The work matrix is now diagonal, so ensure it's positive so its diagonal entries are the singular values ***/
ykuroda 0:13a5d365ba16 900
ykuroda 0:13a5d365ba16 901 for(Index i = 0; i < m_diagSize; ++i)
ykuroda 0:13a5d365ba16 902 {
ykuroda 0:13a5d365ba16 903 RealScalar a = abs(m_workMatrix.coeff(i,i));
ykuroda 0:13a5d365ba16 904 m_singularValues.coeffRef(i) = a;
ykuroda 0:13a5d365ba16 905 if(computeU() && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
ykuroda 0:13a5d365ba16 906 }
ykuroda 0:13a5d365ba16 907
ykuroda 0:13a5d365ba16 908 /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/
ykuroda 0:13a5d365ba16 909
ykuroda 0:13a5d365ba16 910 m_nonzeroSingularValues = m_diagSize;
ykuroda 0:13a5d365ba16 911 for(Index i = 0; i < m_diagSize; i++)
ykuroda 0:13a5d365ba16 912 {
ykuroda 0:13a5d365ba16 913 Index pos;
ykuroda 0:13a5d365ba16 914 RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
ykuroda 0:13a5d365ba16 915 if(maxRemainingSingularValue == RealScalar(0))
ykuroda 0:13a5d365ba16 916 {
ykuroda 0:13a5d365ba16 917 m_nonzeroSingularValues = i;
ykuroda 0:13a5d365ba16 918 break;
ykuroda 0:13a5d365ba16 919 }
ykuroda 0:13a5d365ba16 920 if(pos)
ykuroda 0:13a5d365ba16 921 {
ykuroda 0:13a5d365ba16 922 pos += i;
ykuroda 0:13a5d365ba16 923 std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
ykuroda 0:13a5d365ba16 924 if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i));
ykuroda 0:13a5d365ba16 925 if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i));
ykuroda 0:13a5d365ba16 926 }
ykuroda 0:13a5d365ba16 927 }
ykuroda 0:13a5d365ba16 928
ykuroda 0:13a5d365ba16 929 m_singularValues *= scale;
ykuroda 0:13a5d365ba16 930
ykuroda 0:13a5d365ba16 931 m_isInitialized = true;
ykuroda 0:13a5d365ba16 932 return *this;
ykuroda 0:13a5d365ba16 933 }
ykuroda 0:13a5d365ba16 934
ykuroda 0:13a5d365ba16 935 namespace internal {
ykuroda 0:13a5d365ba16 936 template<typename _MatrixType, int QRPreconditioner, typename Rhs>
ykuroda 0:13a5d365ba16 937 struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
ykuroda 0:13a5d365ba16 938 : solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
ykuroda 0:13a5d365ba16 939 {
ykuroda 0:13a5d365ba16 940 typedef JacobiSVD<_MatrixType, QRPreconditioner> JacobiSVDType;
ykuroda 0:13a5d365ba16 941 EIGEN_MAKE_SOLVE_HELPERS(JacobiSVDType,Rhs)
ykuroda 0:13a5d365ba16 942
ykuroda 0:13a5d365ba16 943 template<typename Dest> void evalTo(Dest& dst) const
ykuroda 0:13a5d365ba16 944 {
ykuroda 0:13a5d365ba16 945 eigen_assert(rhs().rows() == dec().rows());
ykuroda 0:13a5d365ba16 946
ykuroda 0:13a5d365ba16 947 // A = U S V^*
ykuroda 0:13a5d365ba16 948 // So A^{-1} = V S^{-1} U^*
ykuroda 0:13a5d365ba16 949
ykuroda 0:13a5d365ba16 950 Matrix<Scalar, Dynamic, Rhs::ColsAtCompileTime, 0, _MatrixType::MaxRowsAtCompileTime, Rhs::MaxColsAtCompileTime> tmp;
ykuroda 0:13a5d365ba16 951 Index rank = dec().rank();
ykuroda 0:13a5d365ba16 952
ykuroda 0:13a5d365ba16 953 tmp.noalias() = dec().matrixU().leftCols(rank).adjoint() * rhs();
ykuroda 0:13a5d365ba16 954 tmp = dec().singularValues().head(rank).asDiagonal().inverse() * tmp;
ykuroda 0:13a5d365ba16 955 dst = dec().matrixV().leftCols(rank) * tmp;
ykuroda 0:13a5d365ba16 956 }
ykuroda 0:13a5d365ba16 957 };
ykuroda 0:13a5d365ba16 958 } // end namespace internal
ykuroda 0:13a5d365ba16 959
ykuroda 0:13a5d365ba16 960 /** \svd_module
ykuroda 0:13a5d365ba16 961 *
ykuroda 0:13a5d365ba16 962 * \return the singular value decomposition of \c *this computed by two-sided
ykuroda 0:13a5d365ba16 963 * Jacobi transformations.
ykuroda 0:13a5d365ba16 964 *
ykuroda 0:13a5d365ba16 965 * \sa class JacobiSVD
ykuroda 0:13a5d365ba16 966 */
ykuroda 0:13a5d365ba16 967 template<typename Derived>
ykuroda 0:13a5d365ba16 968 JacobiSVD<typename MatrixBase<Derived>::PlainObject>
ykuroda 0:13a5d365ba16 969 MatrixBase<Derived>::jacobiSvd(unsigned int computationOptions) const
ykuroda 0:13a5d365ba16 970 {
ykuroda 0:13a5d365ba16 971 return JacobiSVD<PlainObject>(*this, computationOptions);
ykuroda 0:13a5d365ba16 972 }
ykuroda 0:13a5d365ba16 973
ykuroda 0:13a5d365ba16 974 } // end namespace Eigen
ykuroda 0:13a5d365ba16 975
ykuroda 0:13a5d365ba16 976 #endif // EIGEN_JACOBISVD_H