Spectra
SymEigsSolver.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 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
13 #include <algorithm> // std::min, std::copy
14 #include <stdexcept> // std::invalid_argument
15 
16 #include "Util/SelectionRule.h"
17 #include "Util/CompInfo.h"
18 #include "Util/SimpleRandom.h"
19 #include "LinAlg/UpperHessenbergQR.h"
20 #include "LinAlg/TridiagEigen.h"
21 #include "MatOp/DenseSymMatProd.h"
22 
23 
24 namespace Spectra {
25 
26 
32 
152 template < typename Scalar = double,
153  int SelectionRule = LARGEST_MAGN,
154  typename OpType = DenseSymMatProd<double> >
156 {
157 private:
158  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Matrix;
159  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;
160  typedef Eigen::Array<Scalar, Eigen::Dynamic, 1> Array;
161  typedef Eigen::Array<bool, Eigen::Dynamic, 1> BoolArray;
162  typedef Eigen::Map<Matrix> MapMat;
163  typedef Eigen::Map<Vector> MapVec;
164 
165 protected:
166  OpType* m_op; // object to conduct matrix operation,
167  // e.g. matrix-vector product
168 
169 private:
170  const int m_n; // dimension of matrix A
171 
172 protected:
173  const int m_nev; // number of eigenvalues requested
174 
175 private:
176  const int m_ncv; // number of ritz values
177  int m_nmatop; // number of matrix operations called
178  int m_niter; // number of restarting iterations
179 
180  Matrix m_fac_V; // V matrix in the Arnoldi factorization
181  Matrix m_fac_H; // H matrix in the Arnoldi factorization
182  Vector m_fac_f; // residual in the Arnoldi factorization
183 
184 protected:
185  Vector m_ritz_val; // ritz values
186 
187 private:
188  Matrix m_ritz_vec; // ritz vectors
189  Vector m_ritz_est; // last row of m_ritz_vec
190  BoolArray m_ritz_conv; // indicator of the convergence of ritz values
191  int m_info; // status of the computation
192 
193  const Scalar m_eps; // the machine precision,
194  // e.g. ~= 1e-16 for the "double" type
195  const Scalar m_approx_0; // a number that is approximately zero
196  // m_approx_0 = m_eps^(2/3)
197  // used to test the orthogonality of vectors
198 
199  // Arnoldi factorization starting from step-k
200  void factorize_from(int from_k, int to_m, const Vector& fk)
201  {
202  if(to_m <= from_k) return;
203 
204  m_fac_f.noalias() = fk;
205 
206  Vector w(m_n);
207  Scalar beta = norm(m_fac_f), Hii = 0.0;
208  // Keep the upperleft k x k submatrix of H and set other elements to 0
209  m_fac_H.rightCols(m_ncv - from_k).setZero();
210  m_fac_H.block(from_k, 0, m_ncv - from_k, from_k).setZero();
211  for(int i = from_k; i <= to_m - 1; i++)
212  {
213  bool restart = false;
214  // If beta = 0, then the next V is not full rank
215  // We need to generate a new residual vector that is orthogonal
216  // to the current V, which we call a restart
217  if(beta < m_eps)
218  {
219  SimpleRandom<Scalar> rng(2 * i);
220  m_fac_f.noalias() = rng.random_vec(m_n);
221  // f <- f - V * V' * f, so that f is orthogonal to V
222  MapMat V(m_fac_V.data(), m_n, i); // The first i columns
223  Vector Vf = inner_product(V, m_fac_f);
224  m_fac_f.noalias() -= V * Vf;
225  // beta <- ||f||
226  beta = norm(m_fac_f);
227 
228  restart = true;
229  }
230 
231  // v <- f / ||f||
232  MapVec v(&m_fac_V(0, i), m_n); // The (i+1)-th column
233  v.noalias() = m_fac_f / beta;
234 
235  // Note that H[i+1, i] equals to the unrestarted beta
236  if(restart)
237  m_fac_H(i, i - 1) = 0.0;
238  else
239  m_fac_H(i, i - 1) = beta;
240 
241  // w <- A * v
242  m_op->perform_op(v.data(), w.data());
243  m_nmatop++;
244 
245  Hii = inner_product(v, w);
246  m_fac_H(i - 1, i) = m_fac_H(i, i - 1); // Due to symmetry
247  m_fac_H(i, i) = Hii;
248 
249  // f <- w - V * V' * w = w - H[i+1, i] * V{i} - H[i+1, i+1] * V{i+1}
250  // If restarting, we know that H[i+1, i] = 0
251  if(restart)
252  m_fac_f.noalias() = w - Hii * v;
253  else
254  m_fac_f.noalias() = w - m_fac_H(i, i - 1) * m_fac_V.col(i - 1) - Hii * v;
255 
256  beta = norm(m_fac_f);
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  MapMat V(m_fac_V.data(), m_n, i + 1); // The first (i+1) columns
261  Vector Vf = inner_product(V, m_fac_f);
262  // If not, iteratively correct the residual
263  int count = 0;
264  while(count < 5 && Vf.cwiseAbs().maxCoeff() > m_approx_0 * beta)
265  {
266  // f <- f - V * Vf
267  m_fac_f.noalias() -= V * Vf;
268  // h <- h + Vf
269  m_fac_H(i - 1, i) += Vf[i - 1];
270  m_fac_H(i, i - 1) = m_fac_H(i - 1, i);
271  m_fac_H(i, i) += Vf[i];
272  // beta <- ||f||
273  beta = norm(m_fac_f);
274 
275  Vf.noalias() = inner_product(V, m_fac_f);
276  count++;
277  }
278  }
279  }
280 
281  // Implicitly restarted Arnoldi factorization
282  void restart(int k)
283  {
284  if(k >= m_ncv)
285  return;
286 
287  TridiagQR<Scalar> decomp;
288  Matrix Q = Matrix::Identity(m_ncv, m_ncv);
289 
290  for(int i = k; i < m_ncv; i++)
291  {
292  // QR decomposition of H-mu*I, mu is the shift
293  m_fac_H.diagonal().array() -= m_ritz_val[i];
294  decomp.compute(m_fac_H);
295 
296  // Q -> Q * Qi
297  decomp.apply_YQ(Q);
298  // H -> Q'HQ
299  // Since QR = H - mu * I, we have H = QR + mu * I
300  // and therefore Q'HQ = RQ + mu * I
301  m_fac_H = decomp.matrix_RQ();
302  m_fac_H.diagonal().array() += m_ritz_val[i];
303  }
304  // V -> VQ, only need to update the first k+1 columns
305  // Q has some elements being zero
306  // The first (ncv - k + i) elements of the i-th column of Q are non-zero
307  Matrix Vs(m_n, k + 1);
308  int nnz;
309  for(int i = 0; i < k; i++)
310  {
311  nnz = m_ncv - k + i + 1;
312  MapMat V(m_fac_V.data(), m_n, nnz);
313  MapVec q(&Q(0, i), nnz);
314  Vs.col(i).noalias() = V * q;
315  }
316  Vs.col(k).noalias() = m_fac_V * Q.col(k);
317  m_fac_V.leftCols(k + 1).noalias() = Vs;
318 
319  Vector fk = m_fac_f * Q(m_ncv - 1, k - 1) + m_fac_V.col(k) * m_fac_H(k, k - 1);
320  factorize_from(k, m_ncv, fk);
321  retrieve_ritzpair();
322  }
323 
324  // Calculates the number of converged Ritz values
325  int num_converged(Scalar tol)
326  {
327  // thresh = tol * max(m_approx_0, abs(theta)), theta for ritz value
328  Array thresh = tol * m_ritz_val.head(m_nev).array().abs().max(m_approx_0);
329  Array resid = m_ritz_est.head(m_nev).array().abs() * norm(m_fac_f);
330  // Converged "wanted" ritz values
331  m_ritz_conv = (resid < thresh);
332 
333  return m_ritz_conv.cast<int>().sum();
334  }
335 
336  // Returns the adjusted nev for restarting
337  int nev_adjusted(int nconv)
338  {
339  using std::abs;
340 
341  int nev_new = m_nev;
342  for(int i = m_nev; i < m_ncv; i++)
343  if(abs(m_ritz_est[i]) < m_eps) nev_new++;
344 
345  // Adjust nev_new, according to dsaup2.f line 677~684 in ARPACK
346  nev_new += std::min(nconv, (m_ncv - nev_new) / 2);
347  if(nev_new == 1 && m_ncv >= 6)
348  nev_new = m_ncv / 2;
349  else if(nev_new == 1 && m_ncv > 2)
350  nev_new = 2;
351 
352  if(nev_new > m_ncv - 1)
353  nev_new = m_ncv - 1;
354 
355  return nev_new;
356  }
357 
358  // Retrieves and sorts ritz values and ritz vectors
359  void retrieve_ritzpair()
360  {
361  TridiagEigen<Scalar> decomp(m_fac_H);
362  Vector evals = decomp.eigenvalues();
363  Matrix evecs = decomp.eigenvectors();
364 
365  SortEigenvalue<Scalar, SelectionRule> sorting(evals.data(), evals.size());
366  std::vector<int> ind = sorting.index();
367 
368  // For BOTH_ENDS, the eigenvalues are sorted according
369  // to the LARGEST_ALGE rule, so we need to move those smallest
370  // values to the left
371  // The order would be
372  // Largest => Smallest => 2nd largest => 2nd smallest => ...
373  // We keep this order since the first k values will always be
374  // the wanted collection, no matter k is nev_updated (used in restart())
375  // or is nev (used in sort_ritzpair())
376  if(SelectionRule == BOTH_ENDS)
377  {
378  std::vector<int> ind_copy(ind);
379  for(int i = 0; i < m_ncv; i++)
380  {
381  // If i is even, pick values from the left (large values)
382  // If i is odd, pick values from the right (small values)
383  if(i % 2 == 0)
384  ind[i] = ind_copy[i / 2];
385  else
386  ind[i] = ind_copy[m_ncv - 1 - i / 2];
387  }
388  }
389 
390  // Copy the ritz values and vectors to m_ritz_val and m_ritz_vec, respectively
391  for(int i = 0; i < m_ncv; i++)
392  {
393  m_ritz_val[i] = evals[ind[i]];
394  m_ritz_est[i] = evecs(m_ncv - 1, ind[i]);
395  }
396  for(int i = 0; i < m_nev; i++)
397  {
398  m_ritz_vec.col(i) = evecs.col(ind[i]);
399  }
400  }
401 
402 protected:
403  // In generalized eigenvalue problem Ax=lambda*Bx, define the inner product to be <x, y> = x'By
404  // For regular eigenvalue problems, it is the usual inner product <x, y> = x'y
405  virtual Scalar inner_product(const Vector& x, const Vector& y) { return x.dot(y); }
406  virtual Scalar inner_product(const MapVec& x, const Vector& y) { return x.dot(y); }
407  virtual Vector inner_product(const MapMat& x, const Vector& y) { return x.transpose() * y; }
408 
409  // B-norm of a vector. For regular eigenvalue problems it is simply the L2 norm
410  virtual Scalar norm(const Vector& x) { return x.norm(); }
411 
412  // Sorts the first nev Ritz pairs in the specified order
413  // This is used to return the final results
414  virtual void sort_ritzpair(int sort_rule)
415  {
416  // First make sure that we have a valid index vector
417  SortEigenvalue<Scalar, LARGEST_ALGE> sorting(m_ritz_val.data(), m_nev);
418  std::vector<int> ind = sorting.index();
419 
420  switch(sort_rule)
421  {
422  case LARGEST_ALGE:
423  break;
424  case LARGEST_MAGN:
425  {
426  SortEigenvalue<Scalar, LARGEST_MAGN> sorting(m_ritz_val.data(), m_nev);
427  ind = sorting.index();
428  }
429  break;
430  case SMALLEST_ALGE:
431  {
432  SortEigenvalue<Scalar, SMALLEST_ALGE> sorting(m_ritz_val.data(), m_nev);
433  ind = sorting.index();
434  }
435  break;
436  case SMALLEST_MAGN:
437  {
438  SortEigenvalue<Scalar, SMALLEST_MAGN> 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  Vector new_ritz_val(m_ncv);
447  Matrix 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  SymEigsSolver(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 - 1)
492  throw std::invalid_argument("nev must satisfy 1 <= nev <= n - 1, n is the size of matrix");
493 
494  if(ncv_ <= nev_ || ncv_ > m_n)
495  throw std::invalid_argument("ncv must satisfy nev < ncv <= n, n is the size of matrix");
496  }
497 
501  virtual ~SymEigsSolver() {}
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  m_nmatop = 0;
532  m_niter = 0;
533 
534  // Set the initial vector
535  Vector v(m_n);
536  std::copy(init_resid, init_resid + m_n, v.data());
537  Scalar vnorm = norm(v);
538  if(vnorm < m_eps)
539  throw std::invalid_argument("initial residual vector cannot be zero");
540  v /= vnorm;
541 
542  Vector w(m_n);
543  m_op->perform_op(v.data(), w.data());
544  m_nmatop++;
545 
546  m_fac_H(0, 0) = inner_product(v, w);
547  m_fac_f.noalias() = w - v * m_fac_H(0, 0);
548  m_fac_V.col(0) = v;
549  }
550 
558  void init()
559  {
560  SimpleRandom<Scalar> rng(0);
561  Vector init_resid = rng.random_vec(m_n);
562  init(init_resid.data());
563  }
564 
583  int compute(int maxit = 1000, Scalar tol = 1e-10, int sort_rule = LARGEST_ALGE)
584  {
585  // The m-step Arnoldi factorization
586  factorize_from(1, m_ncv, m_fac_f);
587  retrieve_ritzpair();
588  // Restarting
589  int i, nconv = 0, nev_adj;
590  for(i = 0; i < maxit; i++)
591  {
592  nconv = num_converged(tol);
593  if(nconv >= m_nev)
594  break;
595 
596  nev_adj = nev_adjusted(nconv);
597  restart(nev_adj);
598  }
599  // Sorting results
600  sort_ritzpair(sort_rule);
601 
602  m_niter += i + 1;
603  m_info = (nconv >= m_nev) ? SUCCESSFUL : NOT_CONVERGING;
604 
605  return std::min(m_nev, nconv);
606  }
607 
612  int info() const { return m_info; }
613 
617  int num_iterations() const { return m_niter; }
618 
622  int num_operations() const { return m_nmatop; }
623 
631  Vector eigenvalues() const
632  {
633  int nconv = m_ritz_conv.cast<int>().sum();
634  Vector res(nconv);
635 
636  if(!nconv)
637  return res;
638 
639  int j = 0;
640  for(int i = 0; i < m_nev; i++)
641  {
642  if(m_ritz_conv[i])
643  {
644  res[j] = m_ritz_val[i];
645  j++;
646  }
647  }
648 
649  return res;
650  }
651 
661  Matrix eigenvectors(int nvec) const
662  {
663  int nconv = m_ritz_conv.cast<int>().sum();
664  nvec = std::min(nvec, nconv);
665  Matrix res(m_n, nvec);
666 
667  if(!nvec)
668  return res;
669 
670  Matrix ritz_vec_conv(m_ncv, nvec);
671  int j = 0;
672  for(int i = 0; i < m_nev && j < nvec; i++)
673  {
674  if(m_ritz_conv[i])
675  {
676  ritz_vec_conv.col(j) = m_ritz_vec.col(i);
677  j++;
678  }
679  }
680 
681  res.noalias() = m_fac_V * ritz_vec_conv;
682 
683  return res;
684  }
685 
689  Matrix eigenvectors() const
690  {
691  return eigenvectors(m_nev);
692  }
693 };
694 
695 
696 } // namespace Spectra
697 
698 #endif // SYM_EIGS_SOLVER_H
int num_operations() const
int compute(int maxit=1000, Scalar tol=1e-10, int sort_rule=LARGEST_ALGE)
Matrix eigenvectors(int nvec) const
Computation was successful.
Definition: CompInfo.h:20
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
Matrix eigenvectors() const
SymEigsSolver(OpType *op_, int nev_, int ncv_)