mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-18 14:34:17 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix bug #836: extend SparseQR to support more columns than rows.

Browse Source
This commit is contained in:
Gael Guennebaud 2014-07-01 10:24:46 +02:00
parent d73ee84d37
commit 75e574275c
2 changed files with 84 additions and 65 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

114
Eigen/src/SparseQR/SparseQR.h
View File

@ -2,7 +2,7 @@
// for linear algebra.
//
// Copyright (C) 2012-2013 Desire Nuentsa <desire.nuentsa_wakam@inria.fr>
// Copyright (C) 2012-2013 Gael Guennebaud <gael.guennebaud@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
@ -180,7 +180,7 @@ class SparseQR
y.bottomRows(y.rows()-rank).setZero();
// Apply the column permutation
if (m_perm_c.size()) dest.topRows(cols()) = colsPermutation() * y.topRows(cols());
if (m_perm_c.size()) dest = colsPermutation() * y.topRows(cols());
else dest = y.topRows(cols());
m_info = Success;
@ -286,6 +286,7 @@ void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat)
ord(mat, m_perm_c);
Index n = mat.cols();
Index m = mat.rows();
Index diagSize = (std::min)(m,n);
if (!m_perm_c.size())
{
@ -297,13 +298,13 @@ void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat)
m_outputPerm_c = m_perm_c.inverse();
internal::coletree(mat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
m_R.resize(n, n);
m_Q.resize(m, n);
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(n);
m_hcoeffs.resize(diagSize);
m_analysisIsok = true;
}
@ -323,11 +324,12 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step");
Index m = mat.rows();
Index n = mat.cols();
IndexVector mark(m); 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
bool found_diag;
Index 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 = mat;
m_pmat.uncompress(); // To have the innerNonZeroPtr allocated
@ -339,7 +341,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
m_pmat.innerNonZeroPtr()[p] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i];
}
/* Compute the default threshold, see :
/* 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
*/
@ -347,24 +349,24 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
{
RealScalar max2Norm = 0.0;
for (int j = 0; j < n; j++) max2Norm = (max)(max2Norm, m_pmat.col(j).norm());
m_threshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon();
pivotThreshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon();
}
// Initialize the numerical permutation
m_pivotperm.setIdentity(n);
Index 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 (Index col = 0; col < (std::min)(n,m); ++col)
for (Index col = 0; col < n; ++col)
{
mark.setConstant(-1);
m_R.startVec(col);
m_Q.startVec(col);
mark(nonzeroCol) = col;
Qidx(0) = nonzeroCol;
nzcolR = 0; nzcolQ = 1;
found_diag = col>=m;
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.,
@ -373,7 +375,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
// 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 MatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
{
Index curIdx = nonzeroCol ;
Index curIdx = nonzeroCol;
if(itp) curIdx = itp.row();
if(curIdx == nonzeroCol) found_diag = true;
@ -415,7 +417,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
// Browse all the indexes of R(:,col) in reverse order
for (Index i = nzcolR-1; i >= 0; i--)
{
Index curIdx = m_pivotperm.indices()(Ridx(i));
Index curIdx = Ridx(i);
// Apply the curIdx-th householder vector to the current column (temporarily stored into tval)
Scalar tdot(0);
@ -444,34 +446,37 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
}
}
} // End update current column
// Compute the Householder reflection that eliminate the current column
// FIXME this step should call the Householder module.
Scalar tau;
RealScalar beta;
Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);
RealScalar beta = 0;
// First, the squared norm of Q((col+1):m, col)
RealScalar sqrNorm = 0.;
for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq)));
if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0))
if(nonzeroCol < diagSize)
{
tau = RealScalar(0);
beta = numext::real(c0);
tval(Qidx(0)) = 1;
}
else
{
using std::sqrt;
beta = sqrt(numext::abs2(c0) + sqrNorm);
if(numext::real(c0) >= RealScalar(0))
beta = -beta;
tval(Qidx(0)) = 1;
for (Index itq = 1; itq < nzcolQ; ++itq)
tval(Qidx(itq)) /= (c0 - beta);
tau = numext::conj((beta-c0) / beta);
// 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);
// First, the squared norm of Q((col+1):m, col)
RealScalar sqrNorm = 0.;
for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq)));
if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0))
{
tau = RealScalar(0);
beta = numext::real(c0);
tval(Qidx(0)) = 1;
}
else
{
using std::sqrt;
beta = sqrt(numext::abs2(c0) + sqrNorm);
if(numext::real(c0) >= RealScalar(0))
beta = -beta;
tval(Qidx(0)) = 1;
for (Index itq = 1; itq < nzcolQ; ++itq)
tval(Qidx(itq)) /= (c0 - beta);
tau = numext::conj((beta-c0) / beta);
}
}
// Insert values in R
@ -485,24 +490,25 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
}
}
if(abs(beta) >= m_threshold)
if(nonzeroCol < diagSize && abs(beta) >= pivotThreshold)
{
m_R.insertBackByOuterInner(col, nonzeroCol) = beta;
nonzeroCol++;
// The householder coefficient
m_hcoeffs(col) = tau;
m_hcoeffs(nonzeroCol) = tau;
// Record the householder reflections
for (Index itq = 0; itq < nzcolQ; ++itq)
{
Index iQ = Qidx(itq);
m_Q.insertBackByOuterInnerUnordered(col,iQ) = tval(iQ);
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
m_hcoeffs(col) = Scalar(0);
for (Index j = nonzeroCol; j < n-1; j++)
std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]);
@ -511,6 +517,8 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
}
}
m_hcoeffs.tail(diagSize-nonzeroCol).setZero();
// Finalize the column pointers of the sparse matrices R and Q
m_Q.finalize();
m_Q.makeCompressed();
@ -579,14 +587,16 @@ struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived
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 < n; k++)
for (Index k = 0; k < diagSize; k++)
{
Scalar tau = Scalar(0);
tau = m_qr.m_Q.col(k).dot(res.col(j));
@ -599,10 +609,10 @@ struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived
else
{
eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes");
// Compute res = Q' * other column by column
// Compute res = Q * other column by column
for(Index j = 0; j < res.cols(); j++)
{
for (Index k = n-1; k >=0; k--)
for (Index k = diagSize-1; k >=0; k--)
{
Scalar tau = Scalar(0);
tau = m_qr.m_Q.col(k).dot(res.col(j));
@ -636,7 +646,7 @@ struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<Sp
return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
}
inline Index rows() const { return m_qr.rows(); }
inline Index cols() const { return m_qr.cols(); }
inline Index cols() const { return (std::min)(m_qr.rows(),m_qr.cols()); }
// To use for operations with the transpose of Q
SparseQRMatrixQTransposeReturnType<SparseQRType> transpose() const
{

35
test/sparseqr.cpp
View File

@ -2,24 +2,24 @@
// for linear algebra.
//
// Copyright (C) 2012 Desire Nuentsa Wakam <desire.nuentsa_wakam@inria.fr>
// Copyright (C) 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
#include "sparse.h"
#include <Eigen/SparseQR>
template<typename MatrixType,typename DenseMat>
int generate_sparse_rectangular_problem(MatrixType& A, DenseMat& dA, int maxRows = 300, int maxCols = 300)
int generate_sparse_rectangular_problem(MatrixType& A, DenseMat& dA, int maxRows = 300, int maxCols = 150)
{
eigen_assert(maxRows >= maxCols);
typedef typename MatrixType::Scalar Scalar;
int rows = internal::random<int>(1,maxRows);
int cols = internal::random<int>(1,rows);
int cols = internal::random<int>(1,maxCols);
double density = (std::max)(8./(rows*cols), 0.01);
A.resize(rows,rows);
dA.resize(rows,rows);
A.resize(rows,cols);
dA.resize(rows,cols);
initSparse<Scalar>(density, dA, A,ForceNonZeroDiag);
A.makeCompressed();
int nop = internal::random<int>(0, internal::random<double>(0,1) > 0.5 ? cols/2 : 0);
@ -31,6 +31,13 @@ int generate_sparse_rectangular_problem(MatrixType& A, DenseMat& dA, int maxRows
A.col(j0) = s * A.col(j1);
dA.col(j0) = s * dA.col(j1);
}
// if(rows<cols) {
// A.conservativeResize(cols,cols);
// dA.conservativeResize(cols,cols);
// dA.bottomRows(cols-rows).setZero();
// }
return rows;
}
@ -42,11 +49,10 @@ template<typename Scalar> void test_sparseqr_scalar()
MatrixType A;
DenseMat dA;
DenseVector refX,x,b;
SparseQR<MatrixType, AMDOrdering<int> > solver;
SparseQR<MatrixType, COLAMDOrdering<int> > solver;
generate_sparse_rectangular_problem(A,dA);
int n = A.cols();
b = DenseVector::Random(n);
b = dA * DenseVector::Random(A.cols());
solver.compute(A);
if (solver.info() != Success)
{
@ -60,17 +66,19 @@ template<typename Scalar> void test_sparseqr_scalar()
std::cerr << "sparse QR factorization failed\n";
exit(0);
return;
}
}
VERIFY_IS_APPROX(A * x, b);
//Compare with a dense QR solver
ColPivHouseholderQR<DenseMat> dqr(dA);
refX = dqr.solve(b);
VERIFY_IS_EQUAL(dqr.rank(), solver.rank());
if(solver.rank()<A.cols())
VERIFY((dA * refX - b).norm() * 2 > (A * x - b).norm() );
else
if(solver.rank()==A.cols()) // full rank
VERIFY_IS_APPROX(x, refX);
// else
// VERIFY((dA * refX - b).norm() * 2 > (A * x - b).norm() );
// Compute explicitly the matrix Q
MatrixType Q, QtQ, idM;
@ -88,3 +96,4 @@ void test_sparseqr()
CALL_SUBTEST_2(test_sparseqr_scalar<std::complex<double> >());
}
}