From f14be5e3316fc674a589e801d3aada6f2d306754 Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Wed, 22 Apr 2026 23:42:41 -0400 Subject: [PATCH] Implicit restart in Gamma5BlockLanczos --- .../algorithms/iterative/Gamma5BlockLanczos.h | 491 +++++++++++++++++- 1 file changed, 489 insertions(+), 2 deletions(-) diff --git a/Grid/algorithms/iterative/Gamma5BlockLanczos.h b/Grid/algorithms/iterative/Gamma5BlockLanczos.h index af5e11603..9f2a2afef 100644 --- a/Grid/algorithms/iterative/Gamma5BlockLanczos.h +++ b/Grid/algorithms/iterative/Gamma5BlockLanczos.h @@ -58,6 +58,12 @@ private: int nSteps; // number of completed steps + // Krylov-Schur implicit restart state + CMat Hmat_; // full projected matrix (2*nSteps x 2*nSteps) + bool useFullH_; // true while in Krylov-Schur extension mode + int Ncompressed_; // number of compressed column vectors kept after last KS step + CMat2 Blast_; // last normalization block from lanczosStepFull (for Ritz estimate) + // Output CVec evals_; std::vector evecs_; @@ -68,8 +74,11 @@ public: Gamma5BlockLanczos(LinearOperatorBase& op, GridBase* grid, Gamma5Func g5, RealD tol = 1e-8) - : Linop(op), Grid_(grid), applyGamma5(g5), Tolerance(tol), nSteps(0) - {} + : Linop(op), Grid_(grid), applyGamma5(g5), Tolerance(tol), nSteps(0), + useFullH_(false), Ncompressed_(0) + { + Blast_ = CMat2::Zero(); + } CVec getEvals() { return evals_; } std::vector getEvecs() { return evecs_; } @@ -137,6 +146,206 @@ public: computeRitzPairs(nSteps, Nstop); } + /** + * Restarted γ5-Block Lanczos (explicit restart). + * + * Runs operator() for Nstep steps per pass. After each pass the 2*Nstep + * Ritz pairs are sorted by the chosen RitzFilter and the top Nk are the + * "wanted" set. The seed for the next pass is the normalised equal-weight + * sum of the Nstop best Ritz vectors from the wanted set; this keeps the + * restart in the span of the most-wanted approximate eigenspace. + * + * Convergence is declared when at least Nstop of the top-Nk pairs have + * residual < tolerance. + * + * v0 : initial starting vector + * maxRestarts : maximum number of Lanczos passes + * Nstep : Lanczos steps per pass (produces 2*Nstep Ritz values) + * Nk : size of the "wanted" set to track (must be >= Nstop) + * Nstop : target converged pairs; also the number of vectors summed + * to form the restart seed + * reorthog : full γ5-reorthogonalisation within each pass + * filter : EvalNormSmall → sort by |λ| ascending + * EvalImNormSmall → sort by |Im(λ)| ascending (default) + */ + void restart(const Field& v0, int maxRestarts, int Nstep, int Nk, int Nstop, + bool reorthog = false, RitzFilter filter = EvalImNormSmall) + { + assert(Nk >= Nstop); + Field src(Grid_); + src = v0; + + for (int iter = 0; iter < maxRestarts; iter++) { + std::cout << GridLogMessage + << "Gamma5BlockLanczos: ---- restart " << iter << " ----" << std::endl; + + (*this)(src, Nstep, Nstop, reorthog); + + int nRitz = (int)residuals_.size(); + if (nRitz == 0) { + std::cout << GridLogMessage + << "Gamma5BlockLanczos: restart — no Ritz pairs, stopping." << std::endl; + return; + } + + // Sort all Ritz indices by the chosen filter criterion. + std::vector idx(nRitz); + std::iota(idx.begin(), idx.end(), 0); + switch (filter) { + case EvalNormSmall: + std::sort(idx.begin(), idx.end(), [&](int a, int b){ + return std::abs(evals_(a)) < std::abs(evals_(b)); + }); + break; + case EvalImNormSmall: + default: + std::sort(idx.begin(), idx.end(), [&](int a, int b){ + return std::abs(evals_(a).imag()) < std::abs(evals_(b).imag()); + }); + break; + } + + // Count converged pairs within the top-Nk wanted set. + int nKeep = std::min(Nk, nRitz); + int nconv = 0; + for (int i = 0; i < nKeep; i++) + if (residuals_[idx[i]] < Tolerance) nconv++; + + std::cout << GridLogMessage + << "Gamma5BlockLanczos: restart " << iter + << " Ritz = " << nRitz + << " wanted = " << nKeep + << " converged = " << nconv << " / " << Nstop << std::endl; + for (int i = 0; i < nKeep; i++) + std::cout << GridLogMessage + << " wanted[" << i << "] lambda = " << evals_(idx[i]) + << " |res| = " << residuals_[idx[i]] + << (residuals_[idx[i]] < Tolerance ? " *" : "") << std::endl; + + if (nconv >= Nstop) { + std::cout << GridLogMessage + << "Gamma5BlockLanczos: converged after " << iter + 1 + << " restart(s)." << std::endl; + reorderOutput(idx, nKeep); + return; + } + + // Build restart seed: equal-weight sum of the top Nstop Ritz vectors + // (sorted by filter criterion). Spans the best part of the wanted + // eigenspace and avoids locking onto a single approximate eigenvalue. + int nSeed = std::min(Nstop, nKeep); + src = Zero(); + for (int i = 0; i < nSeed; i++) + src += evecs_[idx[i]]; + + RealD nrm = std::sqrt(norm2(src)); + assert(nrm > 1e-14); + src *= (1.0 / nrm); + + std::cout << GridLogMessage + << "Gamma5BlockLanczos: seed = sum of top " << nSeed + << " Ritz vectors ||seed|| = " << nrm << std::endl; + } + + std::cout << GridLogMessage + << "Gamma5BlockLanczos: max restarts (" << maxRestarts + << ") reached without full convergence." << std::endl; + } + + /** + * Implicitly Restarted Block Lanczos (Krylov-Schur variant). + * + * Each outer cycle: + * 1. Runs Nmax block Lanczos steps (first cycle from v0; subsequently + * extends the Nk-step compressed state). + * 2. Applies a Krylov-Schur restart: Schur-decomposes T_{Nmax}, reorders + * by the chosen filter, and compresses the basis from 2*Nmax to Nk + * columns. The leading Nk×Nk upper-triangular Schur block S_k replaces + * the old block-tridiagonal T. + * 3. Extends from Nk/2 to Nmax steps using full γ5-projection (block + * Arnoldi with the complete Nk-vector compressed basis). + * + * Convergence is declared when Nstop pairs have residual < tolerance. + * + * v0 : initial starting vector + * maxIter : maximum restart cycles (excluding initial run) + * Nmax : Lanczos steps per cycle (builds 2*Nmax Ritz pairs) + * Nk : steps kept after compression (even, 2 ≤ Nk < Nmax, Nk ≥ Nstop) + * Nstop : target converged pairs + * reorthog : γ5-reorthogonalisation in the initial Nmax-step run + * filter : eigenvalue selection criterion (default: EvalImNormSmall) + */ + void implicitRestart(const Field& v0, int maxIter, int Nmax, int Nk, int Nstop, + bool reorthog = false, RitzFilter filter = EvalImNormSmall) + { + assert(Nk % 2 == 0 && Nk >= 2 && Nk < Nmax); + assert(Nk >= Nstop); + + // Initial full block-Lanczos run + (*this)(v0, Nmax, Nstop, reorthog); + + for (int iter = 0; iter < maxIter; iter++) { + std::cout << GridLogMessage + << "Gamma5BlockLanczos::implicitRestart ---- cycle " << iter << " ----" << std::endl; + + // Sort current Ritz pairs by filter, count converged + int nRitz = (int)residuals_.size(); + std::vector idx = sortedIdx(nRitz, filter); + int nKeep = std::min(Nk, nRitz); + int nconv = 0; + for (int i = 0; i < nKeep; i++) + if (residuals_[idx[i]] < Tolerance) nconv++; + + std::cout << GridLogMessage + << " nRitz=" << nRitz + << " nconv=" << nconv << "/" << Nstop << std::endl; + for (int i = 0; i < nKeep; i++) + std::cout << GridLogMessage + << " [" << i << "] lambda=" << evals_(idx[i]) + << " |res|=" << residuals_[idx[i]] + << (residuals_[idx[i]] < Tolerance ? " *" : "") << std::endl; + + if (nconv >= Nstop) { + reorderOutput(idx, nKeep); + std::cout << GridLogMessage + << "Gamma5BlockLanczos::implicitRestart: converged after " + << iter + 1 << " cycle(s)." << std::endl; + useFullH_ = false; + return; + } + + // Krylov-Schur: compress T_{Nmax} to S_{Nk} (upper-triangular Nk×Nk) + krylovSchurCompress(Nk, filter); + + // Extend from Nk/2 to Nmax steps via full-projection block Arnoldi + for (int step = nSteps; step < Nmax; step++) { + bool ok = lanczosStepFull(step); + if (!ok) break; + nSteps = step + 1; + if (Blast_.norm() < Tolerance) { + std::cout << GridLogMessage + << "Gamma5BlockLanczos::implicitRestart: beta < tol at full step " + << step << ", stopping extension." << std::endl; + break; + } + } + + // Ritz pairs from the full projected matrix + computeRitzPairsFull(nSteps, Nstop); + } + + // maxIter exhausted + std::cout << GridLogMessage + << "Gamma5BlockLanczos::implicitRestart: maxIter=" << maxIter + << " reached without full convergence." << std::endl; + int nRitz = (int)residuals_.size(); + if (nRitz > 0) { + std::vector idx = sortedIdx(nRitz, filter); + reorderOutput(idx, std::min(Nk, nRitz)); + } + useFullH_ = false; + } + /** * Verify the block Lanczos decomposition after operator() has run. * @@ -269,6 +478,284 @@ public: } private: + // Return a permutation of [0,n) sorted by the chosen RitzFilter criterion. + std::vector sortedIdx(int n, RitzFilter filter) + { + std::vector idx(n); + std::iota(idx.begin(), idx.end(), 0); + switch (filter) { + case EvalNormSmall: + std::sort(idx.begin(), idx.end(), [&](int a, int b){ + return std::abs(evals_(a)) < std::abs(evals_(b)); + }); + break; + case EvalImNormSmall: + default: + std::sort(idx.begin(), idx.end(), [&](int a, int b){ + return std::abs(evals_(a).imag()) < std::abs(evals_(b).imag()); + }); + break; + } + return idx; + } + + // Krylov-Schur: compress the current nSteps-step block-Lanczos basis to Nk + // steps by Schur-decomposing T_{nSteps}, reordering "wanted" eigenvalues to + // the top, and rotating the field basis accordingly. + // + // After this call: + // basis[0..Nk-1] = Nk compressed Schur vectors (ṽ_j = V_m U[:,j]) + // basis[Nk..Nk+1] = original outer-residual block Q_{m+1} (unchanged) + // G_blocks[0..Nk/2] recomputed from the new pairs + // Hmat_ = leading Nk×Nk block of the Schur form S (upper triangular) + // A/B/C blocks cleared; nSteps = Nk/2; useFullH_ = true + void krylovSchurCompress(int Nk, RitzFilter filter) + { + int m = nSteps; + int dim = 2 * m; + assert(Nk > 0 && Nk % 2 == 0 && Nk < dim); + + // Assemble T_m (block tridiagonal, dim×dim) + CMat Tm = CMat::Zero(dim, dim); + for (int k = 0; k < m; k++) { + Tm.block(2*k, 2*k, 2, 2) = A_blocks[k]; + if (k < m - 1) { + Tm.block(2*k+2, 2*k, 2, 2) = B_blocks[k]; + Tm.block(2*k, 2*k+2, 2, 2) = C_blocks[k+1]; + } + } + + // Complex Schur decomposition with wanted eigenvalues first. + // ComplexSchurDecomposition convention: A = Q† S Q, + // getMatrixQ() = U† (so U = Q† → U.adjoint() = Q) + // getMatrixS() = S (upper triangular) + // New basis: ṽ_j = V_m U[:,j] with U = getMatrixQ().adjoint() + ComplexSchurDecomposition schur(Tm, false, filter); + schur.schurReorder(Nk); + CMat U = schur.getMatrixQ().adjoint(); // rotation matrix (dim×dim unitary) + CMat S = schur.getMatrixS(); // upper-triangular Schur form + + // Build Nk compressed field vectors + keep Q_{m+1} as "next block" + std::vector new_basis; + new_basis.reserve(Nk + 2); + for (int j = 0; j < Nk; j++) { + Field col(Grid_); + col = Zero(); + for (int k = 0; k < dim; k++) + col += basis[k] * U(k, j); + new_basis.push_back(col); + } + // Outer residual block Q_{m+1} unchanged (extension starts from here) + new_basis.push_back(basis[dim]); + new_basis.push_back(basis[dim + 1]); + + basis = new_basis; // Nk + 2 field vectors + + // Recompute G_blocks for each pair of new basis vectors + G_blocks.clear(); + for (int k = 0; k <= Nk / 2; k++) + G_blocks.push_back(gramMatrix(basis[2*k], basis[2*k+1])); + + // Compressed projected matrix = leading Nk×Nk block of S + Hmat_ = S.block(0, 0, Nk, Nk); + + // Clear three-term block storage (not valid after rotation) + A_blocks.clear(); + B_blocks.clear(); + C_blocks.clear(); + + useFullH_ = true; + Ncompressed_ = Nk; + nSteps = Nk / 2; + + std::cout << GridLogMessage + << "Gamma5BlockLanczos::krylovSchurCompress: compressed to Nk=" << Nk + << " Hmat=" << Nk << "x" << Nk << std::endl; + } + + // One block Arnoldi step with full γ5-projection against all previous blocks. + // Used after krylovSchurCompress to extend from Nk/2 to Nmax steps. + // + // Computes all coupling entries T_{j,step} = G_j^{-1} Q_j† γ5 D_W Q_step + // (j=0..step) via classical Gram-Schmidt, stores them in a new column of + // Hmat_, subtracts the projections to form the residual, then LDL†-normalises + // to produce Q_{step+1}. Appends to basis and G_blocks; sets Blast_. + bool lanczosStepFull(int step) + { + const Field& q1 = basis[2*step]; + const Field& q2 = basis[2*step + 1]; + CMat2 Gk = G_blocks[step]; + + int old_dim = 2 * step; + assert(Hmat_.rows() == old_dim && Hmat_.cols() == old_dim); + + // Apply D_W to current block + Field p1(Grid_), p2(Grid_); + Linop.Op(q1, p1); + Linop.Op(q2, p2); + + // Grow Hmat_ by 2 (new column; lower rows left zero by Arnoldi projection) + int new_dim = old_dim + 2; + CMat Hmat_new = CMat::Zero(new_dim, new_dim); + Hmat_new.block(0, 0, old_dim, old_dim) = Hmat_; + + // Coupling to all previous blocks (classical GS: use original p1,p2) + Field r1(Grid_), r2(Grid_); + r1 = p1; r2 = p2; + for (int j = 0; j < step; j++) { + CMat2 Mj = g5InnerBlock(basis[2*j], basis[2*j+1], p1, p2); + CMat2 Hj = invert2x2(G_blocks[j]) * Mj; + Hmat_new.block(2*j, old_dim, 2, 2) = Hj; + r1 -= basis[2*j] * Hj(0,0) + basis[2*j+1] * Hj(1,0); + r2 -= basis[2*j] * Hj(0,1) + basis[2*j+1] * Hj(1,1); + } + + // On-diagonal block A_{step} + CMat2 Mk = g5InnerBlock(q1, q2, p1, p2); + CMat2 Ak = invert2x2(Gk) * Mk; + Hmat_new.block(old_dim, old_dim, 2, 2) = Ak; + r1 -= q1 * Ak(0,0) + q2 * Ak(1,0); + r2 -= q1 * Ak(0,1) + q2 * Ak(1,1); + + // LDL† normalisation of the residual block + CMat2 Gamma_k = gramMatrix(r1, r2); + SAEigen2 es(Gamma_k); + Eigen::Vector2d D = es.eigenvalues(); + CMat2 U2 = es.eigenvectors(); + + for (int j = 0; j < 2; j++) { + if (std::abs(D(j)) < 1e-28) { + Field rj = r1 * U2(0,j) + r2 * U2(1,j); + if (std::sqrt(norm2(rj)) < Tolerance) { + std::cout << GridLogMessage + << "Gamma5BlockLanczos: happy breakdown (full step " << step << ")" << std::endl; + } else { + std::cout << GridLogMessage + << "Gamma5BlockLanczos: serious breakdown (full step " << step + << ") — stopping." << std::endl; + return false; + } + } + } + + CMat2 Gkp1 = CMat2::Zero(); + Gkp1(0,0) = ComplexD((D(0) > 0.0) ? 1.0 : -1.0, 0.0); + Gkp1(1,1) = ComplexD((D(1) > 0.0) ? 1.0 : -1.0, 0.0); + double sqd0 = std::sqrt(std::abs(D(0))); + double sqd1 = std::sqrt(std::abs(D(1))); + + CMat2 Bkp1; + Bkp1.row(0) = U2.col(0).adjoint() * sqd0; + Bkp1.row(1) = U2.col(1).adjoint() * sqd1; + + Field qnew1 = (r1 * U2(0,0) + r2 * U2(1,0)) * (1.0 / sqd0); + Field qnew2 = (r1 * U2(0,1) + r2 * U2(1,1)) * (1.0 / sqd1); + + Hmat_ = Hmat_new; + Blast_ = Bkp1; + G_blocks.push_back(Gkp1); + basis.push_back(qnew1); + basis.push_back(qnew2); + + RealD beta = Bkp1.norm(); + std::cout << GridLogMessage + << "Gamma5BlockLanczos: full step " << step << " beta=" << beta << std::endl; + return true; + } + + // Compute Ritz pairs from the full projected matrix Hmat_ (set after + // krylovSchurCompress + lanczosStepFull calls). Ritz residual estimated + // cheaply as ||Blast_ τ_j|| where τ_j = last 2 entries of eigenvector y_j. + void computeRitzPairsFull(int m, int Nstop) + { + int dim = 2 * m; + assert(Hmat_.rows() == dim && Hmat_.cols() == dim); + + Eigen::ComplexEigenSolver ces(Hmat_); + CVec lambdas = ces.eigenvalues(); + CMat Y = ces.eigenvectors(); + + // Sort by |Im(λ)| ascending (near-real = physical modes first) + std::vector idx(dim); + std::iota(idx.begin(), idx.end(), 0); + std::sort(idx.begin(), idx.end(), [&](int a, int b){ + return std::abs(lambdas(a).imag()) < std::abs(lambdas(b).imag()); + }); + + evals_.resize(dim); + evecs_.clear(); + residuals_.clear(); + + for (int ji = 0; ji < dim; ji++) { + int j = idx[ji]; + evals_(ji) = lambdas(j); + CVec yj = Y.col(j); + + // Ritz vector: ũ_j = sum_k basis[k] * yj[k] (k = 0..dim-1) + Field uj(Grid_); + uj = Zero(); + for (int k = 0; k < dim; k++) + uj += basis[k] * yj(k); + evecs_.push_back(uj); + + // Ritz estimate: || Blast_ τ_j || (τ_j = last 2 entries of y_j) + Eigen::Vector2cd tau(yj(dim - 2), yj(dim - 1)); + RealD res = (Blast_ * tau).norm(); + residuals_.push_back(res); + + std::cout << GridLogMessage + << "Gamma5BlockLanczos (full): Ritz[" << ji << "]" + << " lambda=" << evals_(ji) + << " |res|=" << res << std::endl; + } + + if (doEvalCheck) { + Field w(Grid_); + int nCheck = std::min((int)evecs_.size(), 2 * Nstop); + for (int k = 0; k < nCheck; k++) { + Linop.Op(evecs_[k], w); + ComplexD eval_est = toStdCmplx(innerProduct(evecs_[k], w)); + w -= eval_est * evecs_[k]; + RealD res = std::sqrt(norm2(w)); + std::cout << GridLogMessage + << "Gamma5BlockLanczos: evec[" << k << "]" + << " eval_reported=" << evals_(k) + << " eval_est=" << eval_est + << " ||Av-eval*v||=" << res << std::endl; + } + } + } + + // Reorder evals_/evecs_/residuals_ so the first nKeep entries follow idx[]. + void reorderOutput(const std::vector& idx, int nKeep) + { + CVec evals_new(evals_.size()); + std::vector evecs_new; + std::vector res_new; + evecs_new.reserve(evecs_.size()); + res_new.reserve(residuals_.size()); + + for (int i = 0; i < nKeep; i++) { + evals_new(i) = evals_(idx[i]); + evecs_new.push_back(evecs_[idx[i]]); + res_new.push_back(residuals_[idx[i]]); + } + // append remaining entries not in idx[0..nKeep-1] + std::vector used(evals_.size(), false); + for (int i = 0; i < nKeep; i++) used[idx[i]] = true; + int j = nKeep; + for (int i = 0; i < (int)evals_.size(); i++) { + if (!used[i]) { + evals_new(j++) = evals_(i); + evecs_new.push_back(evecs_[i]); + res_new.push_back(residuals_[i]); + } + } + evals_ = evals_new; + evecs_ = evecs_new; + residuals_ = res_new; + } + // One Lanczos step. On success pushes Q_{step+2} and returns true. bool lanczosStep(int step, bool reorthog) {