Spectra
GenEigsSolver.h
1 // Copyright (C) 2016-2017 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 http://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
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/SelectionRule.h"
18 #include "Util/CompInfo.h"
19 #include "Util/SimpleRandom.h"
20 #include "LinAlg/UpperHessenbergQR.h"
21 #include "LinAlg/UpperHessenbergEigen.h"
22 #include "LinAlg/DoubleShiftQR.h"
23 #include "MatOp/DenseGenMatProd.h"
24 
25 namespace Spectra {
26 
27 
134 template < typename Scalar = double,
135  int SelectionRule = LARGEST_MAGN,
136  typename OpType = DenseGenMatProd<double> >
138 {
139 private:
140  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Matrix;
141  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;
142  typedef Eigen::Array<Scalar, Eigen::Dynamic, 1> Array;
143  typedef Eigen::Array<bool, Eigen::Dynamic, 1> BoolArray;
144  typedef Eigen::Map<Matrix> MapMat;
145  typedef Eigen::Map<Vector> MapVec;
146 
147  typedef std::complex<Scalar> Complex;
148  typedef Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic> ComplexMatrix;
149  typedef Eigen::Matrix<Complex, Eigen::Dynamic, 1> ComplexVector;
150 
151 protected:
152  OpType *m_op; // object to conduct matrix operation,
153  // e.g. matrix-vector product
154  const int m_n; // dimension of matrix A
155  const int m_nev; // number of eigenvalues requested
156 
157 private:
158  const int m_ncv; // number of ritz values
159  int m_nmatop; // number of matrix operations called
160  int m_niter; // number of restarting iterations
161 
162 protected:
163  Matrix m_fac_V; // V matrix in the Arnoldi factorization
164  Matrix m_fac_H; // H matrix in the Arnoldi factorization
165  Vector m_fac_f; // residual in the Arnoldi factorization
166 
167  ComplexVector m_ritz_val; // ritz values
168  ComplexMatrix m_ritz_vec; // ritz vectors
169  ComplexVector m_ritz_est; // last row of m_ritz_vec
170 
171 private:
172  BoolArray m_ritz_conv; // indicator of the convergence of ritz values
173  int m_info; // status of the computation
174 
175  const Scalar m_eps; // the machine precision,
176  // e.g. ~= 1e-16 for the "double" type
177  const Scalar m_approx_0; // a number that is approximately zero
178  // m_eps23 = m_eps^(2/3)
179  // used to test whether a number is complex, and
180  // to test the orthogonality of vectors
181 
182  // Arnoldi factorization starting from step-k
183  void factorize_from(int from_k, int to_m, const Vector& fk)
184  {
185  if(to_m <= from_k) return;
186 
187  m_fac_f = fk;
188 
189  Vector w(m_n);
190  Scalar beta = m_fac_f.norm();
191  // Keep the upperleft k x k submatrix of H and set other elements to 0
192  m_fac_H.rightCols(m_ncv - from_k).setZero();
193  m_fac_H.block(from_k, 0, m_ncv - from_k, from_k).setZero();
194  for(int i = from_k; i <= to_m - 1; i++)
195  {
196  bool restart = false;
197  // If beta = 0, then the next V is not full rank
198  // We need to generate a new residual vector that is orthogonal
199  // to the current V, which we call a restart
200  if(beta < m_eps)
201  {
202  SimpleRandom<Scalar> rng(2 * i);
203  m_fac_f.noalias() = rng.random_vec(m_n);
204  // f <- f - V * V' * f, so that f is orthogonal to V
205  MapMat V(m_fac_V.data(), m_n, i); // The first i columns
206  Vector Vf = V.transpose() * m_fac_f;
207  m_fac_f.noalias() -= V * Vf;
208  // beta <- ||f||
209  beta = m_fac_f.norm();
210 
211  restart = true;
212  }
213 
214  // v <- f / ||f||
215  m_fac_V.col(i).noalias() = m_fac_f / beta; // The (i+1)-th column
216 
217  // Note that H[i+1, i] equals to the unrestarted beta
218  if(restart)
219  m_fac_H(i, i - 1) = 0.0;
220  else
221  m_fac_H(i, i - 1) = beta;
222 
223  // w <- A * v, v = m_fac_V.col(i)
224  m_op->perform_op(&m_fac_V(0, i), w.data());
225  m_nmatop++;
226 
227  // First i+1 columns of V
228  MapMat Vs(m_fac_V.data(), m_n, i + 1);
229  // h = m_fac_H(0:i, i)
230  MapVec h(&m_fac_H(0, i), i + 1);
231  // h <- V' * w
232  h.noalias() = Vs.transpose() * w;
233 
234  // f <- w - V * h
235  m_fac_f.noalias() = w - Vs * h;
236  beta = m_fac_f.norm();
237 
238  if(beta > 0.717 * h.norm())
239  continue;
240 
241  // f/||f|| is going to be the next column of V, so we need to test
242  // whether V' * (f/||f||) ~= 0
243  Vector Vf = Vs.transpose() * m_fac_f;
244  // If not, iteratively correct the residual
245  int count = 0;
246  while(count < 5 && Vf.cwiseAbs().maxCoeff() > m_approx_0 * beta)
247  {
248  // f <- f - V * Vf
249  m_fac_f.noalias() -= Vs * Vf;
250  // h <- h + Vf
251  h.noalias() += Vf;
252  // beta <- ||f||
253  beta = m_fac_f.norm();
254 
255  Vf.noalias() = Vs.transpose() * m_fac_f;
256  count++;
257  }
258  }
259  }
260 
261  // Real Ritz values calculated from UpperHessenbergEigen have exact zero imaginary part
262  // Complex Ritz values have exact conjugate pairs
263  // So we use exact tests here
264  static bool is_complex(const Complex& v) { return v.imag() != Scalar(0); }
265  static bool is_conj(const Complex& v1, const Complex& v2) { return v1 == Eigen::numext::conj(v2); }
266 
267  // Implicitly restarted Arnoldi factorization
268  void restart(int k)
269  {
270  using std::norm;
271 
272  if(k >= m_ncv)
273  return;
274 
275  DoubleShiftQR<Scalar> decomp_ds(m_ncv);
276  UpperHessenbergQR<Scalar> decomp_hb;
277  Matrix Q = Matrix::Identity(m_ncv, m_ncv);
278 
279  for(int i = k; i < m_ncv; i++)
280  {
281  if(is_complex(m_ritz_val[i]) && is_conj(m_ritz_val[i], m_ritz_val[i + 1]))
282  {
283  // H - mu * I = Q1 * R1
284  // H <- R1 * Q1 + mu * I = Q1' * H * Q1
285  // H - conj(mu) * I = Q2 * R2
286  // H <- R2 * Q2 + conj(mu) * I = Q2' * H * Q2
287  //
288  // (H - mu * I) * (H - conj(mu) * I) = Q1 * Q2 * R2 * R1 = Q * R
289  Scalar s = 2 * m_ritz_val[i].real();
290  Scalar t = norm(m_ritz_val[i]);
291 
292  decomp_ds.compute(m_fac_H, s, t);
293 
294  // Q -> Q * Qi
295  decomp_ds.apply_YQ(Q);
296  // H -> Q'HQ
297  // Matrix Q = Matrix::Identity(m_ncv, m_ncv);
298  // decomp_ds.apply_YQ(Q);
299  // m_fac_H = Q.transpose() * m_fac_H * Q;
300  m_fac_H = decomp_ds.matrix_QtHQ();
301 
302  i++;
303  } else {
304  // QR decomposition of H - mu * I, mu is real
305  m_fac_H.diagonal().array() -= m_ritz_val[i].real();
306  decomp_hb.compute(m_fac_H);
307 
308  // Q -> Q * Qi
309  decomp_hb.apply_YQ(Q);
310  // H -> Q'HQ = RQ + mu * I
311  m_fac_H = decomp_hb.matrix_RQ();
312  m_fac_H.diagonal().array() += m_ritz_val[i].real();
313  }
314  }
315  // V -> VQ, only need to update the first k+1 columns
316  // Q has some elements being zero
317  // The first (ncv - k + i) elements of the i-th column of Q are non-zero
318  Matrix Vs(m_n, k + 1);
319  int nnz;
320  for(int i = 0; i < k; i++)
321  {
322  nnz = m_ncv - k + i + 1;
323  MapMat V(m_fac_V.data(), m_n, nnz);
324  MapVec q(&Q(0, i), nnz);
325  Vs.col(i).noalias() = V * q;
326  }
327  Vs.col(k).noalias() = m_fac_V * Q.col(k);
328  m_fac_V.leftCols(k + 1).noalias() = Vs;
329 
330  Vector fk = m_fac_f * Q(m_ncv - 1, k - 1) + m_fac_V.col(k) * m_fac_H(k, k - 1);
331  factorize_from(k, m_ncv, fk);
332  retrieve_ritzpair();
333  }
334 
335  // Calculates the number of converged Ritz values
336  int num_converged(Scalar tol)
337  {
338  // thresh = tol * max(m_approx_0, abs(theta)), theta for ritz value
339  Array thresh = tol * m_ritz_val.head(m_nev).array().abs().max(m_approx_0);
340  Array resid = m_ritz_est.head(m_nev).array().abs() * m_fac_f.norm();
341  // Converged "wanted" ritz values
342  m_ritz_conv = (resid < thresh);
343 
344  return m_ritz_conv.cast<int>().sum();
345  }
346 
347  // Returns the adjusted nev for restarting
348  int nev_adjusted(int nconv)
349  {
350  using std::abs;
351 
352  int nev_new = m_nev;
353  for(int i = m_nev; i < m_ncv; i++)
354  if(abs(m_ritz_est[i]) < m_eps) nev_new++;
355 
356  // Adjust nev_new, according to dnaup2.f line 660~674 in ARPACK
357  nev_new += std::min(nconv, (m_ncv - nev_new) / 2);
358  if(nev_new == 1 && m_ncv >= 6)
359  nev_new = m_ncv / 2;
360  else if(nev_new == 1 && m_ncv > 3)
361  nev_new = 2;
362 
363  if(nev_new > m_ncv - 2)
364  nev_new = m_ncv - 2;
365 
366  // Increase nev by one if ritz_val[nev - 1] and
367  // ritz_val[nev] are conjugate pairs
368  if(is_complex(m_ritz_val[nev_new - 1]) &&
369  is_conj(m_ritz_val[nev_new - 1], m_ritz_val[nev_new]))
370  {
371  nev_new++;
372  }
373 
374  return nev_new;
375  }
376 
377  // Retrieves and sorts ritz values and ritz vectors
378  void retrieve_ritzpair()
379  {
380  UpperHessenbergEigen<Scalar> decomp(m_fac_H);
381  ComplexVector evals = decomp.eigenvalues();
382  ComplexMatrix evecs = decomp.eigenvectors();
383 
384  SortEigenvalue<Complex, SelectionRule> sorting(evals.data(), evals.size());
385  std::vector<int> ind = sorting.index();
386 
387  // Copy the ritz values and vectors to m_ritz_val and m_ritz_vec, respectively
388  for(int i = 0; i < m_ncv; i++)
389  {
390  m_ritz_val[i] = evals[ind[i]];
391  m_ritz_est[i] = evecs(m_ncv - 1, ind[i]);
392  }
393  for(int i = 0; i < m_nev; i++)
394  {
395  m_ritz_vec.col(i) = evecs.col(ind[i]);
396  }
397  }
398 
399 protected:
400  // Sorts the first nev Ritz pairs in the specified order
401  // This is used to return the final results
402  virtual void sort_ritzpair(int sort_rule)
403  {
404  // First make sure that we have a valid index vector
405  SortEigenvalue<Complex, LARGEST_MAGN> sorting(m_ritz_val.data(), m_nev);
406  std::vector<int> ind = sorting.index();
407 
408  switch(sort_rule)
409  {
410  case LARGEST_MAGN:
411  break;
412  case LARGEST_REAL:
413  {
414  SortEigenvalue<Complex, LARGEST_REAL> sorting(m_ritz_val.data(), m_nev);
415  ind = sorting.index();
416  }
417  break;
418  case LARGEST_IMAG:
419  {
420  SortEigenvalue<Complex, LARGEST_IMAG> sorting(m_ritz_val.data(), m_nev);
421  ind = sorting.index();
422  }
423  break;
424  case SMALLEST_MAGN:
425  {
426  SortEigenvalue<Complex, SMALLEST_MAGN> sorting(m_ritz_val.data(), m_nev);
427  ind = sorting.index();
428  }
429  break;
430  case SMALLEST_REAL:
431  {
432  SortEigenvalue<Complex, SMALLEST_REAL> sorting(m_ritz_val.data(), m_nev);
433  ind = sorting.index();
434  }
435  break;
436  case SMALLEST_IMAG:
437  {
438  SortEigenvalue<Complex, SMALLEST_IMAG> sorting(m_ritz_val.data(), m_nev);
439  ind = sorting.index();
440  }
441  break;
442  default:
443  throw std::invalid_argument("unsupported sorting rule");
444  }
445 
446  ComplexVector new_ritz_val(m_ncv);
447  ComplexMatrix new_ritz_vec(m_ncv, m_nev);
448  BoolArray new_ritz_conv(m_nev);
449 
450  for(int i = 0; i < m_nev; i++)
451  {
452  new_ritz_val[i] = m_ritz_val[ind[i]];
453  new_ritz_vec.col(i) = m_ritz_vec.col(ind[i]);
454  new_ritz_conv[i] = m_ritz_conv[ind[i]];
455  }
456 
457  m_ritz_val.swap(new_ritz_val);
458  m_ritz_vec.swap(new_ritz_vec);
459  m_ritz_conv.swap(new_ritz_conv);
460  }
461 
462 public:
480  GenEigsSolver(OpType* op_, int nev_, int ncv_) :
481  m_op(op_),
482  m_n(m_op->rows()),
483  m_nev(nev_),
484  m_ncv(ncv_ > m_n ? m_n : ncv_),
485  m_nmatop(0),
486  m_niter(0),
487  m_info(NOT_COMPUTED),
488  m_eps(Eigen::NumTraits<Scalar>::epsilon()),
489  m_approx_0(Eigen::numext::pow(m_eps, Scalar(2.0) / 3))
490  {
491  if(nev_ < 1 || nev_ > m_n - 2)
492  throw std::invalid_argument("nev must satisfy 1 <= nev <= n - 2, n is the size of matrix");
493 
494  if(ncv_ < nev_ + 2 || ncv_ > m_n)
495  throw std::invalid_argument("ncv must satisfy nev + 2 <= ncv <= n, n is the size of matrix");
496  }
497 
501  virtual ~GenEigsSolver() {}
502 
512  void init(const Scalar* init_resid)
513  {
514  // Reset all matrices/vectors to zero
515  m_fac_V.resize(m_n, m_ncv);
516  m_fac_H.resize(m_ncv, m_ncv);
517  m_fac_f.resize(m_n);
518  m_ritz_val.resize(m_ncv);
519  m_ritz_vec.resize(m_ncv, m_nev);
520  m_ritz_est.resize(m_ncv);
521  m_ritz_conv.resize(m_nev);
522 
523  m_fac_V.setZero();
524  m_fac_H.setZero();
525  m_fac_f.setZero();
526  m_ritz_val.setZero();
527  m_ritz_vec.setZero();
528  m_ritz_est.setZero();
529  m_ritz_conv.setZero();
530 
531  // Set the initial vector
532  Vector v(m_n);
533  std::copy(init_resid, init_resid + m_n, v.data());
534  Scalar vnorm = v.norm();
535  if(vnorm < m_eps)
536  throw std::invalid_argument("initial residual vector cannot be zero");
537  v /= vnorm;
538 
539  Vector w(m_n);
540  m_op->perform_op(v.data(), w.data());
541  m_nmatop++;
542 
543  m_fac_H(0, 0) = v.dot(w);
544  m_fac_f = w - v * m_fac_H(0, 0);
545  m_fac_V.col(0) = v;
546  }
547 
555  void init()
556  {
557  SimpleRandom<Scalar> rng(0);
558  Vector init_resid = rng.random_vec(m_n);
559  init(init_resid.data());
560  }
561 
582  int compute(int maxit = 1000, Scalar tol = 1e-10, int sort_rule = LARGEST_MAGN)
583  {
584  // The m-step Arnoldi factorization
585  factorize_from(1, m_ncv, m_fac_f);
586  retrieve_ritzpair();
587  // Restarting
588  int i, nconv = 0, nev_adj;
589  for(i = 0; i < maxit; i++)
590  {
591  nconv = num_converged(tol);
592  if(nconv >= m_nev)
593  break;
594 
595  nev_adj = nev_adjusted(nconv);
596  restart(nev_adj);
597  }
598  // Sorting results
599  sort_ritzpair(sort_rule);
600 
601  m_niter += i + 1;
602  m_info = (nconv >= m_nev) ? SUCCESSFUL : NOT_CONVERGING;
603 
604  return std::min(m_nev, nconv);
605  }
606 
611  int info() const { return m_info; }
612 
616  int num_iterations() const { return m_niter; }
617 
621  int num_operations() const { return m_nmatop; }
622 
630  ComplexVector eigenvalues() const
631  {
632  int nconv = m_ritz_conv.cast<int>().sum();
633  ComplexVector res(nconv);
634 
635  if(!nconv)
636  return res;
637 
638  int j = 0;
639  for(int i = 0; i < m_nev; i++)
640  {
641  if(m_ritz_conv[i])
642  {
643  res[j] = m_ritz_val[i];
644  j++;
645  }
646  }
647 
648  return res;
649  }
650 
660  ComplexMatrix eigenvectors(int nvec) const
661  {
662  int nconv = m_ritz_conv.cast<int>().sum();
663  nvec = std::min(nvec, nconv);
664  ComplexMatrix res(m_n, nvec);
665 
666  if(!nvec)
667  return res;
668 
669  ComplexMatrix ritz_vec_conv(m_ncv, nvec);
670  int j = 0;
671  for(int i = 0; i < m_nev && j < nvec; i++)
672  {
673  if(m_ritz_conv[i])
674  {
675  ritz_vec_conv.col(j) = m_ritz_vec.col(i);
676  j++;
677  }
678  }
679 
680  res.noalias() = m_fac_V * ritz_vec_conv;
681 
682  return res;
683  }
684 
688  ComplexMatrix eigenvectors() const
689  {
690  return eigenvectors(m_nev);
691  }
692 };
693 
694 
695 } // namespace Spectra
696 
697 #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
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
GenEigsSolver(OpType *op_, int nev_, int ncv_)
ComplexVector eigenvalues() const
int num_operations() const