25 #ifndef EIGEN_JACOBISVD_H
26 #define EIGEN_JACOBISVD_H
33 template<
typename MatrixType,
int QRPreconditioner,
35 struct svd_precondition_2x2_block_to_be_real {};
46 template<
typename MatrixType,
int QRPreconditioner,
int Case>
47 struct qr_preconditioner_should_do_anything
49 enum { a = MatrixType::RowsAtCompileTime !=
Dynamic &&
50 MatrixType::ColsAtCompileTime !=
Dynamic &&
51 MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime,
52 b = MatrixType::RowsAtCompileTime !=
Dynamic &&
53 MatrixType::ColsAtCompileTime !=
Dynamic &&
54 MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime,
61 template<
typename MatrixType,
int QRPreconditioner,
int Case,
62 bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret
63 >
struct qr_preconditioner_impl {};
65 template<
typename MatrixType,
int QRPreconditioner,
int Case>
66 class qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false>
69 typedef typename MatrixType::Index Index;
70 void allocate(
const JacobiSVD<MatrixType, QRPreconditioner>&) {}
71 bool run(JacobiSVD<MatrixType, QRPreconditioner>&,
const MatrixType&)
79 template<
typename MatrixType>
83 typedef typename MatrixType::Index Index;
84 typedef typename MatrixType::Scalar Scalar;
87 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
88 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
90 typedef Matrix<Scalar, 1, RowsAtCompileTime, RowMajor, 1, MaxRowsAtCompileTime> WorkspaceType;
92 void allocate(
const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
94 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
96 m_qr = FullPivHouseholderQR<MatrixType>(svd.rows(), svd.cols());
98 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
101 bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd,
const MatrixType& matrix)
103 if(matrix.rows() > matrix.cols())
105 m_qr.compute(matrix);
106 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).
template triangularView<Upper>();
107 if(svd.m_computeFullU) m_qr.matrixQ().evalTo(svd.m_matrixU, m_workspace);
108 if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
114 FullPivHouseholderQR<MatrixType> m_qr;
115 WorkspaceType m_workspace;
118 template<
typename MatrixType>
122 typedef typename MatrixType::Index Index;
123 typedef typename MatrixType::Scalar Scalar;
126 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
127 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
128 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
129 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
130 Options = MatrixType::Options
132 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
133 TransposeTypeWithSameStorageOrder;
135 void allocate(
const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
137 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
139 m_qr = FullPivHouseholderQR<TransposeTypeWithSameStorageOrder>(svd.cols(), svd.rows());
141 m_adjoint.resize(svd.cols(), svd.rows());
142 if (svd.m_computeFullV) m_workspace.resize(svd.cols());
145 bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd,
const MatrixType& matrix)
147 if(matrix.cols() > matrix.rows())
149 m_adjoint = matrix.adjoint();
150 m_qr.compute(m_adjoint);
151 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).
template triangularView<Upper>().adjoint();
152 if(svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace);
153 if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
159 FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
160 TransposeTypeWithSameStorageOrder m_adjoint;
161 typename internal::plain_row_type<MatrixType>::type m_workspace;
166 template<
typename MatrixType>
170 typedef typename MatrixType::Index Index;
172 void allocate(
const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
174 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
176 m_qr = ColPivHouseholderQR<MatrixType>(svd.rows(), svd.cols());
178 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
179 else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
182 bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd,
const MatrixType& matrix)
184 if(matrix.rows() > matrix.cols())
186 m_qr.compute(matrix);
187 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).
template triangularView<Upper>();
188 if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
189 else if(svd.m_computeThinU)
191 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
192 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
194 if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
201 ColPivHouseholderQR<MatrixType> m_qr;
202 typename internal::plain_col_type<MatrixType>::type m_workspace;
205 template<
typename MatrixType>
209 typedef typename MatrixType::Index Index;
210 typedef typename MatrixType::Scalar Scalar;
213 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
214 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
215 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
216 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
217 Options = MatrixType::Options
220 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
221 TransposeTypeWithSameStorageOrder;
223 void allocate(
const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
225 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
227 m_qr = ColPivHouseholderQR<TransposeTypeWithSameStorageOrder>(svd.cols(), svd.rows());
229 if (svd.m_computeFullV) m_workspace.resize(svd.cols());
230 else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
231 m_adjoint.resize(svd.cols(), svd.rows());
234 bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd,
const MatrixType& matrix)
236 if(matrix.cols() > matrix.rows())
238 m_adjoint = matrix.adjoint();
239 m_qr.compute(m_adjoint);
241 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).
template triangularView<Upper>().adjoint();
242 if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
243 else if(svd.m_computeThinV)
245 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
246 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
248 if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
255 ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
256 TransposeTypeWithSameStorageOrder m_adjoint;
257 typename internal::plain_row_type<MatrixType>::type m_workspace;
262 template<
typename MatrixType>
266 typedef typename MatrixType::Index Index;
268 void allocate(
const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
270 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
272 m_qr = HouseholderQR<MatrixType>(svd.rows(), svd.cols());
274 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
275 else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
278 bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd,
const MatrixType& matrix)
280 if(matrix.rows() > matrix.cols())
282 m_qr.compute(matrix);
283 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).
template triangularView<Upper>();
284 if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
285 else if(svd.m_computeThinU)
287 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
288 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
290 if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
296 HouseholderQR<MatrixType> m_qr;
297 typename internal::plain_col_type<MatrixType>::type m_workspace;
300 template<
typename MatrixType>
304 typedef typename MatrixType::Index Index;
305 typedef typename MatrixType::Scalar Scalar;
308 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
309 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
310 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
311 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
312 Options = MatrixType::Options
315 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
316 TransposeTypeWithSameStorageOrder;
318 void allocate(
const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
320 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
322 m_qr = HouseholderQR<TransposeTypeWithSameStorageOrder>(svd.cols(), svd.rows());
324 if (svd.m_computeFullV) m_workspace.resize(svd.cols());
325 else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
326 m_adjoint.resize(svd.cols(), svd.rows());
329 bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd,
const MatrixType& matrix)
331 if(matrix.cols() > matrix.rows())
333 m_adjoint = matrix.adjoint();
334 m_qr.compute(m_adjoint);
336 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).
template triangularView<Upper>().adjoint();
337 if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
338 else if(svd.m_computeThinV)
340 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
341 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
343 if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
350 HouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
351 TransposeTypeWithSameStorageOrder m_adjoint;
352 typename internal::plain_row_type<MatrixType>::type m_workspace;
360 template<
typename MatrixType,
int QRPreconditioner>
361 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false>
363 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
364 typedef typename SVD::Index Index;
365 static void run(
typename SVD::WorkMatrixType&, SVD&, Index, Index) {}
368 template<
typename MatrixType,
int QRPreconditioner>
369 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
371 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
372 typedef typename MatrixType::Scalar Scalar;
373 typedef typename MatrixType::RealScalar RealScalar;
374 typedef typename SVD::Index Index;
375 static void run(
typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index
p, Index
q)
378 JacobiRotation<Scalar> rot;
379 RealScalar n =
sqrt(
abs2(work_matrix.coeff(p,p)) +
abs2(work_matrix.coeff(q,p)));
382 z =
abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
383 work_matrix.row(p) *= z;
384 if(svd.computeU()) svd.m_matrixU.col(p) *=
conj(z);
385 z =
abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
386 work_matrix.row(q) *= z;
387 if(svd.computeU()) svd.m_matrixU.col(q) *=
conj(z);
391 rot.c() =
conj(work_matrix.coeff(p,p)) / n;
392 rot.s() = work_matrix.coeff(q,p) / n;
393 work_matrix.applyOnTheLeft(p,q,rot);
394 if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint());
395 if(work_matrix.coeff(p,q) != Scalar(0))
397 Scalar z =
abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
398 work_matrix.col(q) *= z;
399 if(svd.computeV()) svd.m_matrixV.col(q) *= z;
401 if(work_matrix.coeff(q,q) != Scalar(0))
403 z =
abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
404 work_matrix.row(q) *= z;
405 if(svd.computeU()) svd.m_matrixU.col(q) *=
conj(z);
411 template<
typename MatrixType,
typename RealScalar,
typename Index>
417 m <<
real(matrix.coeff(p,p)),
real(matrix.coeff(p,q)),
418 real(matrix.coeff(q,p)),
real(matrix.coeff(q,q));
422 if(t == RealScalar(0))
424 rot1.c() = RealScalar(0);
425 rot1.s() = d > RealScalar(0) ? RealScalar(1) : RealScalar(-1);
429 RealScalar u = d / t;
430 rot1.c() = RealScalar(1) /
sqrt(RealScalar(1) +
abs2(u));
431 rot1.s() = rot1.c() * u;
493 template<
typename _MatrixType,
int QRPreconditioner>
class JacobiSVD
498 typedef typename MatrixType::Scalar
Scalar;
500 typedef typename MatrixType::Index
Index;
518 typedef typename internal::plain_row_type<MatrixType>::type
RowType;
519 typedef typename internal::plain_col_type<MatrixType>::type
ColType;
549 allocate(rows, cols, computationOptions);
568 compute(matrix, computationOptions);
606 eigen_assert(
computeU() &&
"This JacobiSVD decomposition didn't compute U. Did you ask for it?");
622 eigen_assert(
computeV() &&
"This JacobiSVD decomposition didn't compute V. Did you ask for it?");
651 template<
typename Rhs>
652 inline const internal::solve_retval<JacobiSVD, Rhs>
656 eigen_assert(
computeU() &&
computeV() &&
"JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice).");
657 return internal::solve_retval<JacobiSVD, Rhs>(*
this, b.derived());
684 template<
typename __MatrixType,
int _QRPreconditioner,
bool _IsComplex>
685 friend struct internal::svd_precondition_2x2_block_to_be_real;
686 template<
typename __MatrixType,
int _QRPreconditioner,
int _Case,
bool _DoAnything>
687 friend struct internal::qr_preconditioner_impl;
689 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows>
m_qr_precond_morecols;
690 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols>
m_qr_precond_morerows;
693 template<
typename MatrixType,
int QRPreconditioner>
701 computationOptions == m_computationOptions)
708 m_isInitialized =
false;
709 m_isAllocated =
true;
710 m_computationOptions = computationOptions;
711 m_computeFullU = (computationOptions &
ComputeFullU) != 0;
712 m_computeThinU = (computationOptions &
ComputeThinU) != 0;
713 m_computeFullV = (computationOptions &
ComputeFullV) != 0;
714 m_computeThinV = (computationOptions &
ComputeThinV) != 0;
715 eigen_assert(!(m_computeFullU && m_computeThinU) &&
"JacobiSVD: you can't ask for both full and thin U");
716 eigen_assert(!(m_computeFullV && m_computeThinV) &&
"JacobiSVD: you can't ask for both full and thin V");
718 "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.");
722 "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
723 "Use the ColPivHouseholderQR preconditioner instead.");
725 m_diagSize = (std::min)(m_rows, m_cols);
726 m_singularValues.resize(m_diagSize);
727 m_matrixU.resize(m_rows, m_computeFullU ? m_rows
728 : m_computeThinU ? m_diagSize
730 m_matrixV.resize(m_cols, m_computeFullV ? m_cols
731 : m_computeThinV ? m_diagSize
733 m_workMatrix.resize(m_diagSize, m_diagSize);
735 if(m_cols>m_rows) m_qr_precond_morecols.allocate(*
this);
736 if(m_rows>m_cols) m_qr_precond_morerows.allocate(*
this);
739 template<
typename MatrixType,
int QRPreconditioner>
740 JacobiSVD<MatrixType, QRPreconditioner>&
743 allocate(matrix.rows(), matrix.cols(), computationOptions);
750 const RealScalar considerAsZero =
RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min();
754 if(!m_qr_precond_morecols.run(*
this, matrix) && !m_qr_precond_morerows.run(*
this, matrix))
756 m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize);
757 if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
758 if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
759 if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
760 if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
765 bool finished =
false;
772 for(
Index p = 1; p < m_diagSize; ++
p)
787 internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *
this, p, q);
792 m_workMatrix.applyOnTheLeft(p,q,j_left);
793 if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.
transpose());
795 m_workMatrix.applyOnTheRight(p,q,j_right);
796 if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right);
804 for(
Index i = 0; i < m_diagSize; ++i)
807 m_singularValues.coeffRef(i) = a;
808 if(computeU() && (a!=
RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
813 m_nonzeroSingularValues = m_diagSize;
814 for(
Index i = 0; i < m_diagSize; i++)
817 RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
818 if(maxRemainingSingularValue ==
RealScalar(0))
820 m_nonzeroSingularValues = i;
826 std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
827 if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i));
828 if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i));
832 m_isInitialized =
true;
837 template<
typename _MatrixType,
int QRPreconditioner,
typename Rhs>
838 struct solve_retval<
JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
839 : solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
844 template<typename Dest>
void evalTo(Dest& dst)
const
851 Index diagSize = (std::min)(dec().rows(), dec().cols());
852 typename JacobiSVDType::SingularValuesType invertedSingVals(diagSize);
854 Index nonzeroSingVals = dec().nonzeroSingularValues();
855 invertedSingVals.head(nonzeroSingVals) = dec().singularValues().head(nonzeroSingVals).array().inverse();
856 invertedSingVals.tail(diagSize - nonzeroSingVals).setZero();
858 dst = dec().matrixV().leftCols(diagSize)
859 * invertedSingVals.asDiagonal()
860 * dec().matrixU().leftCols(diagSize).adjoint()
873 template<
typename Derived>
874 JacobiSVD<typename MatrixBase<Derived>::PlainObject>
882 #endif // EIGEN_JACOBISVD_H