// This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2012-2013 Desire Nuentsa <desire.nuentsa_wakam@inria.fr> // Copyright (C) 2012-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/.
/** * \ingroup SparseQR_Module * \class SparseQR * \brief Sparse left-looking QR factorization with numerical column pivoting * * This class implements a left-looking QR decomposition of sparse matrices * with numerical column pivoting. * When a column has a norm less than a given tolerance * it is implicitly permuted to the end. The QR factorization thus obtained is * given by A*P = Q*R where R is upper triangular or trapezoidal. * * P is the column permutation which is the product of the fill-reducing and the * numerical permutations. Use colsPermutation() to get it. * * Q is the orthogonal matrix represented as products of Householder reflectors. * Use matrixQ() to get an expression and matrixQ().adjoint() to get the adjoint. * You can then apply it to a vector. * * R is the sparse triangular or trapezoidal matrix. The later occurs when A is rank-deficient. * matrixR().topLeftCorner(rank(), rank()) always returns a triangular factor of full rank. * * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<> * \tparam _OrderingType The fill-reducing ordering method. See the \link OrderingMethods_Module * OrderingMethods \endlink module for the list of built-in and external ordering methods. * * \implsparsesolverconcept * * The numerical pivoting strategy and default threshold are the same as in SuiteSparse QR, and * detailed in the following paper: * <i> * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011. * </i> * Even though it is qualified as "rank-revealing", this strategy might fail for some * rank deficient problems. When this class is used to solve linear or least-square problems * it is thus strongly recommended to check the accuracy of the computed solution. If it * failed, it usually helps to increase the threshold with setPivotThreshold. * * \warning The input sparse matrix A must be in compressed mode (see SparseMatrix::makeCompressed()). * \warning For complex matrices matrixQ().transpose() will actually return the adjoint matrix. *
*/ template<typename _MatrixType, typename _OrderingType> class SparseQR : public SparseSolverBase<SparseQR<_MatrixType,_OrderingType> >
{ protected: typedef SparseSolverBase<SparseQR<_MatrixType,_OrderingType> > Base; using Base::m_isInitialized; public: using Base::_solve_impl; typedef _MatrixType MatrixType; typedef _OrderingType OrderingType; typedeftypename MatrixType::Scalar Scalar; typedeftypename MatrixType::RealScalar RealScalar; typedeftypename MatrixType::StorageIndex StorageIndex; typedef SparseMatrix<Scalar,ColMajor,StorageIndex> QRMatrixType; typedef Matrix<StorageIndex, Dynamic, 1> IndexVector; typedef Matrix<Scalar, Dynamic, 1> ScalarVector; typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;
/** Construct a QR factorization of the matrix \a mat. * * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()). * * \sa compute()
*/ explicit SparseQR(const MatrixType& mat) : m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
{
compute(mat);
}
/** Computes the QR factorization of the sparse matrix \a mat. * * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()). * * \sa analyzePattern(), factorize()
*/ void compute(const MatrixType& mat)
{
analyzePattern(mat);
factorize(mat);
} void analyzePattern(const MatrixType& mat); void factorize(const MatrixType& mat);
/** \returns the number of rows of the represented matrix.
*/ inline Index rows() const { return m_pmat.rows(); }
/** \returns the number of columns of the represented matrix.
*/ inline Index cols() const { return m_pmat.cols();}
/** \returns a const reference to the \b sparse upper triangular matrix R of the QR factorization. * \warning The entries of the returned matrix are not sorted. This means that using it in algorithms * expecting sorted entries will fail. This include random coefficient accesses (SpaseMatrix::coeff()), * and coefficient-wise operations. Matrix products and triangular solves are fine though. * * To sort the entries, you can assign it to a row-major matrix, and if a column-major matrix * is required, you can copy it again: * \code * SparseMatrix<double> R = qr.matrixR(); // column-major, not sorted! * SparseMatrix<double,RowMajor> Rr = qr.matrixR(); // row-major, sorted * SparseMatrix<double> Rc = Rr; // column-major, sorted * \endcode
*/ const QRMatrixType& matrixR() const { return m_R; }
/** \returns the number of non linearly dependent columns as determined by the pivoting threshold. * * \sa setPivotThreshold()
*/
Index rank() const
{
eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); return m_nonzeropivots;
}
/** \returns an expression of the matrix Q as products of sparse Householder reflectors. * The common usage of this function is to apply it to a dense matrix or vector * \code * VectorXd B1, B2; * // Initialize B1 * B2 = matrixQ() * B1; * \endcode * * To get a plain SparseMatrix representation of Q: * \code * SparseMatrix<double> Q; * Q = SparseQR<SparseMatrix<double> >(A).matrixQ(); * \endcode * Internally, this call simply performs a sparse product between the matrix Q * and a sparse identity matrix. However, due to the fact that the sparse * reflectors are stored unsorted, two transpositions are needed to sort * them before performing the product.
*/
SparseQRMatrixQReturnType<SparseQR> matrixQ() const
{ return SparseQRMatrixQReturnType<SparseQR>(*this); }
/** \returns a const reference to the column permutation P that was applied to A such that A*P = Q*R * It is the combination of the fill-in reducing permutation and numerical column pivoting.
*/ const PermutationType& colsPermutation() const
{
eigen_assert(m_isInitialized && "Decomposition is not initialized."); return m_outputPerm_c;
}
/** \returns A string describing the type of error. * This method is provided to ease debugging, not to handle errors.
*/
std::string lastErrorMessage() const { return m_lastError; }
/** \internal */ template<typename Rhs, typename Dest> bool _solve_impl(const MatrixBase<Rhs> &B, MatrixBase<Dest> &dest) const
{
eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
Index rank = this->rank();
// Compute Q^* * b; typename Dest::PlainObject y, b;
y = this->matrixQ().adjoint() * B;
b = y;
// Solve with the triangular matrix R
y.resize((std::max<Index>)(cols(),y.rows()),y.cols());
y.topRows(rank) = this->matrixR().topLeftCorner(rank, rank).template triangularView<Upper>().solve(b.topRows(rank));
y.bottomRows(y.rows()-rank).setZero();
// Apply the column permutation if (m_perm_c.size()) dest = colsPermutation() * y.topRows(cols()); else dest = y.topRows(cols());
m_info = Success; returntrue;
}
/** Sets the threshold that is used to determine linearly dependent columns during the factorization. * * In practice, if during the factorization the norm of the column that has to be eliminated is below * this threshold, then the entire column is treated as zero, and it is moved at the end.
*/ void setPivotThreshold(const RealScalar& threshold)
{
m_useDefaultThreshold = false;
m_threshold = threshold;
}
/** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. * * \sa compute()
*/ template<typename Rhs> inlineconst Solve<SparseQR, Rhs> solve(const MatrixBase<Rhs>& B) const
{
eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); return Solve<SparseQR, Rhs>(*this, B.derived());
} template<typename Rhs> inlineconst Solve<SparseQR, Rhs> solve(const SparseMatrixBase<Rhs>& B) const
{
eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); return Solve<SparseQR, Rhs>(*this, B.derived());
}
/** \brief Reports whether previous computation was successful. * * \returns \c Success if computation was successful, * \c NumericalIssue if the QR factorization reports a numerical problem * \c InvalidInput if the input matrix is invalid * * \sa iparm()
*/
ComputationInfo info() const
{
eigen_assert(m_isInitialized && "Decomposition is not initialized."); return m_info;
}
/** \internal */ inlinevoid _sort_matrix_Q()
{ if(this->m_isQSorted) return; // The matrix Q is sorted during the transposition
SparseMatrix<Scalar, RowMajor, Index> mQrm(this->m_Q);
this->m_Q = mQrm;
this->m_isQSorted = true;
}
protected: bool m_analysisIsok; bool m_factorizationIsok; mutable ComputationInfo m_info;
std::string m_lastError;
QRMatrixType m_pmat; // Temporary matrix
QRMatrixType m_R; // The triangular factor matrix
QRMatrixType m_Q; // The orthogonal reflectors
ScalarVector m_hcoeffs; // The Householder coefficients
PermutationType m_perm_c; // Fill-reducing Column permutation
PermutationType m_pivotperm; // The permutation for rank revealing
PermutationType m_outputPerm_c; // The final column permutation
RealScalar m_threshold; // Threshold to determine null Householder reflections bool m_useDefaultThreshold; // Use default threshold
Index m_nonzeropivots; // Number of non zero pivots found
IndexVector m_etree; // Column elimination tree
IndexVector m_firstRowElt; // First element in each row bool m_isQSorted; // whether Q is sorted or not bool m_isEtreeOk; // whether the elimination tree match the initial input matrix
/** \brief Preprocessing step of a QR factorization * * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()). * * In this step, the fill-reducing permutation is computed and applied to the columns of A * and the column elimination tree is computed as well. Only the sparsity pattern of \a mat is exploited. * * \note In this step it is assumed that there is no empty row in the matrix \a mat.
*/ template <typename MatrixType, typename OrderingType> void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat)
{
eigen_assert(mat.isCompressed() && "SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR"); // Copy to a column major matrix if the input is rowmajor typename internal::conditional<MatrixType::IsRowMajor,QRMatrixType,const MatrixType&>::type matCpy(mat); // Compute the column fill reducing ordering
OrderingType ord;
ord(matCpy, m_perm_c);
Index n = mat.cols();
Index m = mat.rows();
Index diagSize = (std::min)(m,n);
if (!m_perm_c.size())
{
m_perm_c.resize(n);
m_perm_c.indices().setLinSpaced(n, 0,StorageIndex(n-1));
}
// Compute the column elimination tree of the permuted matrix
m_outputPerm_c = m_perm_c.inverse();
internal::coletree(matCpy, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
m_isEtreeOk = true;
m_R.resize(m, n);
m_Q.resize(m, diagSize);
// Allocate space for nonzero elements: rough estimation
m_R.reserve(2*mat.nonZeros()); //FIXME Get a more accurate estimation through symbolic factorization with the etree
m_Q.reserve(2*mat.nonZeros());
m_hcoeffs.resize(diagSize);
m_analysisIsok = true;
}
/** \brief Performs the numerical QR factorization of the input matrix * * The function SparseQR::analyzePattern(const MatrixType&) must have been called beforehand with * a matrix having the same sparsity pattern than \a mat. * * \param mat The sparse column-major matrix
*/ template <typename MatrixType, typename OrderingType> void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
{ using std::abs;
eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step");
StorageIndex m = StorageIndex(mat.rows());
StorageIndex n = StorageIndex(mat.cols());
StorageIndex diagSize = (std::min)(m,n);
IndexVector mark((std::max)(m,n)); mark.setConstant(-1); // Record the visited nodes
IndexVector Ridx(n), Qidx(m); // Store temporarily the row indexes for the current column of R and Q
Index nzcolR, nzcolQ; // Number of nonzero for the current column of R and Q
ScalarVector tval(m); // The dense vector used to compute the current column
RealScalar pivotThreshold = m_threshold;
m_pmat.uncompress(); // To have the innerNonZeroPtr allocated
// Apply the fill-in reducing permutation lazily:
{ // If the input is row major, copy the original column indices, // otherwise directly use the input matrix //
IndexVector originalOuterIndicesCpy; const StorageIndex *originalOuterIndices = mat.outerIndexPtr(); if(MatrixType::IsRowMajor)
{
originalOuterIndicesCpy = IndexVector::Map(m_pmat.outerIndexPtr(),n+1);
originalOuterIndices = originalOuterIndicesCpy.data();
}
for (int i = 0; i < n; i++)
{
Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i;
m_pmat.outerIndexPtr()[p] = originalOuterIndices[i];
m_pmat.innerNonZeroPtr()[p] = originalOuterIndices[i+1] - originalOuterIndices[i];
}
}
/* Compute the default threshold as in MatLab, see: * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3
*/ if(m_useDefaultThreshold)
{
RealScalar max2Norm = 0.0; for (int j = 0; j < n; j++) max2Norm = numext::maxi(max2Norm, m_pmat.col(j).norm()); if(max2Norm==RealScalar(0))
max2Norm = RealScalar(1);
pivotThreshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon();
}
// Initialize the numerical permutation
m_pivotperm.setIdentity(n);
StorageIndex nonzeroCol = 0; // Record the number of valid pivots
m_Q.startVec(0);
// Left looking rank-revealing QR factorization: compute a column of R and Q at a time for (StorageIndex col = 0; col < n; ++col)
{
mark.setConstant(-1);
m_R.startVec(col);
mark(nonzeroCol) = col;
Qidx(0) = nonzeroCol;
nzcolR = 0; nzcolQ = 1; bool found_diag = nonzeroCol>=m;
tval.setZero();
// Symbolic factorization: find the nonzero locations of the column k of the factors R and Q, i.e., // all the nodes (with indexes lower than rank) reachable through the column elimination tree (etree) rooted at node k. // Note: if the diagonal entry does not exist, then its contribution must be explicitly added, // thus the trick with found_diag that permits to do one more iteration on the diagonal element if this one has not been found. for (typename QRMatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
{
StorageIndex curIdx = nonzeroCol; if(itp) curIdx = StorageIndex(itp.row()); if(curIdx == nonzeroCol) found_diag = true;
// Get the nonzeros indexes of the current column of R
StorageIndex st = m_firstRowElt(curIdx); // The traversal of the etree starts here if (st < 0 )
{
m_lastError = "Empty row found during numerical factorization";
m_info = InvalidInput; return;
}
// Traverse the etree
Index bi = nzcolR; for (; mark(st) != col; st = m_etree(st))
{
Ridx(nzcolR) = st; // Add this row to the list,
mark(st) = col; // and mark this row as visited
nzcolR++;
}
// Reverse the list to get the topological ordering
Index nt = nzcolR-bi; for(Index i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1));
// Copy the current (curIdx,pcol) value of the input matrix if(itp) tval(curIdx) = itp.value(); else tval(curIdx) = Scalar(0);
// Compute the pattern of Q(:,k) if(curIdx > nonzeroCol && mark(curIdx) != col )
{
Qidx(nzcolQ) = curIdx; // Add this row to the pattern of Q,
mark(curIdx) = col; // and mark it as visited
nzcolQ++;
}
}
// Browse all the indexes of R(:,col) in reverse order for (Index i = nzcolR-1; i >= 0; i--)
{
Index curIdx = Ridx(i);
// Apply the curIdx-th householder vector to the current column (temporarily stored into tval)
Scalar tdot(0);
// First compute q' * tval
tdot = m_Q.col(curIdx).dot(tval);
tdot *= m_hcoeffs(curIdx);
// Then update tval = tval - q * tau // FIXME: tval -= tdot * m_Q.col(curIdx) should amount to the same (need to check/add support for efficient "dense ?= sparse") for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
tval(itq.row()) -= itq.value() * tdot;
// Detect fill-in for the current column of Q if(m_etree(Ridx(i)) == nonzeroCol)
{ for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
{
StorageIndex iQ = StorageIndex(itq.row()); if (mark(iQ) != col)
{
Qidx(nzcolQ++) = iQ; // Add this row to the pattern of Q,
mark(iQ) = col; // and mark it as visited
}
}
}
} // End update current column
Scalar tau = RealScalar(0);
RealScalar beta = 0;
if(nonzeroCol < diagSize)
{ // Compute the Householder reflection that eliminate the current column // FIXME this step should call the Householder module.
Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);
// Insert values in R for (Index i = nzcolR-1; i >= 0; i--)
{
Index curIdx = Ridx(i); if(curIdx < nonzeroCol)
{
m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx);
tval(curIdx) = Scalar(0.);
}
}
if(nonzeroCol < diagSize && abs(beta) >= pivotThreshold)
{
m_R.insertBackByOuterInner(col, nonzeroCol) = beta; // The householder coefficient
m_hcoeffs(nonzeroCol) = tau; // Record the householder reflections for (Index itq = 0; itq < nzcolQ; ++itq)
{
Index iQ = Qidx(itq);
m_Q.insertBackByOuterInnerUnordered(nonzeroCol,iQ) = tval(iQ);
tval(iQ) = Scalar(0.);
}
nonzeroCol++; if(nonzeroCol<diagSize)
m_Q.startVec(nonzeroCol);
} else
{ // Zero pivot found: move implicitly this column to the end for (Index j = nonzeroCol; j < n-1; j++)
std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]);
// Recompute the column elimination tree
internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data());
m_isEtreeOk = false;
}
}
m_hcoeffs.tail(diagSize-nonzeroCol).setZero();
// Finalize the column pointers of the sparse matrices R and Q
m_Q.finalize();
m_Q.makeCompressed();
m_R.finalize();
m_R.makeCompressed();
m_isQSorted = false;
m_nonzeropivots = nonzeroCol;
if(nonzeroCol<n)
{ // Permute the triangular factor to put the 'dead' columns to the end
QRMatrixType tempR(m_R);
m_R = tempR * m_pivotperm;
// Assign to a vector template<typename DesType> void evalTo(DesType& res) const
{
Index m = m_qr.rows();
Index n = m_qr.cols();
Index diagSize = (std::min)(m,n);
res = m_other; if (m_transpose)
{
eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes"); //Compute res = Q' * other column by column for(Index j = 0; j < res.cols(); j++){ for (Index k = 0; k < diagSize; k++)
{
Scalar tau = Scalar(0);
tau = m_qr.m_Q.col(k).dot(res.col(j)); if(tau==Scalar(0)) continue;
tau = tau * m_qr.m_hcoeffs(k);
res.col(j) -= tau * m_qr.m_Q.col(k);
}
}
} else
{
eigen_assert(m_qr.matrixQ().cols() == m_other.rows() && "Non conforming object sizes");
res.conservativeResize(rows(), cols());
// Compute res = Q * other column by column for(Index j = 0; j < res.cols(); j++)
{
Index start_k = internal::is_identity<Derived>::value ? numext::mini(j,diagSize-1) : diagSize-1; for (Index k = start_k; k >=0; k--)
{
Scalar tau = Scalar(0);
tau = m_qr.m_Q.col(k).dot(res.col(j)); if(tau==Scalar(0)) continue;
tau = tau * numext::conj(m_qr.m_hcoeffs(k));
res.col(j) -= tau * m_qr.m_Q.col(k);
}
}
}
}
const SparseQRType& m_qr; const Derived& m_other; bool m_transpose; // TODO this actually means adjoint
};
template<typename SparseQRType> struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<SparseQRType> >
{ typedeftypename SparseQRType::Scalar Scalar; typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; enum {
RowsAtCompileTime = Dynamic,
ColsAtCompileTime = Dynamic
}; explicit SparseQRMatrixQReturnType(const SparseQRType& qr) : m_qr(qr) {} template<typename Derived>
SparseQR_QProduct<SparseQRType, Derived> operator*(const MatrixBase<Derived>& other)
{ return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),false);
} // To use for operations with the adjoint of Q
SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint() const
{ return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
} inline Index rows() const { return m_qr.rows(); } inline Index cols() const { return m_qr.rows(); } // To use for operations with the transpose of Q FIXME this is the same as adjoint at the moment
SparseQRMatrixQTransposeReturnType<SparseQRType> transpose() const
{ return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
} const SparseQRType& m_qr;
};
// TODO this actually represents the adjoint of Q template<typename SparseQRType> struct SparseQRMatrixQTransposeReturnType
{ explicit SparseQRMatrixQTransposeReturnType(const SparseQRType& qr) : m_qr(qr) {} template<typename Derived>
SparseQR_QProduct<SparseQRType,Derived> operator*(const MatrixBase<Derived>& other)
{ return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(), true);
} const SparseQRType& m_qr;
};
Die Informationen auf dieser Webseite wurden
nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit,
noch Qualität der bereit gestellten Informationen zugesichert.
Bemerkung:
Die farbliche Syntaxdarstellung ist noch experimentell.
