SymEigsBase.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 SYM_EIGS_BASE_H
8 #define SYM_EIGS_BASE_H
9 
10 #include <Eigen/Core>
11 #include <vector> // std::vector
12 #include <cmath> // std::abs, std::pow, std::sqrt
13 #include <algorithm> // std::min, std::copy
14 #include <stdexcept> // std::invalid_argument
15 
16 #include "Util/TypeTraits.h"
17 #include "Util/SelectionRule.h"
18 #include "Util/CompInfo.h"
19 #include "Util/SimpleRandom.h"
22 #include "LinAlg/TridiagEigen.h"
23 #include "LinAlg/Lanczos.h"
24 
25 namespace Spectra {
26 
32 
40 template <typename Scalar,
41  int SelectionRule,
42  typename OpType,
43  typename BOpType>
45 {
46 private:
55 
58 
59 protected:
60  // clang-format off
61  OpType* m_op; // object to conduct matrix operation,
62  // e.g. matrix-vector product
63  const Index m_n; // dimension of matrix A
64  const Index m_nev; // number of eigenvalues requested
65  const Index m_ncv; // dimension of Krylov subspace in the Lanczos method
66  Index m_nmatop; // number of matrix operations called
67  Index m_niter; // number of restarting iterations
68 
69  LanczosFac m_fac; // Lanczos factorization
70  Vector m_ritz_val; // Ritz values
71 
72 private:
73  Matrix m_ritz_vec; // Ritz vectors
74  Vector m_ritz_est; // last row of m_ritz_vec, also called the Ritz estimates
75  BoolArray m_ritz_conv; // indicator of the convergence of Ritz values
76  int m_info; // status of the computation
77 
78  const Scalar m_near_0; // a very small value, but 1.0 / m_near_0 does not overflow
79  // ~= 1e-307 for the "double" type
80  const Scalar m_eps; // the machine precision, ~= 1e-16 for the "double" type
81  const Scalar m_eps23; // m_eps^(2/3), used to test the convergence
82  // clang-format on
83 
84  // Implicitly restarted Lanczos factorization
85  void restart(Index k)
86  {
87  if (k >= m_ncv)
88  return;
89 
90  TridiagQR<Scalar> decomp(m_ncv);
91  Matrix Q = Matrix::Identity(m_ncv, m_ncv);
92 
93  for (Index i = k; i < m_ncv; i++)
94  {
95  // QR decomposition of H-mu*I, mu is the shift
96  decomp.compute(m_fac.matrix_H(), m_ritz_val[i]);
97 
98  // Q -> Q * Qi
99  decomp.apply_YQ(Q);
100  // H -> Q'HQ
101  // Since QR = H - mu * I, we have H = QR + mu * I
102  // and therefore Q'HQ = RQ + mu * I
103  m_fac.compress_H(decomp);
104  }
105 
106  m_fac.compress_V(Q);
107  m_fac.factorize_from(k, m_ncv, m_nmatop);
108 
110  }
111 
112  // Calculates the number of converged Ritz values
114  {
115  // thresh = tol * max(m_eps23, abs(theta)), theta for Ritz value
116  Array thresh = tol * m_ritz_val.head(m_nev).array().abs().max(m_eps23);
117  Array resid = m_ritz_est.head(m_nev).array().abs() * m_fac.f_norm();
118  // Converged "wanted" Ritz values
119  m_ritz_conv = (resid < thresh);
120 
121  return m_ritz_conv.cast<Index>().sum();
122  }
123 
124  // Returns the adjusted nev for restarting
125  Index nev_adjusted(Index nconv)
126  {
127  using std::abs;
128 
129  Index nev_new = m_nev;
130  for (Index i = m_nev; i < m_ncv; i++)
131  if (abs(m_ritz_est[i]) < m_near_0)
132  nev_new++;
133 
134  // Adjust nev_new, according to dsaup2.f line 677~684 in ARPACK
135  nev_new += std::min(nconv, (m_ncv - nev_new) / 2);
136  if (nev_new == 1 && m_ncv >= 6)
137  nev_new = m_ncv / 2;
138  else if (nev_new == 1 && m_ncv > 2)
139  nev_new = 2;
140 
141  if (nev_new > m_ncv - 1)
142  nev_new = m_ncv - 1;
143 
144  return nev_new;
145  }
146 
147  // Retrieves and sorts Ritz values and Ritz vectors
149  {
150  TridiagEigen<Scalar> decomp(m_fac.matrix_H());
151  const Vector& evals = decomp.eigenvalues();
152  const Matrix& evecs = decomp.eigenvectors();
153 
154  SortEigenvalue<Scalar, SelectionRule> sorting(evals.data(), evals.size());
155  std::vector<int> ind = sorting.index();
156 
157  // For BOTH_ENDS, the eigenvalues are sorted according
158  // to the LARGEST_ALGE rule, so we need to move those smallest
159  // values to the left
160  // The order would be
161  // Largest => Smallest => 2nd largest => 2nd smallest => ...
162  // We keep this order since the first k values will always be
163  // the wanted collection, no matter k is nev_updated (used in restart())
164  // or is nev (used in sort_ritzpair())
165  if (SelectionRule == BOTH_ENDS)
166  {
167  std::vector<int> ind_copy(ind);
168  for (Index i = 0; i < m_ncv; i++)
169  {
170  // If i is even, pick values from the left (large values)
171  // If i is odd, pick values from the right (small values)
172  if (i % 2 == 0)
173  ind[i] = ind_copy[i / 2];
174  else
175  ind[i] = ind_copy[m_ncv - 1 - i / 2];
176  }
177  }
178 
179  // Copy the Ritz values and vectors to m_ritz_val and m_ritz_vec, respectively
180  for (Index i = 0; i < m_ncv; i++)
181  {
182  m_ritz_val[i] = evals[ind[i]];
183  m_ritz_est[i] = evecs(m_ncv - 1, ind[i]);
184  }
185  for (Index i = 0; i < m_nev; i++)
186  {
187  m_ritz_vec.col(i).noalias() = evecs.col(ind[i]);
188  }
189  }
190 
191 protected:
192  // Sorts the first nev Ritz pairs in the specified order
193  // This is used to return the final results
194  virtual void sort_ritzpair(int sort_rule)
195  {
196  // First make sure that we have a valid index vector
197  SortEigenvalue<Scalar, LARGEST_ALGE> sorting(m_ritz_val.data(), m_nev);
198  std::vector<int> ind = sorting.index();
199 
200  switch (sort_rule)
201  {
202  case LARGEST_ALGE:
203  break;
204  case LARGEST_MAGN:
205  {
206  SortEigenvalue<Scalar, LARGEST_MAGN> sorting(m_ritz_val.data(), m_nev);
207  ind = sorting.index();
208  break;
209  }
210  case SMALLEST_ALGE:
211  {
212  SortEigenvalue<Scalar, SMALLEST_ALGE> sorting(m_ritz_val.data(), m_nev);
213  ind = sorting.index();
214  break;
215  }
216  case SMALLEST_MAGN:
217  {
218  SortEigenvalue<Scalar, SMALLEST_MAGN> sorting(m_ritz_val.data(), m_nev);
219  ind = sorting.index();
220  break;
221  }
222  default:
223  throw std::invalid_argument("unsupported sorting rule");
224  }
225 
226  Vector new_ritz_val(m_ncv);
227  Matrix new_ritz_vec(m_ncv, m_nev);
228  BoolArray new_ritz_conv(m_nev);
229 
230  for (Index i = 0; i < m_nev; i++)
231  {
232  new_ritz_val[i] = m_ritz_val[ind[i]];
233  new_ritz_vec.col(i).noalias() = m_ritz_vec.col(ind[i]);
234  new_ritz_conv[i] = m_ritz_conv[ind[i]];
235  }
236 
237  m_ritz_val.swap(new_ritz_val);
238  m_ritz_vec.swap(new_ritz_vec);
239  m_ritz_conv.swap(new_ritz_conv);
240  }
241 
242 public:
244 
245  SymEigsBase(OpType* op, BOpType* Bop, Index nev, Index ncv) :
246  m_op(op),
247  m_n(m_op->rows()),
248  m_nev(nev),
249  m_ncv(ncv > m_n ? m_n : ncv),
250  m_nmatop(0),
251  m_niter(0),
252  m_fac(ArnoldiOpType(op, Bop), m_ncv),
256  m_eps23(Eigen::numext::pow(m_eps, Scalar(2.0) / 3))
257  {
258  if (nev < 1 || nev > m_n - 1)
259  throw std::invalid_argument("nev must satisfy 1 <= nev <= n - 1, n is the size of matrix");
260 
261  if (ncv <= nev || ncv > m_n)
262  throw std::invalid_argument("ncv must satisfy nev < ncv <= n, n is the size of matrix");
263  }
264 
268  virtual ~SymEigsBase() {}
269 
271 
281  void init(const Scalar* init_resid)
282  {
283  // Reset all matrices/vectors to zero
284  m_ritz_val.resize(m_ncv);
285  m_ritz_vec.resize(m_ncv, m_nev);
286  m_ritz_est.resize(m_ncv);
287  m_ritz_conv.resize(m_nev);
288 
289  m_ritz_val.setZero();
290  m_ritz_vec.setZero();
291  m_ritz_est.setZero();
292  m_ritz_conv.setZero();
293 
294  m_nmatop = 0;
295  m_niter = 0;
296 
297  // Initialize the Lanczos factorization
298  MapConstVec v0(init_resid, m_n);
299  m_fac.init(v0, m_nmatop);
300  }
301 
309  void init()
310  {
311  SimpleRandom<Scalar> rng(0);
312  Vector init_resid = rng.random_vec(m_n);
313  init(init_resid.data());
314  }
315 
334  Index compute(Index maxit = 1000, Scalar tol = 1e-10, int sort_rule = LARGEST_ALGE)
335  {
336  // The m-step Lanczos factorization
337  m_fac.factorize_from(1, m_ncv, m_nmatop);
339  // Restarting
340  Index i, nconv = 0, nev_adj;
341  for (i = 0; i < maxit; i++)
342  {
343  nconv = num_converged(tol);
344  if (nconv >= m_nev)
345  break;
346 
347  nev_adj = nev_adjusted(nconv);
348  restart(nev_adj);
349  }
350  // Sorting results
351  sort_ritzpair(sort_rule);
352 
353  m_niter += i + 1;
354  m_info = (nconv >= m_nev) ? SUCCESSFUL : NOT_CONVERGING;
355 
356  return std::min(m_nev, nconv);
357  }
358 
363  int info() const { return m_info; }
364 
368  Index num_iterations() const { return m_niter; }
369 
373  Index num_operations() const { return m_nmatop; }
374 
382  Vector eigenvalues() const
383  {
384  const Index nconv = m_ritz_conv.cast<Index>().sum();
385  Vector res(nconv);
386 
387  if (!nconv)
388  return res;
389 
390  Index j = 0;
391  for (Index i = 0; i < m_nev; i++)
392  {
393  if (m_ritz_conv[i])
394  {
395  res[j] = m_ritz_val[i];
396  j++;
397  }
398  }
399 
400  return res;
401  }
402 
412  virtual Matrix eigenvectors(Index nvec) const
413  {
414  const Index nconv = m_ritz_conv.cast<Index>().sum();
415  nvec = std::min(nvec, nconv);
416  Matrix res(m_n, nvec);
417 
418  if (!nvec)
419  return res;
420 
421  Matrix ritz_vec_conv(m_ncv, nvec);
422  Index j = 0;
423  for (Index i = 0; i < m_nev && j < nvec; i++)
424  {
425  if (m_ritz_conv[i])
426  {
427  ritz_vec_conv.col(j).noalias() = m_ritz_vec.col(i);
428  j++;
429  }
430  }
431 
432  res.noalias() = m_fac.matrix_V() * ritz_vec_conv;
433 
434  return res;
435  }
436 
440  virtual Matrix eigenvectors() const
441  {
442  return eigenvectors(m_nev);
443  }
444 };
445 
446 } // namespace Spectra
447 
448 #endif // SYM_EIGS_BASE_H
SCALAR Scalar
Definition: bench_gemm.cpp:46
ArnoldiOp< Scalar, OpType, BOpType > ArnoldiOpType
Definition: SymEigsBase.h:56
EIGEN_DEVICE_FUNC Derived & setZero(Index size)
Vector eigenvalues() const
Definition: SymEigsBase.h:382
Eigen::Index Index
Definition: SymEigsBase.h:47
void restart(Index k)
Definition: SymEigsBase.h:85
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar * data() const
#define min(a, b)
Definition: datatypes.h:19
Index num_converged(Scalar tol)
Definition: SymEigsBase.h:113
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
void compute(ConstGenericMatrix &mat, const Scalar &shift=Scalar(0))
const Index m_ncv
Definition: SymEigsBase.h:65
const Vector & eigenvalues() const
Definition: TridiagEigen.h:199
static std::mt19937 rng
Eigen::Array< Scalar, Eigen::Dynamic, 1 > Array
Definition: SymEigsBase.h:50
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void swap(DenseBase< OtherDerived > &other)
const Matrix & matrix_V() const
Definition: Arnoldi.h:96
Computation was successful.
Definition: CompInfo.h:19
void init(MapConstVec &v0, Index &op_counter)
Definition: Arnoldi.h:103
Index num_operations() const
Definition: SymEigsBase.h:373
Select eigenvalues with smallest algebraic value. Only for symmetric eigen solvers.
Definition: SelectionRule.h:51
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > Matrix
Definition: SymEigsBase.h:48
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
virtual void sort_ritzpair(int sort_rule)
Definition: SymEigsBase.h:194
Eigen::Map< Vector > MapVec
Definition: SymEigsBase.h:53
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
virtual Matrix eigenvectors() const
Definition: SymEigsBase.h:440
BoolArray m_ritz_conv
Definition: SymEigsBase.h:75
void compress_H(const TridiagQR< Scalar > &decomp)
Definition: Lanczos.h:158
const Scalar m_eps
Definition: SymEigsBase.h:80
const Index m_nev
Definition: SymEigsBase.h:64
std::vector< int > ind
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Array< double, 1, 3 > e(1./3., 0.5, 2.)
Index num_iterations() const
Definition: SymEigsBase.h:368
const Scalar m_near_0
Definition: SymEigsBase.h:78
const Matrix & matrix_H() const
Definition: Arnoldi.h:97
Eigen::Map< Matrix > MapMat
Definition: SymEigsBase.h:52
Index nev_adjusted(Index nconv)
Definition: SymEigsBase.h:125
void compress_V(const Matrix &Q)
Definition: Arnoldi.h:264
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > Vector
Definition: SymEigsBase.h:49
Index compute(Index maxit=1000, Scalar tol=1e-10, int sort_rule=LARGEST_ALGE)
Definition: SymEigsBase.h:334
Eigen::Array< bool, Eigen::Dynamic, 1 > BoolArray
Definition: SymEigsBase.h:51
EIGEN_DEVICE_FUNC internal::pow_impl< ScalarX, ScalarY >::result_type pow(const ScalarX &x, const ScalarY &y)
Eigen::Map< const Vector > MapConstVec
Definition: SymEigsBase.h:54
void init(const Scalar *init_resid)
Definition: SymEigsBase.h:281
The quaternion class used to represent 3D orientations and rotations.
static const double v0
Scalar f_norm() const
Definition: Arnoldi.h:99
virtual Matrix eigenvectors(Index nvec) const
Definition: SymEigsBase.h:412
const G double tol
Definition: Group.h:86
void apply_YQ(GenericMatrix Y) const
const Scalar m_eps23
Definition: SymEigsBase.h:81
#define abs(x)
Definition: datatypes.h:17
std::ptrdiff_t j
void factorize_from(Index from_k, Index to_m, Index &op_counter)
Definition: Lanczos.h:55
Lanczos< Scalar, ArnoldiOpType > LanczosFac
Definition: SymEigsBase.h:57


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:36:35