Lanczos.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 LANCZOS_H
8 #define LANCZOS_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 "Arnoldi.h"
16 
17 namespace Spectra {
18 
19 // Lanczos factorization A * V = V * H + f * e'
20 // A: n x n
21 // V: n x k
22 // H: k x k
23 // f: n x 1
24 // e: [0, ..., 0, 1]
25 // V and H are allocated of dimension m, so the maximum value of k is m
26 template <typename Scalar, typename ArnoldiOpType>
27 class Lanczos : public Arnoldi<Scalar, ArnoldiOpType>
28 {
29 private:
37 
48 
49 public:
50  Lanczos(const ArnoldiOpType& op, Index m) :
51  Arnoldi<Scalar, ArnoldiOpType>(op, m)
52  {}
53 
54  // Lanczos factorization starting from step-k
55  void factorize_from(Index from_k, Index to_m, Index& op_counter)
56  {
57  using std::sqrt;
58 
59  if (to_m <= from_k)
60  return;
61 
62  if (from_k > m_k)
63  {
64  std::stringstream msg;
65  msg << "Lanczos: from_k (= " << from_k << ") is larger than the current subspace dimension (= " << m_k << ")";
66  throw std::invalid_argument(msg.str());
67  }
68 
69  const Scalar beta_thresh = m_eps * sqrt(Scalar(m_n));
70 
71  // Pre-allocate vectors
72  Vector Vf(to_m);
73  Vector w(m_n);
74 
75  // Keep the upperleft k x k submatrix of H and set other elements to 0
76  m_fac_H.rightCols(m_m - from_k).setZero();
77  m_fac_H.block(from_k, 0, m_m - from_k, from_k).setZero();
78 
79  for (Index i = from_k; i <= to_m - 1; i++)
80  {
81  bool restart = false;
82  // If beta = 0, then the next V is not full rank
83  // We need to generate a new residual vector that is orthogonal
84  // to the current V, which we call a restart
85  if (m_beta < m_near_0)
86  {
87  MapConstMat V(m_fac_V.data(), m_n, i); // The first i columns
88  this->expand_basis(V, 2 * i, m_fac_f, m_beta);
89  restart = true;
90  }
91 
92  // v <- f / ||f||
93  MapVec v(&m_fac_V(0, i), m_n); // The (i+1)-th column
94  v.noalias() = m_fac_f / m_beta;
95 
96  // Note that H[i+1, i] equals to the unrestarted beta
97  m_fac_H(i, i - 1) = restart ? Scalar(0) : m_beta;
98 
99  // w <- A * v
100  m_op.perform_op(v.data(), w.data());
101  op_counter++;
102 
103  // H[i+1, i+1] = <v, w> = v'Bw
104  m_fac_H(i - 1, i) = m_fac_H(i, i - 1); // Due to symmetry
105  m_fac_H(i, i) = m_op.inner_product(v, w);
106 
107  // f <- w - V * V'Bw = w - H[i+1, i] * V{i} - H[i+1, i+1] * V{i+1}
108  // If restarting, we know that H[i+1, i] = 0
109  if (restart)
110  m_fac_f.noalias() = w - m_fac_H(i, i) * v;
111  else
112  m_fac_f.noalias() = w - m_fac_H(i, i - 1) * m_fac_V.col(i - 1) - m_fac_H(i, i) * v;
113 
114  m_beta = m_op.norm(m_fac_f);
115 
116  // f/||f|| is going to be the next column of V, so we need to test
117  // whether V'B(f/||f||) ~= 0
118  const Index i1 = i + 1;
119  MapMat Vs(m_fac_V.data(), m_n, i1); // The first (i+1) columns
120  m_op.trans_product(Vs, m_fac_f, Vf.head(i1));
121  Scalar ortho_err = Vf.head(i1).cwiseAbs().maxCoeff();
122  // If not, iteratively correct the residual
123  int count = 0;
124  while (count < 5 && ortho_err > m_eps * m_beta)
125  {
126  // There is an edge case: when beta=||f|| is close to zero, f mostly consists
127  // of noises of rounding errors, so the test [ortho_err < eps * beta] is very
128  // likely to fail. In particular, if beta=0, then the test is ensured to fail.
129  // Hence when this happens, we force f to be zero, and then restart in the
130  // next iteration.
131  if (m_beta < beta_thresh)
132  {
133  m_fac_f.setZero();
134  m_beta = Scalar(0);
135  break;
136  }
137 
138  // f <- f - V * Vf
139  m_fac_f.noalias() -= Vs * Vf.head(i1);
140  // h <- h + Vf
141  m_fac_H(i - 1, i) += Vf[i - 1];
142  m_fac_H(i, i - 1) = m_fac_H(i - 1, i);
143  m_fac_H(i, i) += Vf[i];
144  // beta <- ||f||
145  m_beta = m_op.norm(m_fac_f);
146 
147  m_op.trans_product(Vs, m_fac_f, Vf.head(i1));
148  ortho_err = Vf.head(i1).cwiseAbs().maxCoeff();
149  count++;
150  }
151  }
152 
153  // Indicate that this is a step-m factorization
154  m_k = to_m;
155  }
156 
157  // Apply H -> Q'HQ, where Q is from a tridiagonal QR decomposition
158  void compress_H(const TridiagQR<Scalar>& decomp)
159  {
160  decomp.matrix_QtHQ(m_fac_H);
161  m_k--;
162  }
163 };
164 
165 } // namespace Spectra
166 
167 #endif // LANCZOS_H
Matrix3f m
Eigen::Map< const Vector > MapConstVec
Definition: Lanczos.h:36
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
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
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
ArnoldiOpType m_op
Definition: Arnoldi.h:44
const Index m_m
Definition: Arnoldi.h:47
const Scalar m_near_0
Definition: Arnoldi.h:55
Eigen::Map< Vector > MapVec
Definition: Lanczos.h:34
Scalar m_beta
Definition: Arnoldi.h:53
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > Vector
Definition: Lanczos.h:32
void compress_H(const TridiagQR< Scalar > &decomp)
Definition: Lanczos.h:158
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Array< int, Dynamic, 1 > v
Vector m_fac_f
Definition: Arnoldi.h:52
Matrix m_fac_V
Definition: Arnoldi.h:50
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > Matrix
Definition: Lanczos.h:31
RowVector3d w
Eigen::Map< Matrix > MapMat
Definition: Lanczos.h:33
Eigen::Map< const Matrix > MapConstMat
Definition: Lanczos.h:35
void matrix_QtHQ(Matrix &dest) const
Lanczos(const ArnoldiOpType &op, Index m)
Definition: Lanczos.h:50
Matrix m_fac_H
Definition: Arnoldi.h:51
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
Eigen::Index Index
Definition: Lanczos.h:30
void factorize_from(Index from_k, Index to_m, Index &op_counter)
Definition: Lanczos.h:55


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:31