Spectra
SymEigsSolver.h
1 // Copyright (C) 2016-2018 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_SOLVER_H
8 #define SYM_EIGS_SOLVER_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"
20 #include "LinAlg/UpperHessenbergQR.h"
21 #include "LinAlg/TridiagEigen.h"
22 #include "MatOp/DenseSymMatProd.h"
23 
24 
25 namespace Spectra {
26 
27 
33 
153 template < typename Scalar = double,
154  int SelectionRule = LARGEST_MAGN,
155  typename OpType = DenseSymMatProd<double> >
157 {
158 private:
159  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Matrix;
160  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;
161  typedef Eigen::Array<Scalar, Eigen::Dynamic, 1> Array;
162  typedef Eigen::Array<bool, Eigen::Dynamic, 1> BoolArray;
163  typedef Eigen::Map<Matrix> MapMat;
164  typedef Eigen::Map<Vector> MapVec;
165 
166 protected:
167  OpType* m_op; // object to conduct matrix operation,
168  // e.g. matrix-vector product
169  const int m_n; // dimension of matrix A
170  const int m_nev; // number of eigenvalues requested
171  const int m_ncv; // dimension of Krylov subspace in the Lanczos method
172  int m_nmatop; // number of matrix operations called
173  int m_niter; // number of restarting iterations
174 
175  Matrix m_fac_V; // V matrix in the Lanczos factorization
176  Matrix m_fac_H; // H matrix in the Lanczos factorization
177  Vector m_fac_f; // residual in the Lanczos factorization
178  Vector m_ritz_val; // Ritz values
179 
180 private:
181  Matrix m_ritz_vec; // Ritz vectors
182  Vector m_ritz_est; // last row of m_ritz_vec, also called the Ritz estimates
183  BoolArray m_ritz_conv; // indicator of the convergence of Ritz values
184  int m_info; // status of the computation
185 
186  const Scalar m_near_0; // a very small value, but 1.0 / m_near_0 does not overflow
187  // ~= 1e-307 for the "double" type
188  const Scalar m_eps; // the machine precision, ~= 1e-16 for the "double" type
189  const Scalar m_eps23; // m_eps^(2/3), used to test the convergence
190 
191  // Given orthonormal basis functions V, find a nonzero vector f such that V'f = 0
192  // Assume that f has been properly allocated
193  void expand_basis(const MapMat& V, const int seed, Vector& f, Scalar& fnorm)
194  {
195  using std::sqrt;
196 
197  const Scalar thresh = m_eps * sqrt(Scalar(m_n));
198  for(int iter = 0; iter < 5; iter++)
199  {
200  // Randomly generate a new vector and orthogonalize it against V
201  SimpleRandom<Scalar> rng(seed + 123 * iter);
202  f.noalias() = rng.random_vec(m_n);
203  // f <- f - V * V' * f, so that f is orthogonal to V
204  Vector Vf = inner_product(V, m_fac_f);
205  f -= V * Vf;
206  // fnorm <- ||f||
207  fnorm = m_fac_f.norm();
208 
209  // If fnorm is too close to zero, we try a new random vector,
210  // otherwise return the result
211  if(fnorm >= thresh)
212  return;
213  }
214  }
215 
216  // Lanczos factorization starting from step-k
217  void factorize_from(int from_k, int to_m, const Vector& fk)
218  {
219  using std::sqrt;
220 
221  if(to_m <= from_k) return;
222 
223  const Scalar beta_thresh = m_eps * sqrt(Scalar(m_n));
224  m_fac_f.noalias() = fk;
225 
226  // Pre-allocate Vf
227  Vector Vf(to_m);
228  Vector w(m_n);
229  Scalar beta = norm(m_fac_f), Hii = Scalar(0);
230  // Keep the upperleft k x k submatrix of H and set other elements to 0
231  m_fac_H.rightCols(m_ncv - from_k).setZero();
232  m_fac_H.block(from_k, 0, m_ncv - from_k, from_k).setZero();
233  for(int i = from_k; i <= to_m - 1; i++)
234  {
235  bool restart = false;
236  // If beta = 0, then the next V is not full rank
237  // We need to generate a new residual vector that is orthogonal
238  // to the current V, which we call a restart
239  if(beta < m_near_0)
240  {
241  MapMat V(m_fac_V.data(), m_n, i); // The first i columns
242  expand_basis(V, 2 * i, m_fac_f, beta);
243  restart = true;
244  }
245 
246  // v <- f / ||f||
247  MapVec v(&m_fac_V(0, i), m_n); // The (i+1)-th column
248  v.noalias() = m_fac_f / beta;
249 
250  // Note that H[i+1, i] equals to the unrestarted beta
251  m_fac_H(i, i - 1) = restart ? Scalar(0) : beta;
252 
253  // w <- A * v
254  m_op->perform_op(v.data(), w.data());
255  m_nmatop++;
256 
257  Hii = inner_product(v, w);
258  m_fac_H(i - 1, i) = m_fac_H(i, i - 1); // Due to symmetry
259  m_fac_H(i, i) = Hii;
260 
261  // f <- w - V * V' * w = w - H[i+1, i] * V{i} - H[i+1, i+1] * V{i+1}
262  // If restarting, we know that H[i+1, i] = 0
263  if(restart)
264  m_fac_f.noalias() = w - Hii * v;
265  else
266  m_fac_f.noalias() = w - m_fac_H(i, i - 1) * m_fac_V.col(i - 1) - Hii * v;
267 
268  beta = norm(m_fac_f);
269 
270  // f/||f|| is going to be the next column of V, so we need to test
271  // whether V' * (f/||f||) ~= 0
272  const int i1 = i + 1;
273  MapMat V(m_fac_V.data(), m_n, i1); // The first (i+1) columns
274  Vf.head(i1).noalias() = inner_product(V, m_fac_f);
275  Scalar ortho_err = Vf.head(i1).cwiseAbs().maxCoeff();
276  // If not, iteratively correct the residual
277  int count = 0;
278  while(count < 5 && ortho_err > m_eps * beta)
279  {
280  // There is an edge case: when beta=||f|| is close to zero, f mostly consists
281  // of noises of rounding errors, so the test [ortho_err < eps * beta] is very
282  // likely to fail. In particular, if beta=0, then the test is ensured to fail.
283  // Hence when this happens, we force f to be zero, and then restart in the
284  // next iteration.
285  if(beta < beta_thresh)
286  {
287  m_fac_f.setZero();
288  beta = Scalar(0);
289  break;
290  }
291 
292  // f <- f - V * Vf
293  m_fac_f.noalias() -= V * Vf.head(i1);
294  // h <- h + Vf
295  m_fac_H(i - 1, i) += Vf[i - 1];
296  m_fac_H(i, i - 1) = m_fac_H(i - 1, i);
297  m_fac_H(i, i) += Vf[i];
298  // beta <- ||f||
299  beta = norm(m_fac_f);
300 
301  Vf.head(i1).noalias() = inner_product(V, m_fac_f);
302  ortho_err = Vf.head(i1).cwiseAbs().maxCoeff();
303  count++;
304  }
305  }
306  }
307 
308  // Implicitly restarted Lanczos factorization
309  void restart(int k)
310  {
311  if(k >= m_ncv)
312  return;
313 
314  TridiagQR<Scalar> decomp;
315  Matrix Q = Matrix::Identity(m_ncv, m_ncv);
316 
317  for(int i = k; i < m_ncv; i++)
318  {
319  // QR decomposition of H-mu*I, mu is the shift
320  m_fac_H.diagonal().array() -= m_ritz_val[i];
321  decomp.compute(m_fac_H);
322 
323  // Q -> Q * Qi
324  decomp.apply_YQ(Q);
325  // H -> Q'HQ
326  // Since QR = H - mu * I, we have H = QR + mu * I
327  // and therefore Q'HQ = RQ + mu * I
328  decomp.matrix_RQ(m_fac_H);
329  m_fac_H.diagonal().array() += m_ritz_val[i];
330  }
331  // V -> VQ, only need to update the first k+1 columns
332  // Q has some elements being zero
333  // The first (ncv - k + i) elements of the i-th column of Q are non-zero
334  Matrix Vs(m_n, k + 1);
335  for(int i = 0; i < k; i++)
336  {
337  const int nnz = m_ncv - k + i + 1;
338  MapVec q(&Q(0, i), nnz);
339  Vs.col(i).noalias() = m_fac_V.leftCols(nnz) * q;
340  }
341  Vs.col(k).noalias() = m_fac_V * Q.col(k);
342  m_fac_V.leftCols(k + 1).noalias() = Vs;
343 
344  const Vector fk = m_fac_f * Q(m_ncv - 1, k - 1) + m_fac_V.col(k) * m_fac_H(k, k - 1);
345  factorize_from(k, m_ncv, fk);
346  retrieve_ritzpair();
347  }
348 
349  // Calculates the number of converged Ritz values
350  int num_converged(Scalar tol)
351  {
352  // thresh = tol * max(m_eps23, abs(theta)), theta for Ritz value
353  Array thresh = tol * m_ritz_val.head(m_nev).array().abs().max(m_eps23);
354  Array resid = m_ritz_est.head(m_nev).array().abs() * norm(m_fac_f);
355  // Converged "wanted" Ritz values
356  m_ritz_conv = (resid < thresh);
357 
358  return m_ritz_conv.cast<int>().sum();
359  }
360 
361  // Returns the adjusted nev for restarting
362  int nev_adjusted(int nconv)
363  {
364  using std::abs;
365 
366  int nev_new = m_nev;
367  for(int i = m_nev; i < m_ncv; i++)
368  if(abs(m_ritz_est[i]) < m_near_0) nev_new++;
369 
370  // Adjust nev_new, according to dsaup2.f line 677~684 in ARPACK
371  nev_new += std::min(nconv, (m_ncv - nev_new) / 2);
372  if(nev_new == 1 && m_ncv >= 6)
373  nev_new = m_ncv / 2;
374  else if(nev_new == 1 && m_ncv > 2)
375  nev_new = 2;
376 
377  if(nev_new > m_ncv - 1)
378  nev_new = m_ncv - 1;
379 
380  return nev_new;
381  }
382 
383  // Retrieves and sorts Ritz values and Ritz vectors
384  void retrieve_ritzpair()
385  {
386  TridiagEigen<Scalar> decomp(m_fac_H);
387  const Vector& evals = decomp.eigenvalues();
388  const Matrix& evecs = decomp.eigenvectors();
389 
390  SortEigenvalue<Scalar, SelectionRule> sorting(evals.data(), evals.size());
391  std::vector<int> ind = sorting.index();
392 
393  // For BOTH_ENDS, the eigenvalues are sorted according
394  // to the LARGEST_ALGE rule, so we need to move those smallest
395  // values to the left
396  // The order would be
397  // Largest => Smallest => 2nd largest => 2nd smallest => ...
398  // We keep this order since the first k values will always be
399  // the wanted collection, no matter k is nev_updated (used in restart())
400  // or is nev (used in sort_ritzpair())
401  if(SelectionRule == BOTH_ENDS)
402  {
403  std::vector<int> ind_copy(ind);
404  for(int i = 0; i < m_ncv; i++)
405  {
406  // If i is even, pick values from the left (large values)
407  // If i is odd, pick values from the right (small values)
408  if(i % 2 == 0)
409  ind[i] = ind_copy[i / 2];
410  else
411  ind[i] = ind_copy[m_ncv - 1 - i / 2];
412  }
413  }
414 
415  // Copy the Ritz values and vectors to m_ritz_val and m_ritz_vec, respectively
416  for(int i = 0; i < m_ncv; i++)
417  {
418  m_ritz_val[i] = evals[ind[i]];
419  m_ritz_est[i] = evecs(m_ncv - 1, ind[i]);
420  }
421  for(int i = 0; i < m_nev; i++)
422  {
423  m_ritz_vec.col(i).noalias() = evecs.col(ind[i]);
424  }
425  }
426 
427 protected:
428  // In generalized eigenvalue problem Ax=lambda*Bx, define the inner product to be <x, y> = x'By
429  // For regular eigenvalue problems, it is the usual inner product <x, y> = x'y
430  virtual Scalar inner_product(const Vector& x, const Vector& y) { return x.dot(y); }
431  virtual Scalar inner_product(const MapVec& x, const Vector& y) { return x.dot(y); }
432  virtual Vector inner_product(const MapMat& x, const Vector& y) { return x.transpose() * y; }
433 
434  // B-norm of a vector. For regular eigenvalue problems it is simply the L2 norm
435  virtual Scalar norm(const Vector& x) { return x.norm(); }
436 
437  // Sorts the first nev Ritz pairs in the specified order
438  // This is used to return the final results
439  virtual void sort_ritzpair(int sort_rule)
440  {
441  // First make sure that we have a valid index vector
442  SortEigenvalue<Scalar, LARGEST_ALGE> sorting(m_ritz_val.data(), m_nev);
443  std::vector<int> ind = sorting.index();
444 
445  switch(sort_rule)
446  {
447  case LARGEST_ALGE:
448  break;
449  case LARGEST_MAGN:
450  {
451  SortEigenvalue<Scalar, LARGEST_MAGN> sorting(m_ritz_val.data(), m_nev);
452  ind = sorting.index();
453  }
454  break;
455  case SMALLEST_ALGE:
456  {
457  SortEigenvalue<Scalar, SMALLEST_ALGE> sorting(m_ritz_val.data(), m_nev);
458  ind = sorting.index();
459  }
460  break;
461  case SMALLEST_MAGN:
462  {
463  SortEigenvalue<Scalar, SMALLEST_MAGN> sorting(m_ritz_val.data(), m_nev);
464  ind = sorting.index();
465  }
466  break;
467  default:
468  throw std::invalid_argument("unsupported sorting rule");
469  }
470 
471  Vector new_ritz_val(m_ncv);
472  Matrix new_ritz_vec(m_ncv, m_nev);
473  BoolArray new_ritz_conv(m_nev);
474 
475  for(int i = 0; i < m_nev; i++)
476  {
477  new_ritz_val[i] = m_ritz_val[ind[i]];
478  new_ritz_vec.col(i).noalias() = m_ritz_vec.col(ind[i]);
479  new_ritz_conv[i] = m_ritz_conv[ind[i]];
480  }
481 
482  m_ritz_val.swap(new_ritz_val);
483  m_ritz_vec.swap(new_ritz_vec);
484  m_ritz_conv.swap(new_ritz_conv);
485  }
486 
487 public:
505  SymEigsSolver(OpType* op, int nev, int ncv) :
506  m_op(op),
507  m_n(m_op->rows()),
508  m_nev(nev),
509  m_ncv(ncv > m_n ? m_n : ncv),
510  m_nmatop(0),
511  m_niter(0),
512  m_info(NOT_COMPUTED),
513  m_near_0(TypeTraits<Scalar>::min() * Scalar(10)),
514  m_eps(Eigen::NumTraits<Scalar>::epsilon()),
515  m_eps23(Eigen::numext::pow(m_eps, Scalar(2.0) / 3))
516  {
517  if(nev < 1 || nev > m_n - 1)
518  throw std::invalid_argument("nev must satisfy 1 <= nev <= n - 1, n is the size of matrix");
519 
520  if(ncv <= nev || ncv > m_n)
521  throw std::invalid_argument("ncv must satisfy nev < ncv <= n, n is the size of matrix");
522  }
523 
527  virtual ~SymEigsSolver() {}
528 
538  void init(const Scalar* init_resid)
539  {
540  // Reset all matrices/vectors to zero
541  m_fac_V.resize(m_n, m_ncv);
542  m_fac_H.resize(m_ncv, m_ncv);
543  m_fac_f.resize(m_n);
544  m_ritz_val.resize(m_ncv);
545  m_ritz_vec.resize(m_ncv, m_nev);
546  m_ritz_est.resize(m_ncv);
547  m_ritz_conv.resize(m_nev);
548 
549  m_fac_V.setZero();
550  m_fac_H.setZero();
551  m_fac_f.setZero();
552  m_ritz_val.setZero();
553  m_ritz_vec.setZero();
554  m_ritz_est.setZero();
555  m_ritz_conv.setZero();
556 
557  m_nmatop = 0;
558  m_niter = 0;
559 
560  // Set the initial vector
561  Vector v(m_n);
562  std::copy(init_resid, init_resid + m_n, v.data());
563  const Scalar vnorm = norm(v);
564  if(vnorm < m_near_0)
565  throw std::invalid_argument("initial residual vector cannot be zero");
566  v /= vnorm;
567 
568  Vector w(m_n);
569  m_op->perform_op(v.data(), w.data());
570  m_nmatop++;
571 
572  m_fac_H(0, 0) = inner_product(v, w);
573  m_fac_f.noalias() = w - v * m_fac_H(0, 0);
574  m_fac_V.col(0).noalias() = v;
575 
576  // In some cases f is zero in exact arithmetics, but due to rounding errors
577  // it may contain tiny fluctuations. When this happens, we force f to be zero.
578  if(m_fac_f.cwiseAbs().maxCoeff() < m_eps)
579  m_fac_f.setZero();
580  }
581 
589  void init()
590  {
591  SimpleRandom<Scalar> rng(0);
592  Vector init_resid = rng.random_vec(m_n);
593  init(init_resid.data());
594  }
595 
614  int compute(int maxit = 1000, Scalar tol = 1e-10, int sort_rule = LARGEST_ALGE)
615  {
616  // The m-step Arnoldi factorization
617  factorize_from(1, m_ncv, m_fac_f);
618  retrieve_ritzpair();
619  // Restarting
620  int i, nconv = 0, nev_adj;
621  for(i = 0; i < maxit; i++)
622  {
623  nconv = num_converged(tol);
624  if(nconv >= m_nev)
625  break;
626 
627  nev_adj = nev_adjusted(nconv);
628  restart(nev_adj);
629  }
630  // Sorting results
631  sort_ritzpair(sort_rule);
632 
633  m_niter += i + 1;
634  m_info = (nconv >= m_nev) ? SUCCESSFUL : NOT_CONVERGING;
635 
636  return std::min(m_nev, nconv);
637  }
638 
643  int info() const { return m_info; }
644 
648  int num_iterations() const { return m_niter; }
649 
653  int num_operations() const { return m_nmatop; }
654 
662  Vector eigenvalues() const
663  {
664  const int nconv = m_ritz_conv.cast<int>().sum();
665  Vector res(nconv);
666 
667  if(!nconv)
668  return res;
669 
670  int j = 0;
671  for(int i = 0; i < m_nev; i++)
672  {
673  if(m_ritz_conv[i])
674  {
675  res[j] = m_ritz_val[i];
676  j++;
677  }
678  }
679 
680  return res;
681  }
682 
692  virtual Matrix eigenvectors(int nvec) const
693  {
694  const int nconv = m_ritz_conv.cast<int>().sum();
695  nvec = std::min(nvec, nconv);
696  Matrix res(m_n, nvec);
697 
698  if(!nvec)
699  return res;
700 
701  Matrix ritz_vec_conv(m_ncv, nvec);
702  int j = 0;
703  for(int i = 0; i < m_nev && j < nvec; i++)
704  {
705  if(m_ritz_conv[i])
706  {
707  ritz_vec_conv.col(j).noalias() = m_ritz_vec.col(i);
708  j++;
709  }
710  }
711 
712  res.noalias() = m_fac_V * ritz_vec_conv;
713 
714  return res;
715  }
716 
720  virtual Matrix eigenvectors() const
721  {
722  return eigenvectors(m_nev);
723  }
724 };
725 
726 
727 } // namespace Spectra
728 
729 #endif // SYM_EIGS_SOLVER_H
int num_operations() const
int compute(int maxit=1000, Scalar tol=1e-10, int sort_rule=LARGEST_ALGE)
Computation was successful.
Definition: CompInfo.h:20
virtual Matrix eigenvectors() const
Select eigenvalues with smallest algebraic value. Only for symmetric eigen solvers.
Definition: SelectionRule.h:52
int num_iterations() const
void init(const Scalar *init_resid)
Vector eigenvalues() const
SymEigsSolver(OpType *op, int nev, int ncv)
virtual Matrix eigenvectors(int nvec) const