| // This file is part of Eigen, a lightweight C++ template library | |
| // for linear algebra. | |
| // | |
| // Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com> | |
| // Copyright (C) 2013-2014 Gael Guennebaud <gael.guennebaud@inria.fr> | |
| // | |
| // This Source Code Form is subject to the terms of the Mozilla | |
| // Public License v. 2.0. If a copy of the MPL was not distributed | |
| // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. | |
| // SPDX-License-Identifier: MPL-2.0 | |
| #ifndef EIGEN_JACOBISVD_H | |
| #define EIGEN_JACOBISVD_H | |
| // IWYU pragma: private | |
| #include "./InternalHeaderCheck.h" | |
| namespace Eigen { | |
| namespace internal { | |
| // forward declaration (needed by ICC) | |
| // the empty body is required by MSVC | |
| template <typename MatrixType, int Options, bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex> | |
| struct svd_precondition_2x2_block_to_be_real {}; | |
| /*** QR preconditioners (R-SVD) | |
| *** | |
| *** Their role is to reduce the problem of computing the SVD to the case of a square matrix. | |
| *** This approach, known as R-SVD, is an optimization for rectangular-enough matrices, and is a requirement for | |
| *** JacobiSVD which by itself is only able to work on square matrices. | |
| ***/ | |
| enum { PreconditionIfMoreColsThanRows, PreconditionIfMoreRowsThanCols }; | |
| template <typename MatrixType, int QRPreconditioner, int Case> | |
| struct qr_preconditioner_should_do_anything | |
| : bool_constant<!((QRPreconditioner == NoQRPreconditioner) || | |
| (Case == PreconditionIfMoreColsThanRows && MatrixType::RowsAtCompileTime != Dynamic && | |
| MatrixType::ColsAtCompileTime != Dynamic && | |
| MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime) || | |
| (Case == PreconditionIfMoreRowsThanCols && MatrixType::RowsAtCompileTime != Dynamic && | |
| MatrixType::ColsAtCompileTime != Dynamic && | |
| MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime))> {}; | |
| template <typename MatrixType, int Options, int QRPreconditioner, int Case, | |
| bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::value> | |
| struct qr_preconditioner_impl {}; | |
| template <typename MatrixType, int Options, int QRPreconditioner, int Case> | |
| class qr_preconditioner_impl<MatrixType, Options, QRPreconditioner, Case, false> { | |
| public: | |
| void allocate(const JacobiSVD<MatrixType, Options>&) {} | |
| template <typename Xpr> | |
| bool run(JacobiSVD<MatrixType, Options>&, const Xpr&) { | |
| return false; | |
| } | |
| }; | |
| /*** preconditioner using FullPivHouseholderQR ***/ | |
| template <typename MatrixType, int Options> | |
| class qr_preconditioner_impl<MatrixType, Options, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, | |
| true> { | |
| public: | |
| using Scalar = typename MatrixType::Scalar; | |
| using SVDType = JacobiSVD<MatrixType, Options>; | |
| enum { WorkspaceSize = MatrixType::RowsAtCompileTime, MaxWorkspaceSize = MatrixType::MaxRowsAtCompileTime }; | |
| using WorkspaceType = Matrix<Scalar, 1, WorkspaceSize, RowMajor, 1, MaxWorkspaceSize>; | |
| void allocate(const SVDType& svd) { | |
| if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) { | |
| internal::destroy_at(&m_qr); | |
| internal::construct_at(&m_qr, svd.rows(), svd.cols()); | |
| } | |
| if (svd.m_computeFullU) m_workspace.resize(svd.rows()); | |
| } | |
| template <typename Xpr> | |
| bool run(SVDType& svd, const Xpr& matrix) { | |
| if (matrix.rows() > matrix.cols()) { | |
| m_qr.compute(matrix); | |
| svd.m_workMatrix = m_qr.matrixQR().block(0, 0, matrix.cols(), matrix.cols()).template triangularView<Upper>(); | |
| if (svd.m_computeFullU) m_qr.matrixQ().evalTo(svd.m_matrixU, m_workspace); | |
| if (svd.computeV()) svd.m_matrixV = m_qr.colsPermutation(); | |
| return true; | |
| } | |
| return false; | |
| } | |
| private: | |
| using QRType = FullPivHouseholderQR<MatrixType>; | |
| QRType m_qr; | |
| WorkspaceType m_workspace; | |
| }; | |
| template <typename MatrixType, int Options> | |
| class qr_preconditioner_impl<MatrixType, Options, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, | |
| true> { | |
| public: | |
| using Scalar = typename MatrixType::Scalar; | |
| using SVDType = JacobiSVD<MatrixType, Options>; | |
| enum { | |
| RowsAtCompileTime = MatrixType::RowsAtCompileTime, | |
| ColsAtCompileTime = MatrixType::ColsAtCompileTime, | |
| MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, | |
| MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, | |
| MatrixOptions = traits<MatrixType>::Options | |
| }; | |
| using TransposeTypeWithSameStorageOrder = | |
| typename internal::make_proper_matrix_type<Scalar, ColsAtCompileTime, RowsAtCompileTime, MatrixOptions, | |
| MaxColsAtCompileTime, MaxRowsAtCompileTime>::type; | |
| void allocate(const SVDType& svd) { | |
| if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) { | |
| internal::destroy_at(&m_qr); | |
| internal::construct_at(&m_qr, svd.cols(), svd.rows()); | |
| } | |
| if (svd.m_computeFullV) m_workspace.resize(svd.cols()); | |
| } | |
| template <typename Xpr> | |
| bool run(SVDType& svd, const Xpr& matrix) { | |
| if (matrix.cols() > matrix.rows()) { | |
| m_qr.compute(matrix.adjoint()); | |
| svd.m_workMatrix = | |
| m_qr.matrixQR().block(0, 0, matrix.rows(), matrix.rows()).template triangularView<Upper>().adjoint(); | |
| if (svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace); | |
| if (svd.computeU()) svd.m_matrixU = m_qr.colsPermutation(); | |
| return true; | |
| } else | |
| return false; | |
| } | |
| private: | |
| using QRType = FullPivHouseholderQR<TransposeTypeWithSameStorageOrder>; | |
| QRType m_qr; | |
| typename plain_row_type<MatrixType>::type m_workspace; | |
| }; | |
| /*** preconditioner using ColPivHouseholderQR ***/ | |
| template <typename MatrixType, int Options> | |
| class qr_preconditioner_impl<MatrixType, Options, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, | |
| true> { | |
| public: | |
| using Scalar = typename MatrixType::Scalar; | |
| using SVDType = JacobiSVD<MatrixType, Options>; | |
| enum { | |
| WorkspaceSize = internal::traits<SVDType>::MatrixUColsAtCompileTime, | |
| MaxWorkspaceSize = internal::traits<SVDType>::MatrixUMaxColsAtCompileTime | |
| }; | |
| using WorkspaceType = Matrix<Scalar, 1, WorkspaceSize, RowMajor, 1, MaxWorkspaceSize>; | |
| void allocate(const SVDType& svd) { | |
| if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) { | |
| internal::destroy_at(&m_qr); | |
| internal::construct_at(&m_qr, svd.rows(), svd.cols()); | |
| } | |
| if (svd.m_computeFullU) | |
| m_workspace.resize(svd.rows()); | |
| else if (svd.m_computeThinU) | |
| m_workspace.resize(svd.cols()); | |
| } | |
| template <typename Xpr> | |
| bool run(SVDType& svd, const Xpr& matrix) { | |
| if (matrix.rows() > matrix.cols()) { | |
| m_qr.compute(matrix); | |
| svd.m_workMatrix = m_qr.matrixQR().block(0, 0, matrix.cols(), matrix.cols()).template triangularView<Upper>(); | |
| if (svd.m_computeFullU) | |
| m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace); | |
| else if (svd.m_computeThinU) { | |
| svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols()); | |
| m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace); | |
| } | |
| if (svd.computeV()) svd.m_matrixV = m_qr.colsPermutation(); | |
| return true; | |
| } | |
| return false; | |
| } | |
| private: | |
| using QRType = ColPivHouseholderQR<MatrixType>; | |
| QRType m_qr; | |
| WorkspaceType m_workspace; | |
| }; | |
| template <typename MatrixType, int Options> | |
| class qr_preconditioner_impl<MatrixType, Options, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, | |
| true> { | |
| public: | |
| using Scalar = typename MatrixType::Scalar; | |
| using SVDType = JacobiSVD<MatrixType, Options>; | |
| enum { | |
| RowsAtCompileTime = MatrixType::RowsAtCompileTime, | |
| ColsAtCompileTime = MatrixType::ColsAtCompileTime, | |
| MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, | |
| MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, | |
| MatrixOptions = internal::traits<MatrixType>::Options, | |
| WorkspaceSize = internal::traits<SVDType>::MatrixVColsAtCompileTime, | |
| MaxWorkspaceSize = internal::traits<SVDType>::MatrixVMaxColsAtCompileTime | |
| }; | |
| using WorkspaceType = Matrix<Scalar, WorkspaceSize, 1, ColMajor, MaxWorkspaceSize, 1>; | |
| using TransposeTypeWithSameStorageOrder = | |
| typename internal::make_proper_matrix_type<Scalar, ColsAtCompileTime, RowsAtCompileTime, MatrixOptions, | |
| MaxColsAtCompileTime, MaxRowsAtCompileTime>::type; | |
| void allocate(const SVDType& svd) { | |
| if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) { | |
| internal::destroy_at(&m_qr); | |
| internal::construct_at(&m_qr, svd.cols(), svd.rows()); | |
| } | |
| if (svd.m_computeFullV) | |
| m_workspace.resize(svd.cols()); | |
| else if (svd.m_computeThinV) | |
| m_workspace.resize(svd.rows()); | |
| } | |
| template <typename Xpr> | |
| bool run(SVDType& svd, const Xpr& matrix) { | |
| if (matrix.cols() > matrix.rows()) { | |
| m_qr.compute(matrix.adjoint()); | |
| svd.m_workMatrix = | |
| m_qr.matrixQR().block(0, 0, matrix.rows(), matrix.rows()).template triangularView<Upper>().adjoint(); | |
| if (svd.m_computeFullV) | |
| m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace); | |
| else if (svd.m_computeThinV) { | |
| svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows()); | |
| m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace); | |
| } | |
| if (svd.computeU()) svd.m_matrixU = m_qr.colsPermutation(); | |
| return true; | |
| } else | |
| return false; | |
| } | |
| private: | |
| using QRType = ColPivHouseholderQR<TransposeTypeWithSameStorageOrder>; | |
| QRType m_qr; | |
| WorkspaceType m_workspace; | |
| }; | |
| /*** preconditioner using HouseholderQR ***/ | |
| template <typename MatrixType, int Options> | |
| class qr_preconditioner_impl<MatrixType, Options, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> { | |
| public: | |
| using Scalar = typename MatrixType::Scalar; | |
| using SVDType = JacobiSVD<MatrixType, Options>; | |
| enum { | |
| WorkspaceSize = internal::traits<SVDType>::MatrixUColsAtCompileTime, | |
| MaxWorkspaceSize = internal::traits<SVDType>::MatrixUMaxColsAtCompileTime | |
| }; | |
| using WorkspaceType = Matrix<Scalar, 1, WorkspaceSize, RowMajor, 1, MaxWorkspaceSize>; | |
| void allocate(const SVDType& svd) { | |
| if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) { | |
| internal::destroy_at(&m_qr); | |
| internal::construct_at(&m_qr, svd.rows(), svd.cols()); | |
| } | |
| if (svd.m_computeFullU) | |
| m_workspace.resize(svd.rows()); | |
| else if (svd.m_computeThinU) | |
| m_workspace.resize(svd.cols()); | |
| } | |
| template <typename Xpr> | |
| bool run(SVDType& svd, const Xpr& matrix) { | |
| if (matrix.rows() > matrix.cols()) { | |
| m_qr.compute(matrix); | |
| svd.m_workMatrix = m_qr.matrixQR().block(0, 0, matrix.cols(), matrix.cols()).template triangularView<Upper>(); | |
| if (svd.m_computeFullU) | |
| m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace); | |
| else if (svd.m_computeThinU) { | |
| svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols()); | |
| m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace); | |
| } | |
| if (svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols()); | |
| return true; | |
| } | |
| return false; | |
| } | |
| private: | |
| using QRType = HouseholderQR<MatrixType>; | |
| QRType m_qr; | |
| WorkspaceType m_workspace; | |
| }; | |
| template <typename MatrixType, int Options> | |
| class qr_preconditioner_impl<MatrixType, Options, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> { | |
| public: | |
| using Scalar = typename MatrixType::Scalar; | |
| using SVDType = JacobiSVD<MatrixType, Options>; | |
| enum { | |
| RowsAtCompileTime = MatrixType::RowsAtCompileTime, | |
| ColsAtCompileTime = MatrixType::ColsAtCompileTime, | |
| MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, | |
| MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, | |
| MatrixOptions = internal::traits<MatrixType>::Options, | |
| WorkspaceSize = internal::traits<SVDType>::MatrixVColsAtCompileTime, | |
| MaxWorkspaceSize = internal::traits<SVDType>::MatrixVMaxColsAtCompileTime | |
| }; | |
| using WorkspaceType = Matrix<Scalar, WorkspaceSize, 1, ColMajor, MaxWorkspaceSize, 1>; | |
| using TransposeTypeWithSameStorageOrder = | |
| typename internal::make_proper_matrix_type<Scalar, ColsAtCompileTime, RowsAtCompileTime, MatrixOptions, | |
| MaxColsAtCompileTime, MaxRowsAtCompileTime>::type; | |
| void allocate(const SVDType& svd) { | |
| if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) { | |
| internal::destroy_at(&m_qr); | |
| internal::construct_at(&m_qr, svd.cols(), svd.rows()); | |
| } | |
| if (svd.m_computeFullV) | |
| m_workspace.resize(svd.cols()); | |
| else if (svd.m_computeThinV) | |
| m_workspace.resize(svd.rows()); | |
| } | |
| template <typename Xpr> | |
| bool run(SVDType& svd, const Xpr& matrix) { | |
| if (matrix.cols() > matrix.rows()) { | |
| m_qr.compute(matrix.adjoint()); | |
| svd.m_workMatrix = | |
| m_qr.matrixQR().block(0, 0, matrix.rows(), matrix.rows()).template triangularView<Upper>().adjoint(); | |
| if (svd.m_computeFullV) | |
| m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace); | |
| else if (svd.m_computeThinV) { | |
| svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows()); | |
| m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace); | |
| } | |
| if (svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows()); | |
| return true; | |
| } else | |
| return false; | |
| } | |
| private: | |
| using QRType = HouseholderQR<TransposeTypeWithSameStorageOrder>; | |
| QRType m_qr; | |
| WorkspaceType m_workspace; | |
| }; | |
| /*** 2x2 SVD implementation | |
| *** | |
| *** JacobiSVD consists in performing a series of 2x2 SVD subproblems | |
| ***/ | |
| template <typename MatrixType, int Options> | |
| struct svd_precondition_2x2_block_to_be_real<MatrixType, Options, false> { | |
| using SVD = JacobiSVD<MatrixType, Options>; | |
| using RealScalar = typename MatrixType::RealScalar; | |
| static bool run(typename SVD::WorkMatrixType&, SVD&, Index, Index, RealScalar&) { return true; } | |
| }; | |
| template <typename MatrixType, int Options> | |
| struct svd_precondition_2x2_block_to_be_real<MatrixType, Options, true> { | |
| using SVD = JacobiSVD<MatrixType, Options>; | |
| using Scalar = typename MatrixType::Scalar; | |
| using RealScalar = typename MatrixType::RealScalar; | |
| static bool run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q, RealScalar& maxDiagEntry) { | |
| using numext::abs; | |
| using numext::sqrt; | |
| Scalar z; | |
| JacobiRotation<Scalar> rot; | |
| RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p, p)) + numext::abs2(work_matrix.coeff(q, p))); | |
| const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)(); | |
| const RealScalar precision = NumTraits<Scalar>::epsilon(); | |
| if (numext::is_exactly_zero(n)) { | |
| // make sure first column is zero | |
| work_matrix.coeffRef(p, p) = work_matrix.coeffRef(q, p) = Scalar(0); | |
| if (abs(numext::imag(work_matrix.coeff(p, q))) > considerAsZero) { | |
| // work_matrix.coeff(p,q) can be zero if work_matrix.coeff(q,p) is not zero but small enough to underflow when | |
| // computing n | |
| z = abs(work_matrix.coeff(p, q)) / work_matrix.coeff(p, q); | |
| work_matrix.row(p) *= z; | |
| if (svd.computeU()) svd.m_matrixU.col(p) *= conj(z); | |
| } | |
| if (abs(numext::imag(work_matrix.coeff(q, q))) > considerAsZero) { | |
| z = abs(work_matrix.coeff(q, q)) / work_matrix.coeff(q, q); | |
| work_matrix.row(q) *= z; | |
| if (svd.computeU()) svd.m_matrixU.col(q) *= conj(z); | |
| } | |
| // otherwise the second row is already zero, so we have nothing to do. | |
| } else { | |
| rot.c() = conj(work_matrix.coeff(p, p)) / n; | |
| rot.s() = work_matrix.coeff(q, p) / n; | |
| work_matrix.applyOnTheLeft(p, q, rot); | |
| if (svd.computeU()) svd.m_matrixU.applyOnTheRight(p, q, rot.adjoint()); | |
| if (abs(numext::imag(work_matrix.coeff(p, q))) > considerAsZero) { | |
| z = abs(work_matrix.coeff(p, q)) / work_matrix.coeff(p, q); | |
| work_matrix.col(q) *= z; | |
| if (svd.computeV()) svd.m_matrixV.col(q) *= z; | |
| } | |
| if (abs(numext::imag(work_matrix.coeff(q, q))) > considerAsZero) { | |
| z = abs(work_matrix.coeff(q, q)) / work_matrix.coeff(q, q); | |
| work_matrix.row(q) *= z; | |
| if (svd.computeU()) svd.m_matrixU.col(q) *= conj(z); | |
| } | |
| } | |
| // update largest diagonal entry | |
| maxDiagEntry = numext::maxi<RealScalar>( | |
| maxDiagEntry, numext::maxi<RealScalar>(abs(work_matrix.coeff(p, p)), abs(work_matrix.coeff(q, q)))); | |
| // and check whether the 2x2 block is already diagonal | |
| RealScalar threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry); | |
| return abs(work_matrix.coeff(p, q)) > threshold || abs(work_matrix.coeff(q, p)) > threshold; | |
| } | |
| }; | |
| template <typename WorkMatrixType, typename MatrixUType, typename MatrixVType> | |
| EIGEN_DONT_INLINE bool jacobi_svd_nonblocking_sweep(WorkMatrixType& work_matrix, MatrixUType& matrix_u, | |
| MatrixVType& matrix_v, bool compute_u, bool compute_v, | |
| typename WorkMatrixType::RealScalar considerAsZero, | |
| typename WorkMatrixType::RealScalar precision, | |
| typename WorkMatrixType::RealScalar& maxDiagEntry) { | |
| using numext::abs; | |
| using numext::sqrt; | |
| using Scalar = typename WorkMatrixType::Scalar; | |
| using RealScalar = typename WorkMatrixType::RealScalar; | |
| const Index n = work_matrix.rows(); | |
| bool notFinished = false; | |
| EIGEN_IF_CONSTEXPR (NumTraits<Scalar>::IsComplex) { | |
| // Complex non-blocking sweep: condition each 2x2 block to be real before diagonalizing. | |
| for (Index p = 1; p < n; ++p) { | |
| for (Index q = 0; q < p; ++q) { | |
| RealScalar threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry); | |
| if (abs(work_matrix.coeff(p, q)) > threshold || abs(work_matrix.coeff(q, p)) > threshold) { | |
| notFinished = true; | |
| Scalar z; | |
| bool doRealSvd = true; | |
| RealScalar nn = sqrt(numext::abs2(work_matrix.coeff(p, p)) + numext::abs2(work_matrix.coeff(q, p))); | |
| if (numext::is_exactly_zero(nn)) { | |
| // Make sure first column is zero. | |
| work_matrix.coeffRef(p, p) = work_matrix.coeffRef(q, p) = Scalar(0); | |
| if (abs(numext::imag(work_matrix.coeff(p, q))) > considerAsZero) { | |
| // work_matrix.coeff(p, q) can be zero if work_matrix.coeff(q, p) is not zero but small enough to | |
| // underflow when computing nn. | |
| z = abs(work_matrix.coeff(p, q)) / work_matrix.coeff(p, q); | |
| work_matrix.row(p) *= z; | |
| if (compute_u) matrix_u.col(p) *= numext::conj(z); | |
| } | |
| if (abs(numext::imag(work_matrix.coeff(q, q))) > considerAsZero) { | |
| z = abs(work_matrix.coeff(q, q)) / work_matrix.coeff(q, q); | |
| work_matrix.row(q) *= z; | |
| if (compute_u) matrix_u.col(q) *= numext::conj(z); | |
| } | |
| } else { | |
| JacobiRotation<Scalar> rot; | |
| rot.c() = numext::conj(work_matrix.coeff(p, p)) / nn; | |
| rot.s() = work_matrix.coeff(q, p) / nn; | |
| work_matrix.applyOnTheLeft(p, q, rot); | |
| if (compute_u) matrix_u.applyOnTheRight(p, q, rot.adjoint()); | |
| if (abs(numext::imag(work_matrix.coeff(p, q))) > considerAsZero) { | |
| z = abs(work_matrix.coeff(p, q)) / work_matrix.coeff(p, q); | |
| work_matrix.col(q) *= z; | |
| if (compute_v) matrix_v.col(q) *= z; | |
| } | |
| if (abs(numext::imag(work_matrix.coeff(q, q))) > considerAsZero) { | |
| z = abs(work_matrix.coeff(q, q)) / work_matrix.coeff(q, q); | |
| work_matrix.row(q) *= z; | |
| if (compute_u) matrix_u.col(q) *= numext::conj(z); | |
| } | |
| } | |
| maxDiagEntry = numext::maxi<RealScalar>( | |
| maxDiagEntry, numext::maxi<RealScalar>(abs(work_matrix.coeff(p, p)), abs(work_matrix.coeff(q, q)))); | |
| threshold = numext::maxi<RealScalar>(considerAsZero, NumTraits<Scalar>::epsilon() * maxDiagEntry); | |
| doRealSvd = abs(work_matrix.coeff(p, q)) > threshold || abs(work_matrix.coeff(q, p)) > threshold; | |
| if (doRealSvd) { | |
| JacobiRotation<RealScalar> j_left, j_right; | |
| internal::real_2x2_jacobi_svd(work_matrix, p, q, &j_left, &j_right); | |
| work_matrix.applyOnTheLeft(p, q, j_left); | |
| if (compute_u) matrix_u.applyOnTheRight(p, q, j_left.transpose()); | |
| work_matrix.applyOnTheRight(p, q, j_right); | |
| if (compute_v) matrix_v.applyOnTheRight(p, q, j_right); | |
| maxDiagEntry = numext::maxi<RealScalar>( | |
| maxDiagEntry, numext::maxi<RealScalar>(abs(work_matrix.coeff(p, p)), abs(work_matrix.coeff(q, q)))); | |
| } | |
| } | |
| } | |
| } | |
| } else { | |
| // Real non-blocking sweep: diagonalize each 2x2 block directly. | |
| RealScalar threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry); | |
| for (Index p = 1; p < n; ++p) { | |
| for (Index q = 0; q < p; ++q) { | |
| if (abs(work_matrix.coeff(p, q)) > threshold || abs(work_matrix.coeff(q, p)) > threshold) { | |
| notFinished = true; | |
| JacobiRotation<RealScalar> j_left, j_right; | |
| internal::real_2x2_jacobi_svd(work_matrix, p, q, &j_left, &j_right); | |
| work_matrix.applyOnTheLeft(p, q, j_left); | |
| if (compute_u) matrix_u.applyOnTheRight(p, q, j_left.transpose()); | |
| work_matrix.applyOnTheRight(p, q, j_right); | |
| if (compute_v) matrix_v.applyOnTheRight(p, q, j_right); | |
| maxDiagEntry = numext::maxi<RealScalar>( | |
| maxDiagEntry, numext::maxi<RealScalar>(abs(work_matrix.coeff(p, p)), abs(work_matrix.coeff(q, q)))); | |
| threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry); | |
| } | |
| } | |
| } | |
| } | |
| return notFinished; | |
| } | |
| template <typename MatrixType_, int Options> | |
| struct traits<JacobiSVD<MatrixType_, Options>> : svd_traits<MatrixType_, Options> { | |
| using MatrixType = MatrixType_; | |
| }; | |
| } // end namespace internal | |
| /** \ingroup SVD_Module | |
| * | |
| * | |
| * \class JacobiSVD | |
| * | |
| * \brief Two-sided Jacobi SVD decomposition of a rectangular matrix | |
| * | |
| * \tparam MatrixType_ the type of the matrix of which we are computing the SVD decomposition | |
| * \tparam Options this optional parameter allows one to specify the type of QR decomposition that will be used | |
| * internally for the R-SVD step for non-square matrices. Additionally, it allows one to specify whether to compute thin | |
| * or full unitaries \a U and \a V. See discussion of possible values below. | |
| * | |
| * SVD decomposition consists in decomposing any n-by-p matrix \a A as a product | |
| * \f[ A = U S V^* \f] | |
| * 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; 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 and right \em singular \em vectors of \a A respectively. | |
| * | |
| * Singular values are always sorted in decreasing order. | |
| * | |
| * 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. | |
| * | |
| * 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 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 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, 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. | |
| * | |
| * Here's an example demonstrating basic usage: | |
| * \include JacobiSVD_basic.cpp | |
| * Output: \verbinclude JacobiSVD_basic.out | |
| * | |
| * This JacobiSVD class is a two-sided Jacobi R-SVD decomposition, ensuring optimal reliability and accuracy. The | |
| * downside is that it's slower than 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 \a p is the greater dimension, meaning that it is still | |
| * of the same order of complexity as the faster bidiagonalizing R-SVD algorithms. In particular, like any R-SVD, it | |
| * takes advantage of non-squareness in that its complexity is only linear in the greater dimension. | |
| * | |
| * If the input matrix has inf or nan coefficients, the result of the computation is undefined, but the computation is | |
| * guaranteed to terminate in finite (and reasonable) time. | |
| * | |
| * The possible QR preconditioners that can be set with Options template parameter are: | |
| * \li ColPivHouseholderQRPreconditioner is the default. In practice it's very safe. It uses column-pivoting QR. | |
| * \li FullPivHouseholderQRPreconditioner, is the safest and slowest. It uses full-pivoting QR. | |
| * Contrary to other QRs, it doesn't allow computing thin unitaries. | |
| * \li HouseholderQRPreconditioner is the fastest, and less safe and accurate than the pivoting variants. It uses | |
| * non-pivoting QR. This is very similar in safety and accuracy to the bidiagonalization process used by bidiagonalizing | |
| * SVD algorithms (since bidiagonalization is inherently non-pivoting). However the resulting SVD is still more reliable | |
| * than bidiagonalizing SVDs because the Jacobi-based iterative process is more reliable than the optimized bidiagonal | |
| * SVD iterations. \li NoQRPreconditioner allows not to use a QR preconditioner at all. This is useful if you know that | |
| * you will only be computing JacobiSVD decompositions of square matrices. Non-square matrices require a QR | |
| * preconditioner. Using this option will result in faster compilation and smaller executable code. It won't | |
| * significantly speed up computation, since JacobiSVD is always checking if QR preconditioning is needed before | |
| * applying it anyway. | |
| * | |
| * One may also use the Options template parameter to specify how the unitaries should be computed. The options are | |
| * #ComputeThinU, #ComputeThinV, #ComputeFullU, #ComputeFullV. It is not possible to request both the thin and full | |
| * versions of a unitary. By default, unitaries will not be computed. | |
| * | |
| * You can set the QRPreconditioner and unitary options together: JacobiSVD<MatrixType, | |
| * ColPivHouseholderQRPreconditioner | ComputeThinU | ComputeFullV> | |
| * | |
| * \sa MatrixBase::jacobiSvd() | |
| */ | |
| template <typename MatrixType_, int Options_> | |
| class JacobiSVD : public SVDBase<JacobiSVD<MatrixType_, Options_>> { | |
| using Base = SVDBase<JacobiSVD>; | |
| public: | |
| using MatrixType = MatrixType_; | |
| using Scalar = typename Base::Scalar; | |
| using RealScalar = typename Base::RealScalar; | |
| enum : int { | |
| Options = Options_, | |
| QRPreconditioner = internal::get_qr_preconditioner(Options), | |
| RowsAtCompileTime = Base::RowsAtCompileTime, | |
| ColsAtCompileTime = Base::ColsAtCompileTime, | |
| DiagSizeAtCompileTime = Base::DiagSizeAtCompileTime, | |
| MaxRowsAtCompileTime = Base::MaxRowsAtCompileTime, | |
| MaxColsAtCompileTime = Base::MaxColsAtCompileTime, | |
| MaxDiagSizeAtCompileTime = Base::MaxDiagSizeAtCompileTime, | |
| MatrixOptions = Base::MatrixOptions | |
| }; | |
| using MatrixUType = typename Base::MatrixUType; | |
| using MatrixVType = typename Base::MatrixVType; | |
| using SingularValuesType = typename Base::SingularValuesType; | |
| using WorkMatrixType = Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime, MatrixOptions, | |
| MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime>; | |
| /** \brief Default Constructor. | |
| * | |
| * The default constructor is useful in cases in which the user intends to | |
| * perform decompositions via JacobiSVD::compute(const MatrixType&). | |
| */ | |
| JacobiSVD() {} | |
| /** \brief Default Constructor with memory preallocation | |
| * | |
| * Like the default constructor but with preallocation of the internal data | |
| * according to the specified problem size and \a Options template parameter. | |
| * | |
| * \sa JacobiSVD() | |
| */ | |
| JacobiSVD(Index rows, Index cols) { allocate(rows, cols, internal::get_computation_options(Options)); } | |
| /** \brief Default Constructor with memory preallocation | |
| * | |
| * Like the default constructor but with preallocation of the internal data | |
| * according to the specified problem size. | |
| * | |
| * One \b cannot request unitaries using both the \a Options template parameter | |
| * and the constructor. If possible, prefer using the \a Options template parameter. | |
| * | |
| * \param rows number of rows for the input matrix | |
| * \param cols number of columns for the input matrix | |
| * \param computationOptions specify whether to compute Thin/Full unitaries U/V | |
| * \sa JacobiSVD() | |
| * | |
| * \deprecated Will be removed in the next major Eigen version. Options should | |
| * be specified in the \a Options template parameter. | |
| */ | |
| EIGEN_DEPRECATED_WITH_REASON("Options should be specified using the class template parameter.") | |
| JacobiSVD(Index rows, Index cols, unsigned int computationOptions) { | |
| internal::check_svd_options_assertions<MatrixType, Options>(computationOptions, rows, cols); | |
| allocate(rows, cols, computationOptions); | |
| } | |
| /** \brief Constructor performing the decomposition of given matrix, using the custom options specified | |
| * with the \a Options template parameter. | |
| * | |
| * \param matrix the matrix to decompose | |
| */ | |
| template <typename Derived> | |
| explicit JacobiSVD(const MatrixBase<Derived>& matrix) { | |
| compute_impl(matrix, internal::get_computation_options(Options)); | |
| } | |
| template <typename Derived> | |
| explicit JacobiSVD(const TriangularBase<Derived>& matrix) { | |
| compute_impl(matrix, internal::get_computation_options(Options)); | |
| } | |
| /** \brief Constructor performing the decomposition of given matrix using specified options | |
| * for computing unitaries. | |
| * | |
| * One \b cannot request unitiaries using both the \a Options template parameter | |
| * and the constructor. If possible, prefer using the \a Options template parameter. | |
| * | |
| * \param matrix the matrix to decompose | |
| * \param computationOptions specify whether to compute Thin/Full unitaries U/V | |
| * | |
| * \deprecated Will be removed in the next major Eigen version. Options should | |
| * be specified in the \a Options template parameter. | |
| */ | |
| // EIGEN_DEPRECATED // TODO(cantonios): re-enable after fixing a few 3p libraries that error on deprecation warnings. | |
| template <typename Derived> | |
| JacobiSVD(const MatrixBase<Derived>& matrix, unsigned int computationOptions) { | |
| internal::check_svd_options_assertions<MatrixType, Options>(computationOptions, matrix.rows(), matrix.cols()); | |
| compute_impl(matrix, computationOptions); | |
| } | |
| /** \brief Method performing the decomposition of given matrix. Computes Thin/Full unitaries U/V if specified | |
| * using the \a Options template parameter or the class constructor. | |
| * | |
| * \param matrix the matrix to decompose | |
| */ | |
| template <typename Derived> | |
| JacobiSVD& compute(const MatrixBase<Derived>& matrix) { | |
| return compute_impl(matrix, m_computationOptions); | |
| } | |
| template <typename Derived> | |
| JacobiSVD& compute(const TriangularBase<Derived>& matrix) { | |
| return compute_impl(matrix, m_computationOptions); | |
| } | |
| /** \brief Method performing the decomposition of given matrix, as specified by | |
| * the `computationOptions` parameter. | |
| * | |
| * \param matrix the matrix to decompose | |
| * \param computationOptions specify whether to compute Thin/Full unitaries U/V | |
| * | |
| * \deprecated Will be removed in the next major Eigen version. Options should | |
| * be specified in the \a Options template parameter. | |
| */ | |
| template <typename Derived> | |
| EIGEN_DEPRECATED_WITH_REASON("Options should be specified using the class template parameter.") | |
| JacobiSVD& compute(const MatrixBase<Derived>& matrix, unsigned int computationOptions) { | |
| internal::check_svd_options_assertions<MatrixType, Options>(m_computationOptions, matrix.rows(), matrix.cols()); | |
| return compute_impl(matrix, computationOptions); | |
| } | |
| using Base::cols; | |
| using Base::computeU; | |
| using Base::computeV; | |
| using Base::diagSize; | |
| using Base::rank; | |
| using Base::rows; | |
| void allocate(Index rows_, Index cols_, unsigned int computationOptions) { | |
| if (Base::allocate(rows_, cols_, computationOptions)) return; | |
| eigen_assert(!(ShouldComputeThinU && int(QRPreconditioner) == int(FullPivHouseholderQRPreconditioner)) && | |
| !(ShouldComputeThinU && int(QRPreconditioner) == int(FullPivHouseholderQRPreconditioner)) && | |
| "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. " | |
| "Use the ColPivHouseholderQR preconditioner instead."); | |
| m_workMatrix.resize(diagSize(), diagSize()); | |
| if (cols() > rows()) m_qr_precond_morecols.allocate(*this); | |
| if (rows() > cols()) m_qr_precond_morerows.allocate(*this); | |
| } | |
| private: | |
| template <typename Derived> | |
| JacobiSVD& compute_impl(const TriangularBase<Derived>& matrix, unsigned int computationOptions); | |
| template <typename Derived> | |
| JacobiSVD& compute_impl(const MatrixBase<Derived>& matrix, unsigned int computationOptions); | |
| // Blocked sweep for the Jacobi SVD (works for both real and complex scalars). | |
| // Extracted into a separate EIGEN_DONT_INLINE method to prevent the blocking | |
| // code from interfering with the compiler's optimization of the non-blocking | |
| // scalar sweep. | |
| EIGEN_DONT_INLINE bool blocked_sweep(RealScalar considerAsZero, RealScalar precision, RealScalar& maxDiagEntry); | |
| protected: | |
| using Base::m_computationOptions; | |
| using Base::m_computeFullU; | |
| using Base::m_computeFullV; | |
| using Base::m_computeThinU; | |
| using Base::m_computeThinV; | |
| using Base::m_info; | |
| using Base::m_isAllocated; | |
| using Base::m_isInitialized; | |
| using Base::m_matrixU; | |
| using Base::m_matrixV; | |
| using Base::m_nonzeroSingularValues; | |
| using Base::m_prescribedThreshold; | |
| using Base::m_singularValues; | |
| using Base::m_usePrescribedThreshold; | |
| using Base::ShouldComputeThinU; | |
| using Base::ShouldComputeThinV; | |
| EIGEN_STATIC_ASSERT(!(ShouldComputeThinU && int(QRPreconditioner) == int(FullPivHouseholderQRPreconditioner)) && | |
| !(ShouldComputeThinU && int(QRPreconditioner) == int(FullPivHouseholderQRPreconditioner)), | |
| "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. " | |
| "Use the ColPivHouseholderQR preconditioner instead.") | |
| template <typename MatrixType__, int Options__, bool IsComplex_> | |
| friend struct internal::svd_precondition_2x2_block_to_be_real; | |
| template <typename MatrixType__, int Options__, int QRPreconditioner_, int Case_, bool DoAnything_> | |
| friend struct internal::qr_preconditioner_impl; | |
| internal::qr_preconditioner_impl<MatrixType, Options, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> | |
| m_qr_precond_morecols; | |
| internal::qr_preconditioner_impl<MatrixType, Options, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> | |
| m_qr_precond_morerows; | |
| WorkMatrixType m_workMatrix; | |
| // Blocking parameters for the Jacobi SVD sweep. | |
| #ifdef EIGEN_JACOBI_SVD_BLOCK_SIZE | |
| static constexpr Index kDefaultBlockSize = EIGEN_JACOBI_SVD_BLOCK_SIZE; | |
| #else | |
| static constexpr Index kDefaultBlockSize = 32; | |
| #endif | |
| // Use the lower of the default block size and static maximum matrix dimensions. | |
| static constexpr Index kBlockSize = internal::min_size_prefer_fixed(kDefaultBlockSize, MaxDiagSizeAtCompileTime); | |
| }; | |
| template <typename MatrixType, int Options> | |
| template <typename Derived> | |
| JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute_impl(const TriangularBase<Derived>& matrix, | |
| unsigned int computationOptions) { | |
| return compute_impl(matrix.toDenseMatrix(), computationOptions); | |
| } | |
| template <typename MatrixType, int Options> | |
| template <typename Derived> | |
| JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute_impl(const MatrixBase<Derived>& matrix, | |
| unsigned int computationOptions) { | |
| EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Derived, MatrixType); | |
| EIGEN_STATIC_ASSERT((std::is_same<typename Derived::Scalar, typename MatrixType::Scalar>::value), | |
| Input matrix must have the same Scalar type as the JacobiSVD object.); | |
| using numext::abs; | |
| allocate(matrix.rows(), matrix.cols(), computationOptions); | |
| // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number | |
| // of iterations, only worsening the precision of U and V as we accumulate more rotations | |
| const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon(); | |
| // limit for denormal numbers to be considered zero in order to avoid infinite loops (see bug 286) | |
| const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)(); | |
| // Scaling factor to reduce over/under-flows | |
| const RealScalar maxCoeff = matrix.cwiseAbs().template maxCoeff<PropagateNaN>(); | |
| if (!(numext::isfinite)(maxCoeff)) { | |
| m_isInitialized = true; | |
| m_info = InvalidInput; | |
| m_nonzeroSingularValues = 0; | |
| m_singularValues.setZero(); | |
| return *this; | |
| } | |
| internal::safe_scaling_factors<RealScalar> factors; | |
| /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */ | |
| if (rows() != cols()) { | |
| factors = | |
| internal::safe_scaling<RealScalar>::with_scaled(matrix.derived(), maxCoeff, [&](const auto& scaledMatrix) { | |
| m_qr_precond_morecols.run(*this, scaledMatrix); | |
| m_qr_precond_morerows.run(*this, scaledMatrix); | |
| }); | |
| } else { | |
| factors = internal::safe_scaling<RealScalar>::scale_to( | |
| m_workMatrix, | |
| matrix.template topLeftCorner<DiagSizeAtCompileTime, DiagSizeAtCompileTime>(diagSize(), diagSize()), maxCoeff); | |
| if (m_computeFullU) m_matrixU.setIdentity(rows(), rows()); | |
| if (m_computeThinU) m_matrixU.setIdentity(rows(), diagSize()); | |
| if (m_computeFullV) m_matrixV.setIdentity(cols(), cols()); | |
| if (m_computeThinV) m_matrixV.setIdentity(cols(), diagSize()); | |
| } | |
| /*** step 2. The main Jacobi SVD iteration. ***/ | |
| RealScalar maxDiagEntry = m_workMatrix.cwiseAbs().diagonal().maxCoeff(); | |
| bool finished = false; | |
| while (!finished) { | |
| finished = true; | |
| { | |
| // Sweep with optional blocking for large matrices. | |
| // Use blocking when the matrix is large enough that individual left rotations | |
| // (strided row operations on column-major data) cause significant cache misses. | |
| // The threshold is derived from the L2 cache size: blocking becomes worthwhile | |
| // when n exceeds sqrt(L2 / 4). We divide by sizeof(float) rather than sizeof(RealScalar) | |
| // because the cache miss pattern depends on the number of columns accessed (one cache | |
| // line per column), not the scalar size. This also makes the threshold appropriately | |
| // more conservative for larger types where GEMM overhead is higher. | |
| const Index n = diagSize(); | |
| #ifdef EIGEN_JACOBI_SVD_BLOCKING_THRESHOLD | |
| const Index blockingThreshold = EIGEN_JACOBI_SVD_BLOCKING_THRESHOLD; | |
| #else | |
| const Index blockingThreshold = | |
| static_cast<Index>(numext::sqrt(static_cast<double>(l2CacheSize() / sizeof(float)))); | |
| #endif | |
| if (n >= blockingThreshold) { | |
| // The blocked sweep is in a separate EIGEN_DONT_INLINE method to prevent | |
| // the blocking code from interfering with the compiler's optimization of | |
| // the non-blocking scalar sweep below. | |
| finished = !blocked_sweep(considerAsZero, precision, maxDiagEntry); | |
| } else | |
| finished = !internal::jacobi_svd_nonblocking_sweep(m_workMatrix, m_matrixU, m_matrixV, computeU(), computeV(), | |
| considerAsZero, precision, maxDiagEntry); | |
| } | |
| } | |
| /*** step 3. The work matrix is now diagonal, so ensure it's positive so its diagonal entries are the singular values | |
| * ***/ | |
| for (Index i = 0; i < diagSize(); ++i) { | |
| // For a complex matrix, some diagonal coefficients might not have been | |
| // treated by svd_precondition_2x2_block_to_be_real, and the imaginary part | |
| // of some diagonal entry might not be null. | |
| bool diagonal_has_imaginary_part = false; | |
| EIGEN_IF_CONSTEXPR (NumTraits<Scalar>::IsComplex) { | |
| diagonal_has_imaginary_part = abs(numext::imag(m_workMatrix.coeff(i, i))) > considerAsZero; | |
| } | |
| if (diagonal_has_imaginary_part) { | |
| RealScalar a = abs(m_workMatrix.coeff(i, i)); | |
| m_singularValues.coeffRef(i) = abs(a); | |
| if (computeU()) m_matrixU.col(i) *= m_workMatrix.coeff(i, i) / a; | |
| } else { | |
| // m_workMatrix.coeff(i,i) is already real, no difficulty: | |
| RealScalar a = numext::real(m_workMatrix.coeff(i, i)); | |
| m_singularValues.coeffRef(i) = abs(a); | |
| if (computeU() && (a < RealScalar(0))) m_matrixU.col(i) = -m_matrixU.col(i); | |
| } | |
| } | |
| internal::safe_scaling<RealScalar>::unscale_in_place(m_singularValues, factors); | |
| /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/ | |
| m_nonzeroSingularValues = diagSize(); | |
| for (Index i = 0; i < diagSize(); i++) { | |
| Index pos; | |
| RealScalar maxRemainingSingularValue = m_singularValues.tail(diagSize() - i).maxCoeff(&pos); | |
| if (numext::is_exactly_zero(maxRemainingSingularValue)) { | |
| m_nonzeroSingularValues = i; | |
| break; | |
| } | |
| if (pos) { | |
| pos += i; | |
| std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos)); | |
| if (computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i)); | |
| if (computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i)); | |
| } | |
| } | |
| m_isInitialized = true; | |
| return *this; | |
| } | |
| // Blocked Jacobi SVD sweep for both real and complex scalar types. For large n, | |
| // applying left rotations (row operations on column-major data) causes cache | |
| // misses due to strided access. To mitigate this, we accumulate kBlockSize left | |
| // rotations into a small dense matrix and apply them via a single GEMM to the | |
| // contiguous row block q..q+kBlockSize-1 and the (possibly distant) row p. | |
| // Right rotations and column scalings act on columns (contiguous in column-major) | |
| // and are applied individually. | |
| // | |
| // For complex types, the 2x2 preconditioning (making the block real) involves | |
| // complex left rotations and row scalings, which are also accumulated into the | |
| // block matrix. Column scalings from preconditioning are applied directly. | |
| // | |
| // The accumulated rotation matrix has lower-triangular structure in its top-left | |
| // kBlockSize x kBlockSize corner, which we exploit with triangularView. | |
| // | |
| // Returns true if any off-diagonal element exceeded the threshold (i.e. sweep | |
| // is not yet converged). | |
| template <typename MatrixType, int Options> | |
| EIGEN_DONT_INLINE bool JacobiSVD<MatrixType, Options>::blocked_sweep(RealScalar considerAsZero, RealScalar precision, | |
| RealScalar& maxDiagEntry) { | |
| using numext::abs; | |
| using numext::sqrt; | |
| const Index n = diagSize(); | |
| RealScalar threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry); | |
| bool notFinished = false; | |
| static constexpr Index kBlockBufferSize = (kBlockSize + 1) * (kBlockSize + 1); | |
| ei_declare_aligned_stack_constructed_variable(Scalar, blockBufferPtr, kBlockBufferSize, 0); | |
| Map<Matrix<Scalar, kBlockSize + 1, kBlockSize + 1, MatrixOptions>, AlignedMax> blockBuffer( | |
| blockBufferPtr, kBlockSize + 1, kBlockSize + 1); | |
| ei_declare_aligned_stack_constructed_variable(Scalar, accumPtr, kBlockBufferSize, 0); | |
| Map<Matrix<Scalar, kBlockSize + 1, kBlockSize + 1, MatrixOptions>, AlignedMax> accum(accumPtr, kBlockSize + 1, | |
| kBlockSize + 1); | |
| for (Index p = 1; p < n; ++p) { | |
| Index q = 0; | |
| // Blocked loop: process kBlockSize pairs (p,q+qq) for qq=0..kBlockSize-1. | |
| // We extract the relevant (kBlockSize+1) x (kBlockSize+1) submatrix of W | |
| // into a small buffer, compute all rotations on the buffer, accumulate the | |
| // left transformations into `accum`, and apply them in one GEMM at the end. | |
| for (; q + kBlockSize <= p; q += kBlockSize) { | |
| // Buffer = [ W(q:q+k, q:q+k) W(q:q+k, p) ] | |
| // [ W(p, q:q+k) W(p, p) ] | |
| blockBuffer.template topLeftCorner<kBlockSize, kBlockSize>() = | |
| m_workMatrix.template block<kBlockSize, kBlockSize>(q, q); | |
| blockBuffer.col(kBlockSize).template head<kBlockSize>() = m_workMatrix.col(p).template segment<kBlockSize>(q); | |
| blockBuffer.row(kBlockSize).template head<kBlockSize>() = m_workMatrix.row(p).template segment<kBlockSize>(q); | |
| blockBuffer(kBlockSize, kBlockSize) = m_workMatrix(p, p); | |
| // Accumulator for left transformations: W <- accum * W. | |
| // After processing qq pairs, accum's top-left kBlockSize x kBlockSize | |
| // block is lower-triangular (each rotation only mixes row qq with row | |
| // kBlockSize, so rows 0..qq-1 are unchanged). | |
| accum.setIdentity(kBlockSize + 1, kBlockSize + 1); | |
| bool blockDirty = false; | |
| for (Index qq = 0; qq < kBlockSize; ++qq) { | |
| if (abs(blockBuffer.coeff(kBlockSize, qq)) > threshold || abs(blockBuffer.coeff(qq, kBlockSize)) > threshold) { | |
| notFinished = true; | |
| blockDirty = true; | |
| // Complex preconditioning: transform the 2x2 block | |
| // [w_pp w_pq] = [buffer(kBlockSize, kBlockSize) buffer(kBlockSize, qq)] | |
| // [w_qp w_qq] [buffer(qq, kBlockSize) buffer(qq, qq) ] | |
| // to have real entries via unitary row/column operations, so | |
| // real_2x2_jacobi_svd can be applied. | |
| // | |
| // Left operations (complex rotation, row scaling by e^{i*theta}) are | |
| // accumulated into `accum` for deferred GEMM application. | |
| // Right operations (column scaling) are applied directly since column | |
| // ops are contiguous in column-major layout. | |
| bool doRealSvd = true; | |
| EIGEN_IF_CONSTEXPR (NumTraits<Scalar>::IsComplex) { | |
| Scalar z; | |
| // nn = ||(w_pp, w_qp)||_2, the norm of the first column of the 2x2 block. | |
| RealScalar nn = sqrt(numext::abs2(blockBuffer.coeff(kBlockSize, kBlockSize)) + | |
| numext::abs2(blockBuffer.coeff(qq, kBlockSize))); | |
| if (numext::is_exactly_zero(nn)) { | |
| // First column is zero => block is already upper triangular. | |
| blockBuffer.coeffRef(kBlockSize, kBlockSize) = Scalar(0); | |
| blockBuffer.coeffRef(qq, kBlockSize) = Scalar(0); | |
| // Scale rows by z = e^{-i*arg(w)} to make remaining entries real. | |
| if (abs(numext::imag(blockBuffer.coeff(kBlockSize, qq))) > considerAsZero) { | |
| z = abs(blockBuffer.coeff(kBlockSize, qq)) / blockBuffer.coeff(kBlockSize, qq); | |
| blockBuffer.row(kBlockSize) *= z; | |
| accum.row(kBlockSize) *= z; // accumulate left op | |
| if (computeU()) m_matrixU.col(p) *= numext::conj(z); | |
| } | |
| if (abs(numext::imag(blockBuffer.coeff(qq, qq))) > considerAsZero) { | |
| z = abs(blockBuffer.coeff(qq, qq)) / blockBuffer.coeff(qq, qq); | |
| blockBuffer.row(qq) *= z; | |
| accum.row(qq) *= z; // accumulate left op | |
| if (computeU()) m_matrixU.col(q + qq) *= numext::conj(z); | |
| } | |
| } else { | |
| // Apply complex Givens rotation to zero out w_qp: | |
| // [c s] [w_pp] [nn] conj(w_pp) w_qp | |
| // [-s c] [w_qp] = [0 ] c = ----------, s = ------ | |
| // nn nn | |
| JacobiRotation<Scalar> rot; | |
| rot.c() = numext::conj(blockBuffer.coeff(kBlockSize, kBlockSize)) / nn; | |
| rot.s() = blockBuffer.coeff(qq, kBlockSize) / nn; | |
| blockBuffer.applyOnTheLeft(kBlockSize, qq, rot); | |
| accum.applyOnTheLeft(kBlockSize, qq, rot); // accumulate left op | |
| if (computeU()) m_matrixU.applyOnTheRight(p, q + qq, rot.adjoint()); | |
| // Scale column qq by z = e^{-i*arg(w_pq)} to make w_pq real. | |
| if (abs(numext::imag(blockBuffer.coeff(kBlockSize, qq))) > considerAsZero) { | |
| z = abs(blockBuffer.coeff(kBlockSize, qq)) / blockBuffer.coeff(kBlockSize, qq); | |
| blockBuffer.col(qq) *= z; | |
| m_workMatrix.col(q + qq) *= z; // right op: apply directly | |
| if (computeV()) m_matrixV.col(q + qq) *= z; | |
| } | |
| // Scale row qq by z = e^{-i*arg(w_qq)} to make w_qq real. | |
| if (abs(numext::imag(blockBuffer.coeff(qq, qq))) > considerAsZero) { | |
| z = abs(blockBuffer.coeff(qq, qq)) / blockBuffer.coeff(qq, qq); | |
| blockBuffer.row(qq) *= z; | |
| accum.row(qq) *= z; // accumulate left op | |
| if (computeU()) m_matrixU.col(q + qq) *= numext::conj(z); | |
| } | |
| } | |
| // Update maxDiagEntry from preconditioning. | |
| maxDiagEntry = numext::maxi<RealScalar>( | |
| maxDiagEntry, numext::maxi<RealScalar>(abs(blockBuffer.coeff(kBlockSize, kBlockSize)), | |
| abs(blockBuffer.coeff(qq, qq)))); | |
| // Check if 2x2 block still needs diagonalizing. | |
| RealScalar precondThreshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry); | |
| doRealSvd = abs(blockBuffer.coeff(kBlockSize, qq)) > precondThreshold || | |
| abs(blockBuffer.coeff(qq, kBlockSize)) > precondThreshold; | |
| } | |
| if (doRealSvd) { | |
| // Compute real 2x2 SVD: buffer_2x2 = j_left * diag * j_right^T. | |
| JacobiRotation<RealScalar> j_left, j_right; | |
| internal::real_2x2_jacobi_svd(blockBuffer, kBlockSize, qq, &j_left, &j_right); | |
| blockBuffer.applyOnTheLeft(kBlockSize, qq, j_left); | |
| blockBuffer.applyOnTheRight(kBlockSize, qq, j_right); | |
| // Accumulate left rotation for deferred GEMM. | |
| accum.applyOnTheLeft(kBlockSize, qq, j_left); | |
| // Right rotation is a column op (contiguous): apply directly. | |
| m_workMatrix.applyOnTheRight(p, q + qq, j_right); | |
| if (computeU()) m_matrixU.applyOnTheRight(p, q + qq, j_left.transpose()); | |
| if (computeV()) m_matrixV.applyOnTheRight(p, q + qq, j_right); | |
| maxDiagEntry = numext::maxi<RealScalar>( | |
| maxDiagEntry, numext::maxi<RealScalar>(abs(blockBuffer.coeff(kBlockSize, kBlockSize)), | |
| abs(blockBuffer.coeff(qq, qq)))); | |
| } | |
| threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry); | |
| } | |
| } | |
| // Apply accumulated left rotations: W <- accum * W, via GEMM. | |
| // When p == q + kBlockSize, all kBlockSize+1 rows are contiguous. | |
| // Otherwise, rows q..q+k-1 and row p are non-adjacent; we split: | |
| // [Mq] [L11 l12] [Mq] | |
| // [Mp] <- [l21 l22] [Mp] | |
| // L11 is lower-triangular (exploited via triangularView). | |
| if (blockDirty) { | |
| if (p == q + kBlockSize) { | |
| m_workMatrix.template middleRows<kBlockSize + 1>(q) = | |
| accum * m_workMatrix.template middleRows<kBlockSize + 1>(q); | |
| } else { | |
| const auto L11 = accum.template topLeftCorner<kBlockSize, kBlockSize>(); | |
| const auto l12 = accum.col(kBlockSize).template head<kBlockSize>(); | |
| const auto l21 = accum.row(kBlockSize).template head<kBlockSize>(); | |
| const Scalar l22 = accum(kBlockSize, kBlockSize); | |
| auto Mq = m_workMatrix.template middleRows<kBlockSize>(q); | |
| auto Mp = m_workMatrix.row(p); | |
| Matrix<Scalar, 1, Dynamic, RowMajor, 1, MaxDiagSizeAtCompileTime> Mp_save = Mp; | |
| Mp.noalias() = l21 * Mq + l22 * Mp_save; | |
| Mq = L11.template triangularView<Lower>() * Mq + l12 * Mp_save; | |
| } | |
| } | |
| } | |
| // Scalar loop for remaining pairs after blocked processing. | |
| for (; q < p; ++q) { | |
| if (abs(m_workMatrix.coeff(p, q)) > threshold || abs(m_workMatrix.coeff(q, p)) > threshold) { | |
| notFinished = true; | |
| bool doRealSvd = true; | |
| EIGEN_IF_CONSTEXPR (NumTraits<Scalar>::IsComplex) { | |
| doRealSvd = internal::svd_precondition_2x2_block_to_be_real<MatrixType, Options>::run(m_workMatrix, *this, p, | |
| q, maxDiagEntry); | |
| } | |
| if (doRealSvd) { | |
| JacobiRotation<RealScalar> j_left, j_right; | |
| internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right); | |
| m_workMatrix.applyOnTheLeft(p, q, j_left); | |
| if (computeU()) m_matrixU.applyOnTheRight(p, q, j_left.transpose()); | |
| m_workMatrix.applyOnTheRight(p, q, j_right); | |
| if (computeV()) m_matrixV.applyOnTheRight(p, q, j_right); | |
| maxDiagEntry = numext::maxi<RealScalar>( | |
| maxDiagEntry, numext::maxi<RealScalar>(abs(m_workMatrix.coeff(p, p)), abs(m_workMatrix.coeff(q, q)))); | |
| } | |
| threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry); | |
| } | |
| } | |
| } | |
| return notFinished; | |
| } | |
| /** \svd_module | |
| * | |
| * \return the singular value decomposition of \c *this computed by two-sided | |
| * Jacobi transformations. | |
| * | |
| * \sa class JacobiSVD | |
| */ | |
| template <typename Derived> | |
| template <int Options> | |
| JacobiSVD<typename MatrixBase<Derived>::PlainObject, Options> MatrixBase<Derived>::jacobiSvd() const { | |
| return JacobiSVD<PlainObject, Options>(*this); | |
| } | |
| template <typename Derived> | |
| template <int Options> | |
| JacobiSVD<typename MatrixBase<Derived>::PlainObject, Options> MatrixBase<Derived>::jacobiSvd( | |
| unsigned int computationOptions) const { | |
| return JacobiSVD<PlainObject, Options>(*this, computationOptions); | |
| } | |
| } // end namespace Eigen | |
| #endif // EIGEN_JACOBISVD_H |