Low-rank eigenvalue solvers for block-sparse matrix product states

1. Problem Statement

The paper addresses the numerical solution of the electronic Schrödinger equation for fermionic quantum systems, specifically focusing on finding the lowest eigenvalues (ground states) and corresponding eigenfunctions.

• Context: The problem is formulated in second quantization, where the Hamiltonian $H$ acts on a high-dimensional tensor space of occupation numbers ($F^K \simeq (\mathbb{R}^2)^{\otimes K}$).
• Constraints:
  • Fermionic Antisymmetry: Wavefunctions must be antisymmetric under particle exchange.
  • Particle Number Conservation: The system is restricted to a fixed number of particles $N$. This imposes a block-sparse structure on the tensor representations, where specific blocks in the Matrix Product State (MPS) cores are zero depending on the number of particles to the left and right of a given orbital.
• Challenges:
  • Curse of Dimensionality: The dimension of the full tensor space grows exponentially with the number of orbitals $K$.
  • Convergence of Iterative Solvers: Standard methods like Density Matrix Renormalization Group (DMRG) often lack rigorous convergence guarantees and can fail to converge due to the global structure of eigenvectors.
  • Rank Control: Iterative solvers in low-rank formats (like MPS) tend to increase tensor ranks during operator application. Controlling these ranks while maintaining accuracy is difficult.
  • Goal: Develop an iterative eigensolver that guarantees convergence in the $H^1$-norm (or equivalent energy norm), maintains the block-sparse structure, and provides near-optimal bounds on the tensor ranks required to achieve a specific accuracy.

2. Methodology

The authors propose a Preconditioned Inverse Iteration (PINVIT) scheme adapted for low-rank tensor formats, combined with adaptive rank truncation.

A. Mathematical Framework

• Representation: Eigenfunctions are approximated using Matrix Product States (MPS) (also known as Tensor Trains). The block-sparse structure arising from particle number conservation is explicitly preserved in the MPS cores.
• Preconditioning:
  • The Hamiltonian $H$ is shifted by $\gamma I$ to ensure positive definiteness ($H_\gamma = H + \gamma I$).
  • A preconditioner $S$ is constructed using exponential sum approximations of $H_\gamma^{-1/2}$.
  • $S$ is represented as a sum of Kronecker products of diagonal matrices, ensuring it preserves the block-sparse structure and has a low MPO (Matrix Product Operator) rank.
  • This leads to a symmetrically preconditioned generalized eigenvalue problem: $\langle SH_\gamma S x, y \rangle = \lambda \langle S^2 x, y \rangle$.

B. The Iterative Algorithm

The solver operates via an inner-outer loop structure:

1. Inner Iteration (Error Reduction):
   • Performs inexact preconditioned inverse iteration steps.
   • Uses an adaptive strategy to compute the residual and the step size $\omega$.
   • The step size is optimized by solving a small 2D (or $2D \times D$ for multiple eigenvalues) eigenvalue problem within the subspace spanned by the current iterate and the residual.
   • Inexactness: Residuals and iterates are computed with controlled errors. Rank truncation is applied to the residual calculation to keep computational costs manageable.
2. Outer Iteration (Complexity Reduction):
   • After a sufficient number of inner steps to reduce the error, the iterate is rank-truncated to a specific tolerance.
   • This truncation ensures the rank does not grow indefinitely and keeps the approximation near the "best possible" rank for the given accuracy.
   • The error tolerance for the outer loop is halved in each step ($\eta_{outer}^{m+1} = \eta_{outer}^m / 2$).

C. Extension to Multiple Eigenvalues

The method is generalized to approximate a subspace of $D$ eigenvectors simultaneously.

• Joint Representation: Instead of separate MPS for each eigenvector, a block tensor train format is used where the eigenvector index is treated as an additional mode.
• Subspace Iteration: The update step involves solving a reduced eigenvalue problem on the subspace spanned by the current basis and the residuals, followed by orthogonalization and joint rank truncation.

3. Key Contributions

1. Rigorous Convergence Analysis:

   • The paper provides a complete convergence analysis for the inexact PINVIT method in the context of block-sparse MPS.
   • It proves that the angle between the approximate eigenspace and the true eigenspace decreases systematically.
   • The analysis sharpens previous results (e.g., from [29]) by adapting them to the specific geometry of the preconditioned problem and the block-sparse structure.

2. Near-Optimal Rank Bounds:

   • A major theoretical contribution is the proof that the ranks of the iterates generated by the algorithm are near-optimal.
   • Specifically, the rank required to achieve an error $\epsilon$ is bounded by a constant factor times the rank of the best possible approximation of that accuracy ($r_{best}$).
   • The bound scales favorably with the number of orbitals $K$ (specifically, the number of inner iterations scales as $O(\log K)$).

3. Block-Sparse Preconditioner Construction:

   • The authors construct a low-rank preconditioner $S$ that respects the particle number conservation (block sparsity).
   • They utilize exponential sums to approximate the inverse square root of the operator, ensuring the preconditioner is a sum of low-rank MPOs, which is crucial for maintaining efficiency.

4. Adaptive Error Control:

   • The algorithm includes a posteriori error estimation based on the residual norm and the spectral gap.
   • It adaptively adjusts the truncation tolerances for residuals and iterates to ensure the theoretical convergence conditions are met without unnecessary computational overhead.

4. Numerical Results

The authors validate the method on model problems involving 1D Coulomb-like potentials with non-smooth interactions.

• Setup: Systems with $N=4$ particles and varying numbers of orbitals ($K=14$ and $K=30$).
• Performance:
  • Convergence: The method exhibits linear convergence of the eigenvalue error, consistent with theoretical predictions.
  • Rank Efficiency: The "inner-outer" scheme significantly outperforms a "naive" approach that only uses inner iterations. The outer truncation prevents rank explosion, keeping the ranks close to the theoretical best approximation.
  • Multiple Eigenvalues: The joint block MPS approach successfully computes the first 4 eigenvalues simultaneously with controlled ranks and convergence.
  • Scalability: The method remains effective as $K$ increases, demonstrating that the logarithmic dependence on $K$ for the number of inner iterations holds in practice.

5. Significance

• Theoretical Rigor: Unlike many tensor-based quantum chemistry methods (like DMRG) which are often heuristic, this work provides mathematical guarantees of convergence and rank bounds.
• Structure Preservation: By explicitly handling the block-sparse structure induced by particle number conservation, the method reduces computational costs and memory usage compared to generic tensor solvers.
• Scalability: The ability to control ranks adaptively while guaranteeing convergence makes this approach promising for large-scale quantum systems where traditional full-space diagonalization is impossible and standard tensor methods may stagnate.
• Generalizability: While focused on second-quantized fermionic systems, the underlying analysis of inexact preconditioned inverse iteration with rank truncation is applicable to other low-rank tensor problems in PDEs and data science.

In summary, this paper bridges the gap between the practical success of tensor network methods in quantum physics and rigorous numerical analysis, offering a robust, convergent, and rank-efficient solver for the electronic Schrödinger equation.