Arnoldi.h
Go to the documentation of this file.
1 // Copyright (C) 2018-2019 Yixuan Qiu <yixuan.qiu@cos.name>
2 //
3 // This Source Code Form is subject to the terms of the Mozilla
4 // Public License v. 2.0. If a copy of the MPL was not distributed
5 // with this file, You can obtain one at https://mozilla.org/MPL/2.0/.
6 
7 #ifndef ARNOLDI_H
8 #define ARNOLDI_H
9 
10 #include <Eigen/Core>
11 #include <cmath> // std::sqrt
12 #include <stdexcept> // std::invalid_argument
13 #include <sstream> // std::stringstream
14 
15 #include "../MatOp/internal/ArnoldiOp.h"
16 #include "../Util/TypeTraits.h"
17 #include "../Util/SimpleRandom.h"
18 #include "UpperHessenbergQR.h"
19 #include "DoubleShiftQR.h"
20 
21 namespace Spectra {
22 
23 // Arnoldi factorization A * V = V * H + f * e'
24 // A: n x n
25 // V: n x k
26 // H: k x k
27 // f: n x 1
28 // e: [0, ..., 0, 1]
29 // V and H are allocated of dimension m, so the maximum value of k is m
30 template <typename Scalar, typename ArnoldiOpType>
31 class Arnoldi
32 {
33 private:
41 
42 protected:
43  // clang-format off
44  ArnoldiOpType m_op; // Operators for the Arnoldi factorization
45 
46  const Index m_n; // dimension of A
47  const Index m_m; // maximum dimension of subspace V
48  Index m_k; // current dimension of subspace V
49 
50  Matrix m_fac_V; // V matrix in the Arnoldi factorization
51  Matrix m_fac_H; // H matrix in the Arnoldi factorization
52  Vector m_fac_f; // residual in the Arnoldi factorization
53  Scalar m_beta; // ||f||, B-norm of f
54 
55  const Scalar m_near_0; // a very small value, but 1.0 / m_near_0 does not overflow
56  // ~= 1e-307 for the "double" type
57  const Scalar m_eps; // the machine precision, ~= 1e-16 for the "double" type
58  // clang-format on
59 
60  // Given orthonormal basis functions V, find a nonzero vector f such that V'Bf = 0
61  // Assume that f has been properly allocated
62  void expand_basis(MapConstMat& V, const Index seed, Vector& f, Scalar& fnorm)
63  {
64  using std::sqrt;
65 
66  const Scalar thresh = m_eps * sqrt(Scalar(m_n));
67  Vector Vf(V.cols());
68  for (Index iter = 0; iter < 5; iter++)
69  {
70  // Randomly generate a new vector and orthogonalize it against V
71  SimpleRandom<Scalar> rng(seed + 123 * iter);
72  f.noalias() = rng.random_vec(m_n);
73  // f <- f - V * V'Bf, so that f is orthogonal to V in B-norm
74  m_op.trans_product(V, f, Vf);
75  f.noalias() -= V * Vf;
76  // fnorm <- ||f||
77  fnorm = m_op.norm(f);
78 
79  // If fnorm is too close to zero, we try a new random vector,
80  // otherwise return the result
81  if (fnorm >= thresh)
82  return;
83  }
84  }
85 
86 public:
87  Arnoldi(const ArnoldiOpType& op, Index m) :
88  m_op(op), m_n(op.rows()), m_m(m), m_k(0),
89  m_near_0(TypeTraits<Scalar>::min() * Scalar(10)),
90  m_eps(Eigen::NumTraits<Scalar>::epsilon())
91  {}
92 
93  virtual ~Arnoldi() {}
94 
95  // Const-reference to internal structures
96  const Matrix& matrix_V() const { return m_fac_V; }
97  const Matrix& matrix_H() const { return m_fac_H; }
98  const Vector& vector_f() const { return m_fac_f; }
99  Scalar f_norm() const { return m_beta; }
100  Index subspace_dim() const { return m_k; }
101 
102  // Initialize with an operator and an initial vector
103  void init(MapConstVec& v0, Index& op_counter)
104  {
105  m_fac_V.resize(m_n, m_m);
106  m_fac_H.resize(m_m, m_m);
107  m_fac_f.resize(m_n);
108  m_fac_H.setZero();
109 
110  // Verify the initial vector
111  const Scalar v0norm = m_op.norm(v0);
112  if (v0norm < m_near_0)
113  throw std::invalid_argument("initial residual vector cannot be zero");
114 
115  // Points to the first column of V
116  MapVec v(m_fac_V.data(), m_n);
117 
118  // Normalize
119  v.noalias() = v0 / v0norm;
120 
121  // Compute H and f
122  Vector w(m_n);
123  m_op.perform_op(v.data(), w.data());
124  op_counter++;
125 
126  m_fac_H(0, 0) = m_op.inner_product(v, w);
127  m_fac_f.noalias() = w - v * m_fac_H(0, 0);
128 
129  // In some cases f is zero in exact arithmetics, but due to rounding errors
130  // it may contain tiny fluctuations. When this happens, we force f to be zero
131  if (m_fac_f.cwiseAbs().maxCoeff() < m_eps)
132  {
133  m_fac_f.setZero();
134  m_beta = Scalar(0);
135  }
136  else
137  {
138  m_beta = m_op.norm(m_fac_f);
139  }
140 
141  // Indicate that this is a step-1 factorization
142  m_k = 1;
143  }
144 
145  // Arnoldi factorization starting from step-k
146  virtual void factorize_from(Index from_k, Index to_m, Index& op_counter)
147  {
148  using std::sqrt;
149 
150  if (to_m <= from_k)
151  return;
152 
153  if (from_k > m_k)
154  {
155  std::stringstream msg;
156  msg << "Arnoldi: from_k (= " << from_k << ") is larger than the current subspace dimension (= " << m_k << ")";
157  throw std::invalid_argument(msg.str());
158  }
159 
160  const Scalar beta_thresh = m_eps * sqrt(Scalar(m_n));
161 
162  // Pre-allocate vectors
163  Vector Vf(to_m);
164  Vector w(m_n);
165 
166  // Keep the upperleft k x k submatrix of H and set other elements to 0
167  m_fac_H.rightCols(m_m - from_k).setZero();
168  m_fac_H.block(from_k, 0, m_m - from_k, from_k).setZero();
169 
170  for (Index i = from_k; i <= to_m - 1; i++)
171  {
172  bool restart = false;
173  // If beta = 0, then the next V is not full rank
174  // We need to generate a new residual vector that is orthogonal
175  // to the current V, which we call a restart
176  if (m_beta < m_near_0)
177  {
178  MapConstMat V(m_fac_V.data(), m_n, i); // The first i columns
179  expand_basis(V, 2 * i, m_fac_f, m_beta);
180  restart = true;
181  }
182 
183  // v <- f / ||f||
184  m_fac_V.col(i).noalias() = m_fac_f / m_beta; // The (i+1)-th column
185 
186  // Note that H[i+1, i] equals to the unrestarted beta
187  m_fac_H(i, i - 1) = restart ? Scalar(0) : m_beta;
188 
189  // w <- A * v, v = m_fac_V.col(i)
190  m_op.perform_op(&m_fac_V(0, i), w.data());
191  op_counter++;
192 
193  const Index i1 = i + 1;
194  // First i+1 columns of V
195  MapConstMat Vs(m_fac_V.data(), m_n, i1);
196  // h = m_fac_H(0:i, i)
197  MapVec h(&m_fac_H(0, i), i1);
198  // h <- V'Bw
199  m_op.trans_product(Vs, w, h);
200 
201  // f <- w - V * h
202  m_fac_f.noalias() = w - Vs * h;
203  m_beta = m_op.norm(m_fac_f);
204 
205  if (m_beta > Scalar(0.717) * m_op.norm(h))
206  continue;
207 
208  // f/||f|| is going to be the next column of V, so we need to test
209  // whether V'B(f/||f||) ~= 0
210  m_op.trans_product(Vs, m_fac_f, Vf.head(i1));
211  Scalar ortho_err = Vf.head(i1).cwiseAbs().maxCoeff();
212  // If not, iteratively correct the residual
213  int count = 0;
214  while (count < 5 && ortho_err > m_eps * m_beta)
215  {
216  // There is an edge case: when beta=||f|| is close to zero, f mostly consists
217  // of noises of rounding errors, so the test [ortho_err < eps * beta] is very
218  // likely to fail. In particular, if beta=0, then the test is ensured to fail.
219  // Hence when this happens, we force f to be zero, and then restart in the
220  // next iteration.
221  if (m_beta < beta_thresh)
222  {
223  m_fac_f.setZero();
224  m_beta = Scalar(0);
225  break;
226  }
227 
228  // f <- f - V * Vf
229  m_fac_f.noalias() -= Vs * Vf.head(i1);
230  // h <- h + Vf
231  h.noalias() += Vf.head(i1);
232  // beta <- ||f||
233  m_beta = m_op.norm(m_fac_f);
234 
235  m_op.trans_product(Vs, m_fac_f, Vf.head(i1));
236  ortho_err = Vf.head(i1).cwiseAbs().maxCoeff();
237  count++;
238  }
239  }
240 
241  // Indicate that this is a step-m factorization
242  m_k = to_m;
243  }
244 
245  // Apply H -> Q'HQ, where Q is from a double shift QR decomposition
246  void compress_H(const DoubleShiftQR<Scalar>& decomp)
247  {
248  decomp.matrix_QtHQ(m_fac_H);
249  m_k -= 2;
250  }
251 
252  // Apply H -> Q'HQ, where Q is from an upper Hessenberg QR decomposition
254  {
255  decomp.matrix_QtHQ(m_fac_H);
256  m_k--;
257  }
258 
259  // Apply V -> VQ and compute the new f.
260  // Should be called after compress_H(), since m_k is updated there.
261  // Only need to update the first k+1 columns of V
262  // The first (m - k + i) elements of the i-th column of Q are non-zero,
263  // and the rest are zero
264  void compress_V(const Matrix& Q)
265  {
266  Matrix Vs(m_n, m_k + 1);
267  for (Index i = 0; i < m_k; i++)
268  {
269  const Index nnz = m_m - m_k + i + 1;
270  MapConstVec q(&Q(0, i), nnz);
271  Vs.col(i).noalias() = m_fac_V.leftCols(nnz) * q;
272  }
273  Vs.col(m_k).noalias() = m_fac_V * Q.col(m_k);
274  m_fac_V.leftCols(m_k + 1).noalias() = Vs;
275 
276  Vector fk = m_fac_f * Q(m_m - 1, m_k - 1) + m_fac_V.col(m_k) * m_fac_H(m_k, m_k - 1);
277  m_fac_f.swap(fk);
278  m_beta = m_op.norm(m_fac_f);
279  }
280 };
281 
282 } // namespace Spectra
283 
284 #endif // ARNOLDI_H
Matrix3f m
const Scalar m_eps
Definition: Arnoldi.h:57
SCALAR Scalar
Definition: bench_gemm.cpp:46
EIGEN_DEVICE_FUNC Derived & setZero(Index size)
const Index m_n
Definition: Arnoldi.h:46
Quaternion Q
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar * data() const
void expand_basis(MapConstMat &V, const Index seed, Vector &f, Scalar &fnorm)
Definition: Arnoldi.h:62
#define min(a, b)
Definition: datatypes.h:19
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > Vector
Definition: Arnoldi.h:36
void compress_H(const UpperHessenbergQR< Scalar > &decomp)
Definition: Arnoldi.h:253
static std::mt19937 rng
ArnoldiOpType m_op
Definition: Arnoldi.h:44
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
const Index m_m
Definition: Arnoldi.h:47
iterator iter(handle obj)
Definition: pytypes.h:2273
const Scalar m_near_0
Definition: Arnoldi.h:55
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void swap(DenseBase< OtherDerived > &other)
const Matrix & matrix_V() const
Definition: Arnoldi.h:96
void init(MapConstVec &v0, Index &op_counter)
Definition: Arnoldi.h:103
void matrix_QtHQ(Matrix &dest) const
static double epsilon
Definition: testRot3.cpp:37
Eigen::Index Index
Definition: Arnoldi.h:34
Scalar m_beta
Definition: Arnoldi.h:53
virtual void matrix_QtHQ(Matrix &dest) const
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
virtual ~Arnoldi()
Definition: Arnoldi.h:93
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Index subspace_dim() const
Definition: Arnoldi.h:100
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > Matrix
Definition: Arnoldi.h:35
Array< int, Dynamic, 1 > v
Vector m_fac_f
Definition: Arnoldi.h:52
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
Eigen::Map< Matrix > MapMat
Definition: Arnoldi.h:37
Matrix m_fac_V
Definition: Arnoldi.h:50
EIGEN_DEVICE_FUNC const Scalar & q
const Matrix & matrix_H() const
Definition: Arnoldi.h:97
RowVector3d w
void compress_V(const Matrix &Q)
Definition: Arnoldi.h:264
Eigen::Map< const Matrix > MapConstMat
Definition: Arnoldi.h:39
const double h
static const double v0
void compress_H(const DoubleShiftQR< Scalar > &decomp)
Definition: Arnoldi.h:246
Eigen::Map< const Vector > MapConstVec
Definition: Arnoldi.h:40
Scalar f_norm() const
Definition: Arnoldi.h:99
const Vector & vector_f() const
Definition: Arnoldi.h:98
Matrix m_fac_H
Definition: Arnoldi.h:51
Eigen::Map< Vector > MapVec
Definition: Arnoldi.h:38
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
virtual void factorize_from(Index from_k, Index to_m, Index &op_counter)
Definition: Arnoldi.h:146
Arnoldi(const ArnoldiOpType &op, Index m)
Definition: Arnoldi.h:87


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:33:53