1 // Written by Anna Araslanova 2 // Modified by Yixuan Qiu 3 // License: MIT 4 5 #ifndef SPECTRA_LOBPCG_SOLVER_H 6 #define SPECTRA_LOBPCG_SOLVER_H 7 8 #include <functional> 9 #include <map> 10 11 #include <Eigen/Core> 12 #include <Eigen/SparseCore> 13 #include <Eigen/Eigenvalues> 14 #include <Eigen/SVD> 15 #include <Eigen/SparseCholesky> 16 17 #include "../SymGEigsSolver.h" 18 19 namespace Spectra { 20 21 /// 22 /// \ingroup EigenSolver 23 /// 24 25 /// *** METHOD 26 /// The class represent the LOBPCG algorithm, which was invented by Andrew Knyazev 27 /// Theoretical background of the procedure can be found in the articles below 28 /// - Knyazev, A.V., 2001. Toward the optimal preconditioned eigensolver : Locally optimal block preconditioned conjugate gradient method.SIAM journal on scientific computing, 23(2), pp.517 - 541. 29 /// - Knyazev, A.V., Argentati, M.E., Lashuk, I. and Ovtchinnikov, E.E., 2007. Block locally optimal preconditioned eigenvalue xolvers(BLOPEX) in HYPRE and PETSc.SIAM Journal on Scientific Computing, 29(5), pp.2224 - 2239. 30 /// 31 /// *** CONDITIONS OF USE 32 /// Locally Optimal Block Preconditioned Conjugate Gradient(LOBPCG) is a method for finding the M smallest eigenvalues 33 /// and eigenvectors of a large symmetric positive definite generalized eigenvalue problem 34 /// \f$Ax=\lambda Bx,\f$ 35 /// where \f$A_{NxN}\f$ is a symmetric matrix, \f$B\f$ is symmetric and positive - definite. \f$A and B\f$ are also assumed large and sparse 36 /// \f$\textit{M}\f$ should be \f$\<< textit{N}\f$ (at least \f$\textit{5M} < \textit{N} \f$) 37 /// 38 /// *** ARGUMENTS 39 /// Eigen::SparseMatrix<long double> A; // N*N - Ax = lambda*Bx, lrage and sparse 40 /// Eigen::SparseMatrix<long double> X; // N*M - initial approximations to eigenvectors (random in general case) 41 /// Spectra::LOBPCGSolver<long double> solver(A, X); 42 /// *Eigen::SparseMatrix<long double> B; // N*N - Ax = lambda*Bx, sparse, positive definite 43 /// solver.setConstraints(B); 44 /// *Eigen::SparseMatrix<long double> Y; // N*K - constraints, already found eigenvectors 45 /// solver.setB(B); 46 /// *Eigen::SparseMatrix<long double> T; // N*N - preconditioner ~ A^-1 47 /// solver.setPreconditioner(T); 48 /// 49 /// *** OUTCOMES 50 /// solver.solve(); // compute eigenpairs // void 51 /// solver.info(); // state of converjance // int 52 /// solver.residuals(); // get residuals to evaluate biases // Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> 53 /// solver.eigenvalues(); // get eigenvalues // Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> 54 /// solver.eigenvectors(); // get eigenvectors // Eigen::Matrix<Scalar, Eigen::Dynamic, 1> 55 /// 56 /// *** EXAMPLE 57 /// \code{.cpp} 58 /// #include <Spectra/contrib/SymSparseEigsSolverLOBPCG.h> 59 /// 60 /// // random A 61 /// Matrix a; 62 /// a = (Matrix::Random(10, 10).array() > 0.6).cast<long double>() * Matrix::Random(10, 10).array() * 5; 63 /// a = Matrix((a).triangularView<Eigen::Lower>()) + Matrix((a).triangularView<Eigen::Lower>()).transpose(); 64 /// for (int i = 0; i < 10; i++) 65 /// a(i, i) = i + 0.5; 66 /// std::cout << a << "\n"; 67 /// Eigen::SparseMatrix<long double> A(a.sparseView()); 68 /// // random X 69 /// Eigen::Matrix<long double, 10, 2> x; 70 /// x = Matrix::Random(10, 2).array(); 71 /// Eigen::SparseMatrix<long double> X(x.sparseView()); 72 /// // solve Ax = lambda*x 73 /// Spectra::LOBPCGSolver<long double> solver(A, X); 74 /// solver.compute(10, 1e-4); // 10 iterations, L2_tolerance = 1e-4*N 75 /// std::cout << "info\n" << solver.info() << std::endl; 76 /// std::cout << "eigenvalues\n" << solver.eigenvalues() << std::endl; 77 /// std::cout << "eigenvectors\n" << solver.eigenvectors() << std::endl; 78 /// std::cout << "residuals\n" << solver.residuals() << std::endl; 79 /// \endcode 80 /// 81 82 template <typename Scalar = long double> 83 class LOBPCGSolver 84 { 85 private: 86 typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Matrix; 87 typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector; 88 89 typedef std::complex<Scalar> Complex; 90 typedef Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic> ComplexMatrix; 91 typedef Eigen::Matrix<Complex, Eigen::Dynamic, 1> ComplexVector; 92 93 typedef Eigen::SparseMatrix<Scalar> SparseMatrix; 94 typedef Eigen::SparseMatrix<Complex> SparseComplexMatrix; 95 96 const int m_n; // dimension of matrix A 97 const int m_nev; // number of eigenvalues requested 98 SparseMatrix A, X; 99 SparseMatrix m_Y, m_B, m_preconditioner; 100 bool flag_with_constraints, flag_with_B, flag_with_preconditioner; 101 102 public: 103 SparseMatrix m_residuals; 104 Matrix m_evectors; 105 Vector m_evalues; 106 int m_info; 107 108 private: 109 // B-orthonormalize matrix M 110 int orthogonalizeInPlace(SparseMatrix& M, SparseMatrix& B, 111 SparseMatrix& true_BM, bool has_true_BM = false) 112 { 113 SparseMatrix BM; 114 115 if (has_true_BM == false) 116 { 117 if (flag_with_B) 118 { 119 BM = B * M; 120 } 121 else 122 { 123 BM = M; 124 } 125 } 126 else 127 { 128 BM = true_BM; 129 } 130 131 Eigen::SimplicialLDLT<SparseMatrix> chol_MBM(M.transpose() * BM); 132 133 if (chol_MBM.info() != Eigen::Success) 134 { 135 // LDLT decomposition fail 136 m_info = chol_MBM.info(); 137 return chol_MBM.info(); 138 } 139 140 SparseComplexMatrix Upper_MBM = chol_MBM.matrixU().template cast<Complex>(); 141 ComplexVector D_MBM_vec = chol_MBM.vectorD().template cast<Complex>(); 142 143 D_MBM_vec = D_MBM_vec.cwiseSqrt(); 144 145 for (int i = 0; i < D_MBM_vec.rows(); i++) 146 { 147 D_MBM_vec(i) = Complex(1.0, 0.0) / D_MBM_vec(i); 148 } 149 150 SparseComplexMatrix D_MBM_mat(D_MBM_vec.asDiagonal()); 151 152 SparseComplexMatrix U_inv(Upper_MBM.rows(), Upper_MBM.cols()); 153 U_inv.setIdentity(); 154 Upper_MBM.template triangularView<Eigen::Upper>().solveInPlace(U_inv); 155 156 SparseComplexMatrix right_product = U_inv * D_MBM_mat; 157 M = M * right_product.real(); 158 if (flag_with_B) 159 { 160 true_BM = B * M; 161 } 162 else 163 { 164 true_BM = M; 165 } 166 167 return Eigen::Success; 168 } 169 applyConstraintsInPlace(SparseMatrix & X,SparseMatrix & Y,SparseMatrix & B)170 void applyConstraintsInPlace(SparseMatrix& X, SparseMatrix& Y, 171 SparseMatrix& B) 172 { 173 SparseMatrix BY; 174 if (flag_with_B) 175 { 176 BY = B * Y; 177 } 178 else 179 { 180 BY = Y; 181 } 182 183 SparseMatrix YBY = Y.transpose() * BY; 184 SparseMatrix BYX = BY.transpose() * X; 185 186 SparseMatrix YBY_XYX = (Matrix(YBY).bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(Matrix(BYX))).sparseView(); 187 X = X - Y * YBY_XYX; 188 } 189 190 /* 191 return 192 'AB 193 CD' 194 */ stack_4_matricies(Matrix A,Matrix B,Matrix C,Matrix D)195 Matrix stack_4_matricies(Matrix A, Matrix B, 196 Matrix C, Matrix D) 197 { 198 Matrix result(A.rows() + C.rows(), A.cols() + B.cols()); 199 result.topLeftCorner(A.rows(), A.cols()) = A; 200 result.topRightCorner(B.rows(), B.cols()) = B; 201 result.bottomLeftCorner(C.rows(), C.cols()) = C; 202 result.bottomRightCorner(D.rows(), D.cols()) = D; 203 return result; 204 } 205 stack_9_matricies(Matrix A,Matrix B,Matrix C,Matrix D,Matrix E,Matrix F,Matrix G,Matrix H,Matrix I)206 Matrix stack_9_matricies(Matrix A, Matrix B, Matrix C, 207 Matrix D, Matrix E, Matrix F, 208 Matrix G, Matrix H, Matrix I) 209 { 210 Matrix result(A.rows() + D.rows() + G.rows(), A.cols() + B.cols() + C.cols()); 211 result.block(0, 0, A.rows(), A.cols()) = A; 212 result.block(0, A.cols(), B.rows(), B.cols()) = B; 213 result.block(0, A.cols() + B.cols(), C.rows(), C.cols()) = C; 214 result.block(A.rows(), 0, D.rows(), D.cols()) = D; 215 result.block(A.rows(), A.cols(), E.rows(), E.cols()) = E; 216 result.block(A.rows(), A.cols() + B.cols(), F.rows(), F.cols()) = F; 217 result.block(A.rows() + D.rows(), 0, G.rows(), G.cols()) = G; 218 result.block(A.rows() + D.rows(), A.cols(), H.rows(), H.cols()) = H; 219 result.block(A.rows() + D.rows(), A.cols() + B.cols(), I.rows(), I.cols()) = I; 220 221 return result; 222 } 223 sort_epairs(Vector & evalues,Matrix & evectors,SortRule SelectionRule)224 void sort_epairs(Vector& evalues, Matrix& evectors, SortRule SelectionRule) 225 { 226 std::function<bool(Scalar, Scalar)> cmp; 227 if (SelectionRule == SortRule::SmallestAlge) 228 cmp = std::less<Scalar>{}; 229 else 230 cmp = std::greater<Scalar>{}; 231 232 std::map<Scalar, Vector, decltype(cmp)> epairs(cmp); 233 for (int i = 0; i < m_evectors.cols(); ++i) 234 epairs.insert(std::make_pair(evalues(i), evectors.col(i))); 235 236 int i = 0; 237 for (auto& epair : epairs) 238 { 239 evectors.col(i) = epair.second; 240 evalues(i) = epair.first; 241 i++; 242 } 243 } 244 removeColumns(SparseMatrix & matrix,std::vector<int> & colToRemove)245 void removeColumns(SparseMatrix& matrix, std::vector<int>& colToRemove) 246 { 247 // remove columns through matrix multiplication 248 SparseMatrix new_matrix(matrix.cols(), matrix.cols() - int(colToRemove.size())); 249 int iCol = 0; 250 std::vector<Eigen::Triplet<Scalar>> tripletList; 251 tripletList.reserve(matrix.cols() - int(colToRemove.size())); 252 253 for (int iRow = 0; iRow < matrix.cols(); iRow++) 254 { 255 if (std::find(colToRemove.begin(), colToRemove.end(), iRow) == colToRemove.end()) 256 { 257 tripletList.push_back(Eigen::Triplet<Scalar>(iRow, iCol, 1)); 258 iCol++; 259 } 260 } 261 262 new_matrix.setFromTriplets(tripletList.begin(), tripletList.end()); 263 matrix = matrix * new_matrix; 264 } 265 checkConvergence_getBlocksize(SparseMatrix & m_residuals,Scalar tolerance_L2,std::vector<int> & columnsToDelete)266 int checkConvergence_getBlocksize(SparseMatrix& m_residuals, Scalar tolerance_L2, std::vector<int>& columnsToDelete) 267 { 268 // square roots from sum of squares by column 269 int BlockSize = m_nev; 270 Scalar sum, buffer; 271 272 for (int iCol = 0; iCol < m_nev; iCol++) 273 { 274 sum = 0; 275 for (int iRow = 0; iRow < m_n; iRow++) 276 { 277 buffer = m_residuals.coeff(iRow, iCol); 278 sum += buffer * buffer; 279 } 280 281 if (sqrt(sum) < tolerance_L2) 282 { 283 BlockSize--; 284 columnsToDelete.push_back(iCol); 285 } 286 } 287 return BlockSize; 288 } 289 290 public: LOBPCGSolver(const SparseMatrix & A,const SparseMatrix X)291 LOBPCGSolver(const SparseMatrix& A, const SparseMatrix X) : 292 m_n(A.rows()), 293 m_nev(X.cols()), 294 A(A), 295 X(X), 296 flag_with_constraints(false), 297 flag_with_B(false), 298 flag_with_preconditioner(false), 299 m_info(Eigen::InvalidInput) 300 { 301 if (A.rows() != X.rows() || A.rows() != A.cols()) 302 throw std::invalid_argument("Wrong size"); 303 304 //if (m_n < 5* m_nev) 305 // throw std::invalid_argument("The problem size is small compared to the block size. Use standard eigensolver"); 306 } 307 setConstraints(const SparseMatrix & Y)308 void setConstraints(const SparseMatrix& Y) 309 { 310 m_Y = Y; 311 flag_with_constraints = true; 312 } 313 setB(const SparseMatrix & B)314 void setB(const SparseMatrix& B) 315 { 316 if (B.rows() != A.rows() || B.cols() != A.cols()) 317 throw std::invalid_argument("Wrong size"); 318 m_B = B; 319 flag_with_B = true; 320 } 321 setPreconditioner(const SparseMatrix & preconditioner)322 void setPreconditioner(const SparseMatrix& preconditioner) 323 { 324 m_preconditioner = preconditioner; 325 flag_with_preconditioner = true; 326 } 327 328 void compute(int maxit = 10, Scalar tol_div_n = 1e-7) 329 { 330 Scalar tolerance_L2 = tol_div_n * m_n; 331 int BlockSize; 332 int max_iter = std::min(m_n, maxit); 333 334 SparseMatrix directions, AX, AR, BX, AD, ADD, DD, BDD, BD, XAD, RAD, DAD, XBD, RBD, BR, sparse_eVecX, sparse_eVecR, sparse_eVecD, inverse_matrix; 335 Matrix XAR, RAR, XBR, gramA, gramB, eVecX, eVecR, eVecD; 336 std::vector<int> columnsToDelete; 337 338 if (flag_with_constraints) 339 { 340 // Apply the constraints Y to X 341 applyConstraintsInPlace(X, m_Y, m_B); 342 } 343 344 // Make initial vectors orthonormal 345 // implicit BX declaration 346 if (orthogonalizeInPlace(X, m_B, BX) != Eigen::Success) 347 { 348 max_iter = 0; 349 } 350 351 AX = A * X; 352 // Solve the following NxN eigenvalue problem for all N eigenvalues and -vectors: 353 // first approximation via a dense problem 354 Eigen::EigenSolver<Matrix> eigs(Matrix(X.transpose() * AX)); 355 356 if (eigs.info() != Eigen::Success) 357 { 358 m_info = eigs.info(); 359 max_iter = 0; 360 } 361 else 362 { 363 m_evalues = eigs.eigenvalues().real(); 364 m_evectors = eigs.eigenvectors().real(); 365 sort_epairs(m_evalues, m_evectors, SortRule::SmallestAlge); 366 sparse_eVecX = m_evectors.sparseView(); 367 368 X = X * sparse_eVecX; 369 AX = AX * sparse_eVecX; 370 BX = BX * sparse_eVecX; 371 } 372 373 for (int iter_num = 0; iter_num < max_iter; iter_num++) 374 { 375 m_residuals.resize(m_n, m_nev); 376 for (int i = 0; i < m_nev; i++) 377 { 378 m_residuals.col(i) = AX.col(i) - m_evalues(i) * BX.col(i); 379 } 380 BlockSize = checkConvergence_getBlocksize(m_residuals, tolerance_L2, columnsToDelete); 381 382 if (BlockSize == 0) 383 { 384 m_info = Eigen::Success; 385 break; 386 } 387 388 // substitution of the original active mask 389 if (columnsToDelete.size() > 0) 390 { 391 removeColumns(m_residuals, columnsToDelete); 392 if (iter_num > 0) 393 { 394 removeColumns(directions, columnsToDelete); 395 removeColumns(AD, columnsToDelete); 396 removeColumns(BD, columnsToDelete); 397 } 398 columnsToDelete.clear(); // for next iteration 399 } 400 401 if (flag_with_preconditioner) 402 { 403 // Apply the preconditioner to the residuals 404 m_residuals = m_preconditioner * m_residuals; 405 } 406 407 if (flag_with_constraints) 408 { 409 // Apply the constraints Y to residuals 410 applyConstraintsInPlace(m_residuals, m_Y, m_B); 411 } 412 413 if (orthogonalizeInPlace(m_residuals, m_B, BR) != Eigen::Success) 414 { 415 break; 416 } 417 AR = A * m_residuals; 418 419 // Orthonormalize conjugate directions 420 if (iter_num > 0) 421 { 422 if (orthogonalizeInPlace(directions, m_B, BD, true) != Eigen::Success) 423 { 424 break; 425 } 426 AD = A * directions; 427 } 428 429 // Perform the Rayleigh Ritz Procedure 430 XAR = Matrix(X.transpose() * AR); 431 RAR = Matrix(m_residuals.transpose() * AR); 432 XBR = Matrix(X.transpose() * BR); 433 434 if (iter_num > 0) 435 { 436 XAD = X.transpose() * AD; 437 RAD = m_residuals.transpose() * AD; 438 DAD = directions.transpose() * AD; 439 XBD = X.transpose() * BD; 440 RBD = m_residuals.transpose() * BD; 441 442 gramA = stack_9_matricies(m_evalues.asDiagonal(), XAR, XAD, XAR.transpose(), RAR, RAD, XAD.transpose(), RAD.transpose(), DAD.transpose()); 443 gramB = stack_9_matricies(Matrix::Identity(m_nev, m_nev), XBR, XBD, XBR.transpose(), Matrix::Identity(BlockSize, BlockSize), RBD, XBD.transpose(), RBD.transpose(), Matrix::Identity(BlockSize, BlockSize)); 444 } 445 else 446 { 447 gramA = stack_4_matricies(m_evalues.asDiagonal(), XAR, XAR.transpose(), RAR); 448 gramB = stack_4_matricies(Matrix::Identity(m_nev, m_nev), XBR, XBR.transpose(), Matrix::Identity(BlockSize, BlockSize)); 449 } 450 451 // Calculate the lowest/largest m eigenpairs; Solve the generalized eigenvalue problem. 452 DenseSymMatProd<Scalar> Aop(gramA); 453 DenseCholesky<Scalar> Bop(gramB); 454 455 SymGEigsSolver<DenseSymMatProd<Scalar>, DenseCholesky<Scalar>, GEigsMode::Cholesky> 456 geigs(Aop, Bop, m_nev, (std::min)(10, int(gramA.rows()) - 1)); 457 458 geigs.init(); 459 geigs.compute(SortRule::SmallestAlge); 460 461 // Mat evecs 462 if (geigs.info() == CompInfo::Successful) 463 { 464 m_evalues = geigs.eigenvalues(); 465 m_evectors = geigs.eigenvectors(); 466 sort_epairs(m_evalues, m_evectors, SortRule::SmallestAlge); 467 } 468 else 469 { 470 // Problem With General EgenVec 471 m_info = Eigen::NoConvergence; 472 break; 473 } 474 475 // Compute Ritz vectors 476 if (iter_num > 0) 477 { 478 eVecX = m_evectors.block(0, 0, m_nev, m_nev); 479 eVecR = m_evectors.block(m_nev, 0, BlockSize, m_nev); 480 eVecD = m_evectors.block(m_nev + BlockSize, 0, BlockSize, m_nev); 481 482 sparse_eVecX = eVecX.sparseView(); 483 sparse_eVecR = eVecR.sparseView(); 484 sparse_eVecD = eVecD.sparseView(); 485 486 DD = m_residuals * sparse_eVecR; // new conjugate directions 487 ADD = AR * sparse_eVecR; 488 BDD = BR * sparse_eVecR; 489 490 DD = DD + directions * sparse_eVecD; 491 ADD = ADD + AD * sparse_eVecD; 492 BDD = BDD + BD * sparse_eVecD; 493 } 494 else 495 { 496 eVecX = m_evectors.block(0, 0, m_nev, m_nev); 497 eVecR = m_evectors.block(m_nev, 0, BlockSize, m_nev); 498 499 sparse_eVecX = eVecX.sparseView(); 500 sparse_eVecR = eVecR.sparseView(); 501 502 DD = m_residuals * sparse_eVecR; 503 ADD = AR * sparse_eVecR; 504 BDD = BR * sparse_eVecR; 505 } 506 507 X = X * sparse_eVecX + DD; 508 AX = AX * sparse_eVecX + ADD; 509 BX = BX * sparse_eVecX + BDD; 510 511 directions = DD; 512 AD = ADD; 513 BD = BDD; 514 515 } // iteration loop 516 517 // calculate last residuals 518 m_residuals.resize(m_n, m_nev); 519 for (int i = 0; i < m_nev; i++) 520 { 521 m_residuals.col(i) = AX.col(i) - m_evalues(i) * BX.col(i); 522 } 523 BlockSize = checkConvergence_getBlocksize(m_residuals, tolerance_L2, columnsToDelete); 524 525 if (BlockSize == 0) 526 { 527 m_info = Eigen::Success; 528 } 529 } // compute 530 eigenvalues()531 Vector eigenvalues() 532 { 533 return m_evalues; 534 } 535 eigenvectors()536 Matrix eigenvectors() 537 { 538 return m_evectors; 539 } 540 residuals()541 Matrix residuals() 542 { 543 return Matrix(m_residuals); 544 } 545 info()546 int info() { return m_info; } 547 }; 548 549 } // namespace Spectra 550 551 #endif // SPECTRA_LOBPCG_SOLVER_H 552