1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
|
// This file is part of Eigen, a lightweight C++ template library
// 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 = 150)
{
eigen_assert(maxRows >= maxCols);
typedef typename MatrixType::Scalar Scalar;
int rows = internal::random<int>(1,maxRows);
int cols = internal::random<int>(1,maxCols);
double density = (std::max)(8./(rows*cols), 0.01);
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);
for(int k=0; k<nop; ++k)
{
int j0 = internal::random<int>(0,cols-1);
int j1 = internal::random<int>(0,cols-1);
Scalar s = internal::random<Scalar>();
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;
}
template<typename Scalar> void test_sparseqr_scalar()
{
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef SparseMatrix<Scalar,ColMajor> MatrixType;
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMat;
typedef Matrix<Scalar,Dynamic,1> DenseVector;
MatrixType A;
DenseMat dA;
DenseVector refX,x,b;
SparseQR<MatrixType, COLAMDOrdering<int> > solver;
generate_sparse_rectangular_problem(A,dA);
b = dA * DenseVector::Random(A.cols());
solver.compute(A);
// Q should be MxM
VERIFY_IS_EQUAL(solver.matrixQ().rows(), A.rows());
VERIFY_IS_EQUAL(solver.matrixQ().cols(), A.rows());
// R should be MxN
VERIFY_IS_EQUAL(solver.matrixR().rows(), A.rows());
VERIFY_IS_EQUAL(solver.matrixR().cols(), A.cols());
// Q and R can be multiplied
DenseMat recoveredA = solver.matrixQ()
* DenseMat(solver.matrixR().template triangularView<Upper>())
* solver.colsPermutation().transpose();
VERIFY_IS_EQUAL(recoveredA.rows(), A.rows());
VERIFY_IS_EQUAL(recoveredA.cols(), A.cols());
// and in the full rank case the original matrix is recovered
if (solver.rank() == A.cols())
{
VERIFY_IS_APPROX(A, recoveredA);
}
if(internal::random<float>(0,1)>0.5f)
solver.factorize(A); // this checks that calling analyzePattern is not needed if the pattern do not change.
if (solver.info() != Success)
{
std::cerr << "sparse QR factorization failed\n";
exit(0);
return;
}
x = solver.solve(b);
if (solver.info() != Success)
{
std::cerr << "sparse QR factorization failed\n";
exit(0);
return;
}
// Compare with a dense QR solver
ColPivHouseholderQR<DenseMat> dqr(dA);
refX = dqr.solve(b);
bool rank_deficient = A.cols()>A.rows() || dqr.rank()<A.cols();
if(rank_deficient)
{
// rank deficient problem -> we might have to increase the threshold
// to get a correct solution.
RealScalar th = RealScalar(20)*dA.colwise().norm().maxCoeff()*(A.rows()+A.cols()) * NumTraits<RealScalar>::epsilon();
for(Index k=0; (k<16) && !test_isApprox(A*x,b); ++k)
{
th *= RealScalar(10);
solver.setPivotThreshold(th);
solver.compute(A);
x = solver.solve(b);
}
}
VERIFY_IS_APPROX(A * x, b);
// For rank deficient problem, the estimated rank might
// be slightly off, so let's only raise a warning in such cases.
if(rank_deficient) ++g_test_level;
VERIFY_IS_EQUAL(solver.rank(), dqr.rank());
if(rank_deficient) --g_test_level;
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;
Q = solver.matrixQ();
//Check ||Q' * Q - I ||
QtQ = Q * Q.adjoint();
idM.resize(Q.rows(), Q.rows()); idM.setIdentity();
VERIFY(idM.isApprox(QtQ));
// Q to dense
DenseMat dQ;
dQ = solver.matrixQ();
VERIFY_IS_APPROX(Q, dQ);
}
EIGEN_DECLARE_TEST(sparseqr)
{
for(int i=0; i<g_repeat; ++i)
{
CALL_SUBTEST_1(test_sparseqr_scalar<double>());
CALL_SUBTEST_2(test_sparseqr_scalar<std::complex<double> >());
}
}
|
