sparseqr.cpp
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2012 Desire Nuentsa Wakam <desire.nuentsa_wakam@inria.fr>
5 // Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 #include "sparse.h"
10 #include <Eigen/SparseQR>
11 
12 template<typename MatrixType,typename DenseMat>
13 int generate_sparse_rectangular_problem(MatrixType& A, DenseMat& dA, int maxRows = 300, int maxCols = 150)
14 {
15  eigen_assert(maxRows >= maxCols);
16  typedef typename MatrixType::Scalar Scalar;
17  int rows = internal::random<int>(1,maxRows);
18  int cols = internal::random<int>(1,maxCols);
19  double density = (std::max)(8./(rows*cols), 0.01);
20 
21  A.resize(rows,cols);
22  dA.resize(rows,cols);
23  initSparse<Scalar>(density, dA, A,ForceNonZeroDiag);
24  A.makeCompressed();
25  int nop = internal::random<int>(0, internal::random<double>(0,1) > 0.5 ? cols/2 : 0);
26  for(int k=0; k<nop; ++k)
27  {
28  int j0 = internal::random<int>(0,cols-1);
29  int j1 = internal::random<int>(0,cols-1);
30  Scalar s = internal::random<Scalar>();
31  A.col(j0) = s * A.col(j1);
32  dA.col(j0) = s * dA.col(j1);
33  }
34 
35 // if(rows<cols) {
36 // A.conservativeResize(cols,cols);
37 // dA.conservativeResize(cols,cols);
38 // dA.bottomRows(cols-rows).setZero();
39 // }
40 
41  return rows;
42 }
43 
44 template<typename Scalar> void test_sparseqr_scalar()
45 {
46  typedef typename NumTraits<Scalar>::Real RealScalar;
48  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMat;
50  MatrixType A;
51  DenseMat dA;
52  DenseVector refX,x,b;
55 
56  b = dA * DenseVector::Random(A.cols());
57  solver.compute(A);
58 
59  // Q should be MxM
60  VERIFY_IS_EQUAL(solver.matrixQ().rows(), A.rows());
61  VERIFY_IS_EQUAL(solver.matrixQ().cols(), A.rows());
62 
63  // R should be MxN
64  VERIFY_IS_EQUAL(solver.matrixR().rows(), A.rows());
65  VERIFY_IS_EQUAL(solver.matrixR().cols(), A.cols());
66 
67  // Q and R can be multiplied
68  DenseMat recoveredA = solver.matrixQ()
69  * DenseMat(solver.matrixR().template triangularView<Upper>())
70  * solver.colsPermutation().transpose();
71  VERIFY_IS_EQUAL(recoveredA.rows(), A.rows());
72  VERIFY_IS_EQUAL(recoveredA.cols(), A.cols());
73 
74  // and in the full rank case the original matrix is recovered
75  if (solver.rank() == A.cols())
76  {
77  VERIFY_IS_APPROX(A, recoveredA);
78  }
79 
80  if(internal::random<float>(0,1)>0.5f)
81  solver.factorize(A); // this checks that calling analyzePattern is not needed if the pattern do not change.
82  if (solver.info() != Success)
83  {
84  std::cerr << "sparse QR factorization failed\n";
85  exit(0);
86  return;
87  }
88  x = solver.solve(b);
89  if (solver.info() != Success)
90  {
91  std::cerr << "sparse QR factorization failed\n";
92  exit(0);
93  return;
94  }
95 
96  // Compare with a dense QR solver
98  refX = dqr.solve(b);
99 
100  bool rank_deficient = A.cols()>A.rows() || dqr.rank()<A.cols();
101  if(rank_deficient)
102  {
103  // rank deficient problem -> we might have to increase the threshold
104  // to get a correct solution.
105  RealScalar th = RealScalar(20)*dA.colwise().norm().maxCoeff()*(A.rows()+A.cols()) * NumTraits<RealScalar>::epsilon();
106  for(Index k=0; (k<16) && !test_isApprox(A*x,b); ++k)
107  {
108  th *= RealScalar(10);
109  solver.setPivotThreshold(th);
110  solver.compute(A);
111  x = solver.solve(b);
112  }
113  }
114 
115  VERIFY_IS_APPROX(A * x, b);
116 
117  // For rank deficient problem, the estimated rank might
118  // be slightly off, so let's only raise a warning in such cases.
119  if(rank_deficient) ++g_test_level;
120  VERIFY_IS_EQUAL(solver.rank(), dqr.rank());
121  if(rank_deficient) --g_test_level;
122 
123  if(solver.rank()==A.cols()) // full rank
124  VERIFY_IS_APPROX(x, refX);
125 // else
126 // VERIFY((dA * refX - b).norm() * 2 > (A * x - b).norm() );
127 
128  // Compute explicitly the matrix Q
129  MatrixType Q, QtQ, idM;
130  Q = solver.matrixQ();
131  //Check ||Q' * Q - I ||
132  QtQ = Q * Q.adjoint();
133  idM.resize(Q.rows(), Q.rows()); idM.setIdentity();
134  VERIFY(idM.isApprox(QtQ));
135 
136  // Q to dense
137  DenseMat dQ;
138  dQ = solver.matrixQ();
139  VERIFY_IS_APPROX(Q, dQ);
140 }
142 {
143  for(int i=0; i<g_repeat; ++i)
144  {
145  CALL_SUBTEST_1(test_sparseqr_scalar<double>());
146  CALL_SUBTEST_2(test_sparseqr_scalar<std::complex<double> >());
147  }
148 }
149 
Eigen::SparseMatrix< Scalar, ColMajor >
s
RealScalar s
Definition: level1_cplx_impl.h:126
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
Eigen::ColPivHouseholderQR::rank
Index rank() const
Definition: ColPivHouseholderQR.h:255
VERIFY_IS_EQUAL
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:386
b
Scalar * b
Definition: benchVecAdd.cpp:17
eigen_assert
#define eigen_assert(x)
Definition: Macros.h:1037
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
generate_sparse_rectangular_problem
int generate_sparse_rectangular_problem(MatrixType &A, DenseMat &dA, int maxRows=300, int maxCols=150)
Definition: sparseqr.cpp:13
Eigen::Success
@ Success
Definition: Constants.h:442
Q
Quaternion Q
Definition: testQuaternion.cpp:27
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
solver
BiCGSTAB< SparseMatrix< double > > solver
Definition: BiCGSTAB_simple.cpp:5
A
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:48
test_isApprox
bool test_isApprox(const AnnoyingScalar &a, const AnnoyingScalar &b)
Definition: AnnoyingScalar.h:159
CALL_SUBTEST_1
#define CALL_SUBTEST_1(FUNC)
Definition: split_test_helper.h:4
density
Definition: testGaussianConditional.cpp:127
Eigen::g_repeat
static int g_repeat
Definition: main.h:169
Eigen::ColPivHouseholderQR
Householder rank-revealing QR decomposition of a matrix with column-pivoting.
Definition: ForwardDeclarations.h:274
CALL_SUBTEST_2
#define CALL_SUBTEST_2(FUNC)
Definition: split_test_helper.h:10
DenseVector
Matrix< Scalar, Dynamic, 1 > DenseVector
Definition: BenchSparseUtil.h:24
tree::f
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
Definition: testExpression.cpp:218
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:15
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
j0
double j0(double x)
Definition: j0.c:185
Eigen::Quaternion
The quaternion class used to represent 3D orientations and rotations.
Definition: ForwardDeclarations.h:293
test_sparseqr_scalar
void test_sparseqr_scalar()
Definition: sparseqr.cpp:44
EIGEN_DECLARE_TEST
EIGEN_DECLARE_TEST(sparseqr)
Definition: sparseqr.cpp:141
ForceNonZeroDiag
@ ForceNonZeroDiag
Definition: sparse.h:37
Eigen::Matrix< Scalar, Dynamic, Dynamic >
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
Eigen::g_test_level
static int g_test_level
Definition: main.h:168
max
#define max(a, b)
Definition: datatypes.h:20
j1
double j1(double x)
Definition: j1.c:174
Eigen::SparseQR
Sparse left-looking QR factorization with numerical column pivoting.
Definition: SparseQR.h:16
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:232
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
sparse.h
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
VERIFY
#define VERIFY(a)
Definition: main.h:380
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74


gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:04:44