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