# Solving eigenvalue problems in high dimensions using contour integration and Tensor Train format

Olivier Coulaud, Luc Giraud, Martina Iannacito, Mohammad Issa

► To cite this version:

Olivier Coulaud, Luc Giraud, Martina Iannacito, Mohammad Issa. Solving eigenvalue problems in high dimensions using contour integration and Tensor Train format. RR-9586, Inria. 2025, pp.19. hal-05017425v3

HAL Id: hal-05017425

https://inria.hal.science/hal-05017425v3

Submitted on 23 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Inria

# Solving eigenvalue problems in high dimensions using contour integration and Tensor Train format

Olivier Coulaud, Luc Giraud, Martina Iannacito, Mohammad Issa

## RESEARCH REPORT

Nº 9586

April 2025

Project-Team Concace
<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

Solving eigenvalue problems in high
dimensions using contour integration and
Tensor Train format

Olivier Coulaud *, Luc Giraud*, Martina Iannacito†,
Mohammad Issa‡

Project-Team Concace

Research Report n° 9586 — April 2025 — 19 pages

**Abstract:** In high-dimensional settings, solving eigenvalue problems is hindered by the curse of dimensionality, particularly when only a subset of eigenpairs within a prescribed spectral interval is sought. In this work, we investigate an adaptation of the FEAST algorithm, originally developed for symmetric eigenproblems based on contour integration, to computations where both operators and vectors are represented in the Tensor Train (TT) format. This representation drastically reduces memory and computational demands. We introduce an adaptive scheme for determining the projection subspace dimension by incorporating a rank-revealing Modified Gram-Schmidt procedure with pivoting tailored to TT-vectors. A perturbation-based analysis provides explicit bounds on the attainable residual accuracy, from which we derive a robust stopping criterion for the proposed TT-FEAST algorithm. Moreover, we design a continuation strategy that gradually refines convergence and rounding tolerances to effectively control memory growth during iterations. To demonstrate the effectiveness of TT-FEAST as a viable alternative to existing high-dimensional eigensolvers when a few eigenvalues are required, we present numerical experiments on problems up to twelve dimensions, including the Laplacian and a vibrational Hamiltonian operator.

**Key-words:** eigenvalue problem, contour integral, tensor train, attainable accuracy

<footnote>* Inria, Concace, France</footnote>

<footnote>† Dipartimento di Matematica and (AM)2, Alma Mater Studiorum Università di Bologna</footnote>

<footnote>‡ Université Paris-Saclay, CEA, List, F-91120, Palaiseau, France</footnote>
# Résolution des problèmes de valeurs propres en haute dimension à l'aide de l'intégration des contours et du format Tensor Train

**Résumé** : En grande dimension, la résolution d'un problème de valeurs propres rencontre des difficultés liées à la malédiction de la dimension, en particulier lorsqu'il s'agit de ne rechercher qu'un sous-ensemble de paires propres dans un intervalle donné. Dans ce travail, nous étudions une adaptation de l'algorithme FEAST, initialement développé pour les problèmes propres symétriques basés sur l'intégration de contour, à des calculs dans lesquels les opérateurs et les vecteurs sont représentés au format Tensor Train (TT). Cette représentation permet de réduire considérablement les besoins en mémoire et en calcul. Nous proposons un schéma adaptatif pour déterminer la dimension du sous-espace de projection, en incorporant une procédure de Gram-Schmidt modifiée qui trouve le rang, avec un pivotage adapté aux vecteurs TT. Une analyse basée sur les perturbations fournit une borne explicite sur la précision atteignable des résidus, à partir desquelles nous dérivons un critère d'arrêt robuste pour l'algorithme TT-FEAST. Nous concevons également une stratégie de continuation qui affine progressivement les tolérances de convergence et de compression afin de contrôler efficacement la croissance de la mémoire pendant les itérations. Pour démontrer l'efficacité de TT-FEAST en tant qu'alternative viable aux solveurs de valeurs propres en grande dimension existants lorsque quelques valeurs propres sont requises, nous présentons des expériences numériques sur des problèmes allant jusqu'à douze dimensions, incluant le Laplacien et un opérateur hamiltonien vibrationnel.

**Mots-clés** : problème aux valeurs propres, intégrale de contour, Train de Tenseur, précision atteignable
# Contents

<table><tr><td>1 Introduction</td><td>3</td></tr><tr><td>2 The FEAST algorithm in a Tensor Train computation framework</td><td>4</td></tr><tr><td>2.1 FEAST in matrix computation</td><td>4</td></tr><tr><td>2.2 FEAST in Tensor Train format</td><td>5</td></tr><tr><td>3 New challenges due to TT calculation</td><td>8</td></tr><tr><td>3.1 Determining the ε-rank and range of the projector</td><td>8</td></tr><tr><td>3.2 Some insight from perturbation theory analysis and stopping criterion</td><td>9</td></tr><tr><td>3.3 Memory growth management policy using a continuation procedure</td><td>11</td></tr><tr><td>4 Numerical results</td><td>11</td></tr><tr><td>4.1 Numerical test cases</td><td>11</td></tr><tr><td>4.2 Impact of the number of Gauss points for the contour integration</td><td>12</td></tr><tr><td>4.3 Attainable accuracy versus rounding threshold</td><td>12</td></tr><tr><td>4.4 Auto-tuning of the projection space dimension</td><td>14</td></tr><tr><td>4.5 Performance and advantages of ε<sub>round</sub> the continuation procedure</td><td>16</td></tr><tr><td>5 Concluding remarks</td><td>18</td></tr></table>

# 1 Introduction

We consider the solution of eigenvalue problems of the form

$$Ax = \lambda x, \qquad (1)$$

where **A** is a symmetric linear operator and **x** is a vector defined on a finite dimensional space over $\mathbb{R}$. We seek eigenpairs whose eigenvalues lie within a specific region of the spectrum. In the classical framework, **A** can be represented as a real or complex matrix.

In many scientific applications, particularly in Schrödinger-like eigenvalue problems related to vibrational spectra in quantum chemistry [3, 14], the *d*-dimensional space where eigenvectors are sought is the Cartesian product of one-dimensional spaces. For a molecule with N atoms, the dimension *d* is 3N-6, which leads to a high-dimensional space, even for small molecules. Consequently, a classical approach results in matrix and vector sizes that are prohibitively large for usual matrix computations. Indeed, a *d*-dimensional space discretized by *n* points in each direction results in a vector of size *n*<sup>*d*</sup>, which quickly leads to prohibitive storage and computational costs. Various algorithms have been developed for calculating eigenvalues at the extremes of the spectrum in both matrix and tensor formats [2, 3, 5, 12, 14, 20].

In order to mitigate the curse of dimensionality, various approaches have been proposed for the computation of eigenpairs using the tensor train (TT) format. In [5], an approach based on alternating minimization of the block Rayleigh quotient is presented for calculating minimal eigenpairs of Hermitian TT-matrices. This method performs well when computing a small number of eigenpairs, though efficiency decreases when many eigenvalues are required due to rank growth in the block TT structure. To address rank growth, [19] proposes an alternating projection method on the manifold of fixed TT-ranks based on Riemannian optimization, with subspace acceleration using the Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) method. In [20], two concepts are introduced for calculating vibrational spectra of molecules: a manifold preconditioner for a tensor version of LOBPCG, and a tensor version of Manifold Projected Simultaneous Inverse Iteration (MP SII). Both approaches assume the solution belongs to a specific low-parametric manifold.

In this work, we investigate how to adapt the FEAST algorithm to the TT-format, with the aim of providing an additional option to compute eigenpairs associated with eigenvalues in a specific region of the spectrum in high-dimensional settings. The FEAST algorithm [18] is a well-established contour integration method for computing eigenvalues of Hermitian matrices within a given closed domain of the complex plane. It iteratively refines the eigenspace by solving a Rayleigh-Ritz problem. The subspace
projector that approximates the eigenspace corresponding to the desired eigenpairs is constructed using contour integration of the resolvent, an approach inspired by density matrix representations in quantum mechanics. This method provides robust and efficient convergence, particularly when computing multiple eigenvalues simultaneously and lends itself naturally to parallel implementation [10].

Our objective is to show how the FEAST solver for the eigenvalue problem $Ax = \lambda x$ can be extended to the situation, where both $A$ and the eigenvectors $x$ are represented in TT-format. The paper is organized as follows. In Section 2, we present the FEAST method in the matrix case, we recall the key concepts of the TT framework, and we introduce the extension of FEAST to the TT-format. In Section 3, we handle the main challenges due to the TT computation. We provide an estimation procedure of the dimension of the search space, and derive a bound on the attainable accuracy as well as a robust stopping criterion supported by some perturbation theory insights. Additionally, we propose a robust algorithm that allows control over the tensor ranks, and consequently the memory consumption, particularly during the first iterations of the algorithm. Finally, in Section 4, we present numerical results for the Laplacian, and more challenging problems involving the computation of part of the infrared spectrum in 6 and 12 dimensions.

## 2 The FEAST algorithm in a Tensor Train computation framework

For clarity in notation, we use lowercase letters from the Greek alphabet, such as $\alpha$, to denote scalars, bold lowercase letters from the Latin alphabet, such as $x$, to denote vectors, and bold uppercase letters from the Latin alphabet, such as $A$, to denote matrices. The uppercase calligraphic letters, such as $\mathcal{G}$, will be used to denote a tensor, either in dense format or in TT-vector format, while the bold uppercase calligraphic letters, such as $A$, will denote a TT-matrix.

However, in the context of TT decomposition, we adopt the standard convention of using Latin letters for scalar quantities related to the decomposition, such as $r_i$ for the TT-ranks and $n_i$ for the mode sizes. This choice corresponds to classical notation in the literature, ensuring clarity when discussing the decomposition structure.

### 2.1 FEAST in matrix computation

In this section, we briefly introduce the key components of the FEAST algorithm [18, 22], which can be viewed as a subspace iteration eigensolver preconditioned by an approximate spectral projection acting as a filter. We describe the algorithm in its classical matrix computation framework for solving the symmetric eigenvalue problem:

$$Ax = \lambda x, \qquad (2)$$

where $A = A^T \in \mathbb{R}^{n \times n}$ and $(\lambda, x)$ is an eigenpair with $\lambda \in \mathbb{R}$ and $x \in \mathbb{R}^n$.

To compute a subset of $m$ eigenpairs for eigenvalues within a prescribed close region $\Gamma$, the Cauchy theorem [1] and complex contour integration are used to compute the projection onto the desired $m$-invariant subspace. Let $(x_i)_{i=1,...,m}$ denote the eigenvectors that form an orthogonal eigenbasis of the desired invariant subspace. The orthogonal projection of a vector $y$ onto $\text{span}((x_i)_{i=1}^m)$ can be written in terms of the contour integral [21] as

$$\mathbf{P}_{\Gamma} y = (\mathbf{X}_m \mathbf{X}_m^T) y = \frac{1}{2\pi i} \oint_{\Gamma} (z\mathbf{I} - \mathbf{A})^{-1} y dz, \qquad (3)$$

where $\mathbf{X}_m = [\mathbf{x}_1, \dots, \mathbf{x}_m] \in \mathbb{R}^{n \times m}$.

The projection defined by Equation (3) provides a way to compute a basis $\mathbf{W} = [\mathbf{w}_1, \dots, \mathbf{w}_m]$ of the desired invariant subspace using $m$ linearly independent vectors $(y_i)_{i=1,...,m}$, such that $\mathbf{w}_i = \mathbf{P}_{\Gamma} y_i$. For practical computational reasons, the contour $\Gamma$ is usually a circle centered on $\gamma$ that is the middle of the interval on which the eigenvalues are sought. Because the integral cannot be computed exactly, it is approximated using quadrature rules. We exploit the symmetry of the contour so that the set of integration points is chosen in the strict upper half of the complex plane, as depicted in Figure 1. The numerical ingredients involved in the quadrature rule are the $p$ weights $\omega_k$ and associated integration points $c_k$.
<figure><img src="image_4.png" /><figcaption>Figure 1: Part of the complex plane with eigenvalues on the real axis, where the dashed blue half-circle curve with center $\gamma$ surrounds the desired eigenvalues. This curve is used to approximate the spectral projector via contour integration through the $p$ quadrature points $c_i$.</figcaption></figure>

This choice has two main advantages:

1. It reduces the computation cost for the projection calculation, as it can be performed by solving only $p$ linear systems:

$$ \frac{1}{2\pi i} \oint_{\Gamma} (z\mathbf{I} - \mathbf{A})^{-1} \mathbf{y} \, dz \approx 2 \sum_{k=1}^{p} \mathrm{Re} \, (\omega_k (c_k \mathbf{I} - \mathbf{A})^{-1} \mathbf{y}), $$

where $\mathrm{Re}(\cdot)$ denotes the function that extracts the real part of a complex number. Since the purpose of the projection is to compute vectors in $\mathrm{span}((\mathbf{x}_i)_{i=1}^m)$, the constant 2 that appears on the right-hand side of the above equality will be discarded in the sequel. Furthermore, because the projection is not exact, it will be used as a preconditioner for a subspace iteration scheme.

2. The matrices that appear in the projection calculation cannot be singular, as the integration points are complex and thus cannot coincide with an eigenvalue.

With these components, the FEAST algorithm iteratively computes a sequence of bases $\mathbf{W}$, which are used to compute approximations of the $m$ desired eigenpairs via a Rayleigh-Ritz procedure. The algorithm defines two reduced matrices: $\mathbf{A}_{\mathbf{W}} = \mathbf{W}^T \mathbf{A} \mathbf{W}$ and $\mathbf{B}_{\mathbf{W}} = \mathbf{W}^T \mathbf{W}$, and then solves the reduced generalized eigenvalue problem:

$$ \mathbf{A}_{\mathbf{W}} \mathbf{\Phi} = \mathbf{B}_{\mathbf{W}} \mathbf{\Phi} \mathbf{\Theta}, \quad (4) $$

where $\mathbf{\Phi} = [\Phi_1, \dots, \Phi_m]$ and $\mathbf{\Theta} = \mathrm{diag}(\theta_1, \dots, \theta_m)$. The computed Ritz-pairs $(\theta_i, \mathbf{W}\mathbf{\Phi}_i)$ are the approximations of the sought eigenpairs $(\lambda_i, \mathbf{x}_i)$.

The main steps of the FEAST algorithm are summarized in the Algorithm 1, where the selected stopping criterion at iteration $j$ is based on a scaled residual, that gives an upper bound on the relative distance to the eigenvalues [21, Corollary 3.3]. For a threshold value $\varepsilon$, the algorithm terminates when

$$ \max_{i=1,\dots,m} \frac{\|\mathbf{A}\mathbf{x}_i^{(j)} - \lambda_i^{(j)}\mathbf{x}_i\|_2}{|\lambda_i^{(j)}|} \leq \varepsilon. $$

In practice, a key factor influencing the effectiveness of the FEAST algorithm is the determination of the number of eigenvalues $m$ within the selected contour $\Gamma$. For symmetric matrices, an accurate and robust stochastic estimator [4] is available, which estimates the number of eigenvalues inside the contour using a set of vectors whose entries are drawn from a normal random distribution. However, the absence of such an effective estimator in the context of low-rank tensors presents a significant challenge, which will be addressed in the design of the algorithm described in the next section.

## 2.2 FEAST in Tensor Train format

**Tensor Train Notations** Every tensor $\mathcal{G} \in \mathbb{R}^{n_1 \times \cdots \times n_d}$ with order $d$ and dimension $n_k$ for each mode $k$, where $k \in \{1, \dots, d\}$ admits a Tensor Train (TT) decomposition (see [17]), which is expressed for each entry of the tensor as:

$$ \mathcal{G}(i_1, \dots, i_d) = \mathbf{G}_1(i_1) \mathbf{G}_2(i_2) \cdots \mathbf{G}_d(i_d), \quad (5) $$
<caption>**Algorithm 1: FEAST algorithm**</caption>

**Input:** A matrix **A**, the interval $[\lambda_l, \lambda_u]$ where the *m* eigenvalues are sought, and *m* random vectors that form the columns of **Y**.
**Output:** Eigenpairs $(\mathbf{\Lambda}, \mathbf{X})$ of **A** inside $[\lambda_l, \lambda_u]$.
**while** not converged **do**
1 &nbsp;&nbsp;&nbsp;&nbsp;Approximate projection onto the invariant subspace $\mathbf{W} = \mathbf{P}_{\Gamma}\mathbf{Y}$ using Equation (3)
2 &nbsp;&nbsp;&nbsp;&nbsp;Form the reduced matrices $\mathbf{A}_{\mathbf{W}} = \mathbf{W}^T\mathbf{A}\mathbf{W}$ and $\mathbf{B}_{\mathbf{W}} = \mathbf{W}^T\mathbf{W}$ in $\mathbb{R}^{m \times m}$
3 &nbsp;&nbsp;&nbsp;&nbsp;Solve the reduced eigenvalue problem $\mathbf{A}_{\mathbf{W}}\mathbf{\Phi} = \mathbf{B}_{\mathbf{W}}\mathbf{\Phi}\mathbf{\Theta}$
4 &nbsp;&nbsp;&nbsp;&nbsp;Select the eigenpairs $(\mathbf{\Lambda}, \mathbf{X}) = \{(\theta_i, \mathbf{W}\mathbf{\Phi}_i)\}$ for $\theta_i \in [\lambda_l, \lambda_u]$
5 &nbsp;&nbsp;&nbsp;&nbsp;Set $\mathbf{Y} = \mathbf{X}$
**end**
**return** $\mathbf{\Lambda}$, and $\mathbf{X}$

where $\mathbf{G}_k(i_k)$ is an $r_{k-1} \times r_k$ matrix for all $i_k \in [1, n_k]$, so that $\mathbf{G}_k$ is a third-order tensor of size $r_{k-1} \times n_k \times r_k$, called the $k^{th}$ TT-core of the decomposition. The vector of ranks, $\mathbf{r} = (r_0, ..., r_d)$, is the TT-rank of the tensor with the conditions $r_0 = r_d = 1$. The product between each pair of TT-cores is the classical matrix-matrix product. Note that while $\mathbf{G}_1(i_1)$ and $\mathbf{G}_d(i_d)$ are vectors, we simplify the notation by treating them as matrices. Such a decomposition, as defined by Equation (5), is referred to as a TT-vector in the following. The memory footprint required to store a tensor in TT-format is $O(\text{ndr}^2)$, where $n = \max_{k=1,...,d}(n_k)$ and $r = \max_{k=1,...,d}(r_k)$. This is significantly smaller than the $O(n^d)$ cost of a dense representation, as the TT storage scales linearly with the tensor order.

A multilinear operator acting on tensor space can be represented by a tensor $\mathbf{A} \in \mathbb{R}^{(n_1 \times \cdots \times n_d) \times (m_1 \times \cdots \times m_d)}$. The format defined in Equation (5) extends naturally to the so-called TT-matrix, which represents such a tensor. The TT-matrix format of $\mathbf{A}$ is written as:

$$ \mathbf{A}(i_1, \cdots, i_d, j_1, \cdots, j_d) = \mathbf{A}_1(i_1, j_1)\mathbf{A}_2(i_2, j_2) \cdots \mathbf{A}_d(i_d, j_d), $$

where $\mathbf{A}_k(i_k, j_k)$ is a $r_{k-1} \times r_k$ matrix for $k = 1, \dots, d$, with $r_0 = r_d = 1$, so that $\mathbf{A}_k$ is a fourth-order tensor of size $r_{k-1} \times n_k \times n_k \times r_k$, called the $k^{th}$ TT-core of the decomposition.

Given a TT-matrix $\mathbf{A} \in \mathbb{R}^{(n_1 \times \cdots \times n_d) \times (m_1 \times \cdots \times m_d)}$ with TT-rank $\mathbf{r}$ and a TT-vector $\mathcal{X} \in \mathbb{R}^{m_1 \times \cdots \times m_d}$ with TT-rank $\mathbf{s}$, their product is a contraction over the modes, which defines a TT-vector $\mathcal{Y} = \mathbf{A}\mathcal{X} \in \mathbb{R}^{n_1 \times \cdots \times n_d}$. The $k^{th}$ core of $\mathcal{Y}$ is a third-order tensor of size at most $r_{k-1}s_{k-1} \times n_k \times r_k s_k$.

Controlling the TT-rank growth in arithmetic operations within the TT-format, such as the addition of two TT-vectors or the TT-matrix/TT-vector product, is crucial. The resulting rank of such operations is not necessarily optimal, therefore, re-compressing the resulting tensor to reduce its rank is important (see Algorithm 2, TT-rounding, in [17]). In the following, the **round** method refers to this rounding algorithm, applied with a given relative accuracy denoted by $\varepsilon_{\text{round}}$. An important property of this rounding procedure, that we use in our analysis, is the relative distance between the input and output TT-vector. More precisely if $\bar{\mathcal{X}} = \text{round}(\mathcal{X}, \varepsilon_{\text{round}})$ then the numerical error in Frobenius norm of the **round** procedure ensures that

$$ \|\bar{\mathcal{X}} - \mathcal{X}\|_{\mathrm{F}} = \delta_{\mathcal{X}} \|\mathcal{X}\|_{\mathrm{F}} \quad \text{with} \quad \delta_{\mathcal{X}} \le \varepsilon_{\text{round}}. $$

For a more detailed description of TT-decomposition and arithmetic operations on such formats, we refer to [16, 17].

We should mention the work described in [5] for the solution of eigenproblems that also relies on set of *m* TT-vectors, where the authors introduce the block TT-vector formalism. The *m* TT-vectors are combined into one TT-vector of dimension *d* + 1, where the additional dimension corresponds to the number of vectors, *m*. In the world of matrices, this is equivalent to storing a set of vectors as a matrix. The main drawback of this approach is that it merges all TT-vectors into a single entity, thereby enforcing the same low-rank structure on each of them. In other words, the block storage format forces all TT-vectors to share the same core, consequently rank, which may be higher than necessary for individual vectors. This can lead to larger memory usage and increased computational costs.

We thus prefer to introduce the notation $[\mathcal{X}_1, \dots, \mathcal{X}_m]$ for an ordered set of *m* TT-vectors, allowing each TT-vector to maintain its own rank.
**TT-FEAST algorithm** The eigenvalue problem, defined on a vector space of real tensors of order $d$, is expressed as:

$$ \mathcal{A}\mathcal{X} = \lambda\mathcal{X}, \qquad (6) $$

where $\mathcal{A} \in \mathbb{R}^{(n_1 \times \dots \times n_d) \times (n_1 \times \dots \times n_d)}$ is a symmetric TT-matrix and $\mathcal{X} \in \mathbb{R}^{n_1 \times \dots \times n_d}$ is a TT-vector. In practical computation, for a large number of modes, $d$, or large size in each mode, $n_i$, storing the full tensor and tensor matrix in memory becomes infeasible. Therefore, the numerical solution technique operates directly on a low-rank tensor approximation, similar in concept to Algorithm 1, which works numerically with floating-point numbers, even though the problem is mathematically formulated in real numbers.

The sketch of the TT-FEAST algorithm is shown in Algorithm 2, where an additional step (Line 5) is introduced to compute a well-conditioned basis of the search space, for reasons that will be clarified in Section 3.1. Conceptually, in the TT-format, the structure of the algorithm is a straightforward extension

<table><caption>**Algorithm 2: TT-FEAST algorithm**</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>Input:</td><td>A TT-matrix A, the interval [λ<sub>ℓ</sub>, λ<sub>u</sub>] where the m eigenvalues are sought, and m random vectors Y = [Y<sub>1</sub>, ..., Y<sub>m</sub>]</td></tr><tr><td>Output:</td><td>Eigenpairs (Λ, X) of A inside [λ<sub>ℓ</sub>, λ<sub>u</sub>] of center γ</td></tr><tr><td>while not converged do</td><td></td></tr><tr><td>1</td><td>Approximate projection onto the invariant subspace</td></tr><tr><td>2</td><td>for i ← 0 to m do</td></tr><tr><td>3</td><td>for k ← 0 to p do</td></tr><tr><td>4</td><td>Solve (c<sub>k</sub>I − A)Z = Y<sub>i</sub> using AMEN with convergence threshold ε<sub>Amen</sub> = ε<sub>round</sub></td></tr><tr><td>5</td><td>W<sub>i</sub> = round(W<sub>i</sub> + Re(w<sub>k</sub>(c<sub>k</sub> − γ)Z), ε<sub>round</sub>)</td></tr><tr><td>6</td><td>end</td></tr><tr><td>7</td><td>end</td></tr><tr><td>8</td><td>Orthogonalize [W<sub>1</sub>, ..., W<sub>m</sub>] using a Gram Schmidt variant to compute Q = [Q<sub>1</sub>, ..., Q<sub>m</sub>]</td></tr><tr><td>9</td><td>Form the reduced matrices A<sub>Q</sub> = Q<sup>T</sup>AQ and B<sub>Q</sub> = Q<sup>T</sup>Q</td></tr><tr><td>10</td><td>Solve the reduced generalized eigenvalue problem A<sub>Q</sub>Φ = B<sub>Q</sub>ΦΘ</td></tr><tr><td>11</td><td>Select the eigenpairs (Λ, X) = {(θ<sub>i</sub>, QΦ<sub>i</sub>)} for θ<sub>i</sub> ∈ [λ<sub>ℓ</sub>, λ<sub>u</sub>]</td></tr><tr><td>12</td><td>Set Y = X</td></tr><tr><td>13</td><td>end</td></tr><tr><td>14</td><td>return Λ, and X</td></tr></tbody></table>

of the one in classical matrix computations but novel difficulties arise. To solve the high-dimensional linear system in Line 3, we use the Alternating Minimal ENergy (AMEN) method [6]. It succeeds in "bridging the optimization in tensor formats and the classical iterative methods of the numerical linear algebra" [6]. In practice, AMEN alternates a local optimization phase on a single TT-core by forming and solving a matrix linear system, and an expansion phase for the considered TT-core, recycling the residual information.

The key difference is that, whenever a TT-vector is computed or updated, a rounding step is performed to control and prevent an unaffordable growth of the TT-ranks. The rounding operations performed when forming the matrix involved in the reduced generalized problem have an impact on the attainable accuracy, as discussed in Section 3.2. Another challenge lies in the lack of a reliable estimator for the number of eigenvalues within a given contour. The estimator originally proposed in the matrix computation framework [4] is impractical when the vectors are replaced by TT-vectors. Notice that this specific estimator can be computed using dense tensors generated from inputs having identical characteristics. Nevertheless, this approach becomes unfeasible in high-dimensional cases because storing dense tensors require excessive memory. To the best of our knowledge, it is not currently possible to generate such tensors directly in a low-rank format, and consequently, neither the estimator in low-rank format. Thus, an alternative estimation procedure must be designed. Section 3.1 explains how to estimate the correct number of eigenvalues during each iteration. Moreover, we propose a continuation procedure, described in Section 3.3, to limit the memory overhead caused by the use of random TT-vectors to span the initial search space.
# 3 New challenges due to TT calculation

## 3.1 Determining the ε-rank and range of the projector

As previously discussed, a crucial parameter in the FEAST algorithm is the dimension of the search subspace defined by the projection. This dimension should be sufficiently large to fully capture the targeted invariant subspace, yet not so large as to incur unnecessary computations.

In matrix computations, the mathematical rank of the projector is exactly $m$, the number of eigenvalues inside the contour $\Gamma$, but the numerical rank of its approximation [8] obtained through the contour integral quadrature is typically slightly larger. When no estimator is available, the governing idea to determine $m$ is to start with an initial guess value $m_0$, generate $m_0$ random vectors $\mathbf{Y} = [\mathbf{y}_1, \cdots, \mathbf{y}_{m_0}]$, and compute the $\varepsilon$-numerical rank [8] of $\mathbf{W} = [\mathbf{w}_1, \cdots, \mathbf{w}_{m_0}]$, denoted $r_{\varepsilon}$, where $\mathbf{w}_i = \mathbf{P}_{\Gamma}(\mathbf{y}_i)$. If $r_{\varepsilon} = m_0$, it is very likely that the desired space span{$\mathbf{w}_1, \cdots, \mathbf{w}_{m_0}$} is a subspace of the sought invariant subspace, meaning that $m_0$ should be increased and new directions introduced. Conversely if $r_{\varepsilon} \ll m_0$, $m_0$ is likely to be decreased to $r_{\varepsilon}$ to best focus on the invariant space. For a given set of vectors gathered in a matrix $\mathbf{W} \in \mathbb{R}^{n \times m_0}$, its $\varepsilon$-numerical rank [8] is defined as:

$$r_{\varepsilon} = \text{rank}(\mathbf{W}, \varepsilon_{\text{rank}}) = \min_{\|\mathbf{W} - \mathbf{Z}\|_2 \le \varepsilon_{\text{rank}}} \text{rank}(\mathbf{Z}).$$

The $\varepsilon$-rank is nicely characterized in terms of the singular values of $\mathbf{W}$, as singular values indicate how close a given matrix is to a matrix of lower rank. If $(\sigma_i)_{i=1,...,m_0}$ are the singular values of $\mathbf{W}$, we have

$$\sigma_1 \ge \cdots \ge \sigma_{r_{\varepsilon}} \ge \varepsilon_{\text{rank}} \ge \sigma_{r_{\varepsilon}+1} \ge \cdots \ge \sigma_{m_0}.$$

This effectively provides a way to determine the numerical rank of the matrix $\mathbf{W}$.

Computing the numerical rank can be implemented via a rank-revealing QR factorization. If Modified Gram-Schmidt with column pivoting (MGSp) is used to compute a QR factorization of $\mathbf{W}$, it can be shown [13, p. 512-513] that the diagonal entries of the $R$ factor, which are positive, satisfy $\sigma_i \le r_{i,i}$. In practice, $r_{i,i}$ is an approximation of $\sigma_i$, so that

$$\frac{\sigma_i}{\sigma_1} \approx \frac{r_{i,i}}{r_{1,1}},$$

is an estimate of the inverse of the condition number of the first $i$ vectors processed by MGSp. Consequently MGSp can be used to implement a rank revealing QR factorization and to compute an estimate of the condition number of the set of input vectors. In finite precision calculation, to detect the rank deficiency, one usually sets $\varepsilon = u\|\mathbf{Z}\|_2$, where $u$ is the unit round-off of the working arithmetic. As stated in [8, p. 276], “However, if the general level of relative error in the data is larger than $u$, then the $\varepsilon$ should be correspondingly bigger, e.g., $\varepsilon = 10^{-2}\|\mathbf{Z}\|_{\infty}$ if the entries in $\mathbf{Z}$ are correct to two digits.” In the TT-context, because the $\delta$-rounding guarantees a relative error bound in terms of Frobenius norm we will rather select

$$\varepsilon_{\text{rank}} = \delta \|\mathbf{Z}\|_{\text{F}} \approx \delta r_{1,1}.$$

To translate these ideas from matrix computation to the TT-context, we propose to compute an estimate of the numerical rank and range of the set of TT-vectors {$\mathcal{W}_1, \cdots, \mathcal{W}_{m_0}$} using MGSp. This enables us to define, for a set of TT-vectors, what could be referred to as estimates of their singular values and their numerical rank. In addition, it provides an orthonormal basis {$\mathcal{Q}_1, \cdots, \mathcal{Q}_{r_{\varepsilon}}$} of the numerical range of the space spanned by {$\mathcal{W}_1, \cdots, \mathcal{W}_{m_0}$}.

In the TT-FEAST algorithm, the set of well-conditioned TT-vectors {$\mathcal{Q}_1, \cdots, \mathcal{Q}_{r_{\varepsilon}}$} is used to build the matrices $\mathbf{A}_{\mathcal{Q}}$ and $\mathbf{B}_{\mathcal{Q}}$, which are involved in the reduced generalized eigenvalue problem. Furthermore, we vary the numerical rank threshold $\varepsilon_{\text{rank}}$ to make it less stringent during the first iterations to maintain a sufficiently large search space. This helps remove the initial randomness and allows the convergence. Once the projector convergences, the threshold is increased to discard the directions that are not in, or close to, the searched invariant subspace (i.e., eigenvectors associated with eigenvalues outside the contour but close to its boundary).

We display the Modified Gram Schmidt procedure with TT-vector pivoting in Algorithm 3. The ratio $r_{m,m}/r_{1,1}$ serves as an estimate of the condition number of the $\mathcal{W}_j$ used to construct the well-conditioned (more precisely, almost orthonormal, due to rounding) set of TT-vectors {$\mathcal{Q}_1, \cdots, \mathcal{Q}_m$}.
**Algorithm 3:** $\mathcal{Q}, r_{m,m} = \text{MGSp}(\mathcal{W}, \varepsilon_{\text{rank}})$ - Right looking Modified Gram-Schmidt procedure with column pivoting to compute the numerical rank and an orthogonal basis of the numerical range

**Input:** Set of $m_0$ TT-vectors $\mathcal{W} = [\mathcal{W}_1, \mathcal{W}_2, \dots, \mathcal{W}_{m_0}]$, a threshold $\varepsilon_{\text{rank}}$, and the rounding threshold $\varepsilon_{\text{round}}$

**Output:** $m$ almost-orthogonal TT-vectors $\mathcal{Q} = [\mathcal{Q}_1, \mathcal{Q}_2, \dots, \mathcal{Q}_m]$ and $\kappa(\mathcal{W}) = r_{1,1}/r_{m,m}$

1. Initialize $\mathcal{Q} = \mathcal{W}$

2. for $i = 1, \dots, m_0$ do

3. Compute the index $i_{max}$ and the value of the maximum norm of the TT-vectors norm; that is, $\max_{j=i}^{m_0} \|\mathcal{Q}_j\|_{\mathcal{F}}$

4. Perform the column pivoting swap($\mathcal{Q}_i, \mathcal{Q}_{i_{max}}$)

5. $r_{i,i} = \|\mathcal{Q}_i\|_{\mathcal{F}}$

6. if $r_{i,i} \varepsilon_{\text{round}} < \varepsilon_{\text{rank}}$ then

7. $m=i$

8. Break

9. end

10. $\mathcal{Q}_i = 1/r_{i,i} \mathcal{Q}_i$

11. for $j = i + 1, \dots, m_0$ do

12. $r_{i,j} = \langle \mathcal{Q}_i, \mathcal{Q}_j \rangle$

14. $\mathcal{Q}_j = \text{round}(\mathcal{Q}_j - r_{i,j} \mathcal{Q}_i, \varepsilon_{\text{round}})$

15. end

16. end

17. return $\mathcal{Q} = [\mathcal{Q}_1, \dots, \mathcal{Q}_m]$, and $\kappa(\mathcal{W}) = r_{1,1}/r_{m,m}$

## 3.2 Some insight from perturbation theory analysis and stopping criterion

In this section, we derive an upper bound on the residual norm when the Rayleigh-Ritz procedure is applied with a non-orthonormal basis and a perturbed matrix-vector product. This setting is similar to that encountered in our TT-FEAST algorithm, where the basis consists of well-conditioned, normalized TT-vectors $[\mathcal{Q}_1, \dots, \mathcal{Q}_{r_\varepsilon}]$.

This basis is used to construct the matrices $\mathbf{A}_{\mathcal{Q}}$ and $\mathbf{B}_{\mathcal{Q}}$, which appear in the reduced generalized eigenvalue problem. Accordingly, their entries are defined as

$$ \mathbf{A}_{\mathcal{Q}}(i,j) = \langle \mathcal{Q}_i, \text{round}(\mathbf{A}\mathcal{Q}_j, \varepsilon_{\text{round}}) \rangle, \quad \mathbf{B}_{\mathcal{Q}}(i,j) = \langle \mathcal{Q}_i, \mathcal{Q}_j \rangle. $$

For the sake of simplicity of notation, the bound is established in the classical matrix framework where we assume that the effect of the rounding introduces a relative normwise perturbation of the vector and denote $\bar{\mathbf{u}}$ the rounded value of $\mathbf{u}$, so that $\bar{\mathbf{u}} = \mathbf{u} + \delta_{\mathbf{u}}$ with $\|\delta_{\mathbf{u}}\|_2 \le \delta \|\mathbf{u}\|_2$. The columns of the matrix $\mathbf{Q} \in \mathbb{R}^{n \times m}$ is the set of normalized vectors computed by the Modified Gram Schmidt procedure with pivoting. The rounding, that is applied on the individual vector $\mathbf{A}\mathbf{q}_j$ to form the columns of $\mathbf{A}\mathcal{Q}$, introduces a perturbation $\delta_{\mathbf{A}\mathbf{q}_j}$ so that $\overline{\mathbf{A}\mathbf{q}_j} = \mathbf{A}\mathbf{q}_j + \delta_{\mathbf{A}\mathbf{q}_j}$ with

$$ \|\delta_{\mathbf{A}\mathbf{q}_j}\|_2 \le \delta \|\mathbf{A}\mathbf{q}_j\|_2 \le \delta \|\mathbf{A}\|_2, \quad (7) $$

because $\|\mathbf{q}_i\|_2 = 1$. The matrices involved in the reduced generalized eigenvalue problems are

$$ \mathbf{A}_{\mathcal{Q}} = \mathbf{Q}^T \overline{\mathbf{A}} \mathbf{Q} = \mathbf{Q}^T \mathbf{A} \mathbf{Q} + \mathbf{Q}^T [\delta_{\mathbf{A}\mathbf{q}_1}, \dots, \delta_{\mathbf{A}\mathbf{q}_m}] = \mathbf{Q}^T \mathbf{A} \mathbf{Q} + \mathbf{C} \quad (8) $$

where $\mathbf{C} = \mathbf{Q}^T \mathbf{\Delta}_{\mathbf{A}\mathcal{Q}}$ with $\mathbf{\Delta}_{\mathbf{A}\mathcal{Q}} = [\delta_{\mathbf{A}\mathbf{q}_1}, \dots, \delta_{\mathbf{A}\mathbf{q}_m}]$, and

$$ \mathbf{B}_{\mathcal{Q}} = \mathbf{Q}^T \mathbf{Q}. $$

Let us denote $(\lambda, \mathbf{g})$ the Ritz pairs of the reduced generalized eigenvalue problem

$$ \mathbf{A}_{\mathcal{Q}} \mathbf{g} = \lambda_{\mathcal{Q}} \mathbf{B}_{\mathcal{Q}} \mathbf{g} \quad (9) $$
so that the approximated pairs of the original problems are $(\lambda_Q, Qg)$ with $\|g\|_2 = 1$. Let us now examine the residual of these approximations with respect to the original problem $Au = \lambda u$. We have:

$$r = AQg - \lambda_Q Qg.$$

Let us multiply the residual by $Q^T$ and use Equation (8) to get

$$Q^T r = (A_Q - C)g - \lambda_Q B_Q g = -Cg.$$

If we take the norm of the two sides of the equality and use properties of the 2-norm of linear operators we have

$$\sigma_{\min}(Q)\|r\|_2 \le \|Q^T r\|_2 = \|Cg\|_2 \le \|C\|_2. \quad (10)$$

We now seek to derive a bound for $\|C\|_2$ with respect to $\delta$. By definition of $C$, we have

$$\|C\|_2 \le \|Q^T\|_2 \|\Delta_{AQ}\|_2.$$

For any vector $y \in \mathbb{R}^m$, we can write

$$\Delta_{AQ} y = [\delta_{Aq_1}, \dots, \delta_{Aq_m}] y = \sum_{j=1}^{m} y_j \delta_{Aq_j}.$$

Therefore,

$$\|\Delta_{AQ} y\|_2 \le \sum_{j=1}^{m} \|y_j \delta_{Aq_j}\|_2 \le \delta \|A\|_2 \sum_{j=1}^{m} |y_j| \le \delta \sqrt{m} \|A\|_2 \|y\|_2$$

where the Cauchy-Schwarz inequality is used in the last step. It follows that the 2-norm of $\Delta_{AQ}$ satisfies

$$\|\Delta_{AQ}\|_2 = \max_{\|y\|_2=1} \{\|\Delta_{AQ} y\|_2\} \le \delta \sqrt{m} \|A\|_2.$$

Finally, we obtain the bound

$$\|C\|_2 \le \|Q^T\|_2 \|\Delta_{AQ}\|_2 \le \delta \sqrt{m} \|Q\|_2 \|A\|_2. \quad (11)$$

Let us summarize the results of the bounds (10) and (11) in the following proposition.

**Proposition 1** The residual norm of the eigenpair $(\tilde{\lambda}, \tilde{x})$, obtained from a subspace spanned by the $m$ normalized columns $q_i$ of the matrix $Q$ via the Rayleigh-Ritz procedure, when using a matrix-vector product affected by a normwise perturbation $\delta$, satisfies:

$$\|A\tilde{x} - \tilde{\lambda}\tilde{x}\|_2 \le \delta \sqrt{m} \kappa(Q) \|A\|_2,$$

where $\kappa(Q)$ denotes the 2-norm condition number of $Q$.

In the context of TT computations, the results of Proposition 1 provide insight into the attainable accuracy of the eigenpairs computed by TT-FEAST. Since the dimension of the search space is adapted to approximate the numerical range of the projection, we typically have

$$\kappa(Q) \approx \mathcal{O}(1).$$

Moreover, once TT-FEAST successfully identifies an invariant subspace associated with the eigenvalues in the interval $[\lambda_\ell, \lambda_u]$, the spectral norm $\|A\|_2$ restricted to this subspace is bounded above by $\lambda_u$. Consequently, a practical rule of thumb for estimating the attainable accuracy of TT-FEAST is

$$\|A\chi - \lambda\chi\|_F \le \varepsilon_{\text{round}} \sqrt{m} \lambda_u,$$

as a rounding is applied to the computed TT-vectors when forming the Rayleigh quotient matrix in the TT-calculation.

To enhance the robustness of the stopping criterion, this residual-based measure is complemented by monitoring the variation in the sum of the $m_0$ Ritz values between two consecutive iterations. This variation should remain below a prescribed threshold, enabling the detection of any new eigenvalue entering the target interval between two successive iterations, as proposed in [18]. At iteration $k$ of the TT-FEAST algorithm, this criterion takes the form:

$$|\text{trace}(\Lambda^{(k)}) - \text{trace}(\Lambda^{(k-1)})| \le \varepsilon_{\text{trace}} \quad \text{and} \quad \max_{i=1,\dots,m_0} \|A\chi_i^{(k)} - \lambda_i^{(k)}\chi_i^{(k)}\|_F \le \varepsilon_{\text{round}} \sqrt{m} \lambda_u. \quad (12)$$
## 3.3 Memory growth management policy using a continuation procedure

At the beginning of TT-FEAST, when random initial TT-vectors (i.e., random cores of rank 1) are used, the TT-rank of some of the involved TT-vectors can potentially become significantly large before the convergence starts to take place. We observed that this phenomenon arises particularly during the solution of the linear systems using AMEN in the projector construction process (Step 3 of Algorithm 2). This occurs before convergence of TT-FEAST to the invariant subspace is achieved and can lead to excessive memory consumption, making the algorithm infeasible for large-scale problems in high dimensions (even on twelve dimensions, in our largest test case). To address this issue, we introduce a continuation procedure for the convergence threshold, which in turn influences the TT-rounding threshold.

The proposed process starts with a sufficiently large threshold $\epsilon_{round} = \epsilon_{Amen}$ to ensure that the TT-ranks of the rounded TT-vectors remain small until the eigenpairs converge with this low-accuracy threshold. These approximations are then used as the initial guess for the next outer TT-FEAST iteration, where both the TT-rounding ($\epsilon_{round}$) and the convergence threshold ($\epsilon_{round}\sqrt{m\lambda_u}$) are gradually reduced. The continuation procedure is detailed in Algorithm 4.

<caption>**Algorithm 4:** Continuation loop for the convergence threshold</caption>

Initialize $\epsilon_{cont} = [\epsilon_1, \dots, \epsilon_s]$ with $\epsilon_i > \epsilon_{i+1}$, and $\mathcal{X}$ as a set of $m_0$ random TT-vectors
for $\epsilon$ in $\epsilon_{cont}$ do
$\mathcal{Y} = \mathcal{X}$
$\Lambda, \mathcal{X} = \text{TT-FEAST}(\lambda_\ell, \lambda_u, \mathcal{A}, \mathcal{Y}, \text{MaxIter}, \epsilon)$
end
return $\Lambda$ and $\mathcal{X}$

# 4 Numerical results

## 4.1 Numerical test cases

The algorithm is implemented using the ttpy library, the Python version of TT-Toolbox [15]. The AMEN solver [6] implemented in the ttpy library is used to solve the linear system in TT-format in line 3 of Algorithm 2. The numerical tests were performed using Intel Python version 3.7.9.

To illustrate the performance of our eigensolver, we consider two test cases. The first is the academic problem of the d-dimensional Laplacian eigenproblem, for which the eigenpairs are known analytically. This allows us to assess performance metrics such as convergence speed. The second example comes from a collaboration with the computational chemistry community and involves the vibrational Hamiltonian [14] of two molecules, namely H₂CO and CH₃CN.

**The d-dimensional Laplacian** To gain insight into the numerical properties of the algorithm, we consider the d-dimensional Laplacian, for which we know the analytic expression of the eigenvalues and eigenvectors, the latter being of low rank. We consider the eigenvalue problem $\mathcal{L}\mathcal{X} = \lambda\mathcal{X}$ with Dirichlet boundary conditions, where $\mathcal{L} = -\Delta$ is the Laplacian operator discretized on a grid $n_1 \times n_2 \times \dots \times n_d$ in tensor form as

$$ \mathcal{L} = -\Delta = \mathbf{D}_{n_1} \otimes \mathbb{I}_{n_2} \otimes \mathbb{I}_{n_3} \otimes \dots \otimes \mathbb{I}_{n_d} + \mathbb{I}_{n_1} \otimes \mathbf{D}_{n_2} \otimes \mathbb{I}_{n_3} \otimes \dots \otimes \mathbb{I}_{n_d} + \dots + \mathbb{I}_{n_1} \otimes \mathbb{I}_{n_2} \otimes \mathbb{I}_{n_3} \otimes \dots \otimes \mathbf{D}_{n_d}, $$

where $\mathbb{I}_{n_i}$ is the $n_i \times n_i$ identity matrix, and $\mathbf{D}_{n_i}$ is tridiagonal(-1, 2, -1) is the finite difference scheme of the one-dimensional Laplace operator on a regular grid with $n_i$ points. The operator $\mathcal{L}$ is then represented in TT-format with TT-ranks equal to two [9].

**The vibrational Hamiltonian** We apply our method to a more challenging physical problem: the calculation of the infrared spectrum of the Formaldehyde molecule, H₂CO (leading to a 6-dimensional
representation) and Acetonitrile molecule, CH₃CN (leading to a 12-dimensional representation). This requires solving for the eigenvalues of the vibrational Hamiltonian operator H:

$$H = H_0 + V,$$

where H₀ represents the harmonic Hamiltonian, and V is the potential energy surface operator as described in [3, 7]. The problem is then discretized using the eigenvectors of the operator H₀, which are products of 1-d Hermite functions. Using a variational approach, we reformulate the problem as a classical eigenvalue problem in tensor format. For more details on this approach, the reader is referred to [7, 11, 14], in particular for the coefficients of the harmonic Hamiltonian and the potential energy. The maximal TT-rank for the TT representation of the Hamiltonian operator of the molecules H₂CO and CH₃CN for a rounding of 10⁻⁹, is 14 and 26, respectively.

By examining both a simpler, analytically tractable problem (the Laplacian) and a higher-dimensional, application-driven system (the vibrational Hamiltonian), we aim to demonstrate the numerical capabilities of our TT-based eigensolver. In the following, we present detailed results on the accuracy, convergence, and computational performance observed in these two settings.

## 4.2 Impact of the number of Gauss points for the contour integration

In this first experiment, we examine how the number of Gauss points p used to discretize the complex contour affects both the convergence rate and accuracy of the TT-FEAST algorithm. We focus on the four-dimensional Laplacian with grid sizes n = [10, 12, 9, 11]. Let us consider the interval [180.0, 194.0], as we know explicitly the spectrum, it contains m = 9 eigenpairs. At each FEAST iteration, we solve pm₀ (with m₀ = m) linear systems in the TT format, keeping both the AMEN solver tolerance and the tensor rounding threshold at εAmen = εround = 10⁻⁸. The number of iterations of TT-FEAST is set to 6 and m = m₀ does not change during the iterations. Figure 2 shows the evolution of the errors on the eigenvalues (Figure 2a), eigenvectors (Figure 2b), and the maximum of the residual norms (Figure 2c), over six iterations for various number of Gauss points p. The eigenvalue error is the maximum of relative errors, that is

$$err_{\lambda} = \max_{i} \frac{|\lambda_{i} - \lambda_{i}^{*}|}{\lambda_{i}^{*}},$$

where λᵢ* is the iᵗʰ theoretical eigenvalue. The quality of the eigenvectors is assessed by measuring the sine of the angle between the true and computed eigenvectors as follows

$$\sin(\chi, \chi^*) = \sqrt{|1 - \langle \chi, \chi^* \rangle|^2},$$

where χ* is the theoretical eigenvectors.

We observe that increasing the number of discretization points, p, along the contour generally enables a faster convergence, as a finer contour discretization yields a better preconditioner. Indeed, when p ≥ 8, the method typically converges in about two or three TT-FEAST iterations. On the other hand, a larger p increases the computational expense, since each iteration must handle more linear systems. However, it does not affect the quality of the solution. Consequently, there is a trade-off between fast convergence and high iteration costs. As in the matrix case, we notice that the eigenvalues (Figure 2a) converge quadratically compared to the eigenvectors (Figure 2b) for a given accuracy level, reflecting the consequence of the Rayleigh quotient procedure.

Based on these observations, choosing p = 8 provides a good trade-off between the precision of the calculations and the calculation time (number of systems to solve). Unless otherwise stated, we therefore adopt p = 8 in subsequent experiments.

The εAmen parameter is not important as long as it is not too far from εround. Unless explicitly stated otherwise, we set the accuracy of the AMEN solver, εAmen, to εround.

## 4.3 Attainable accuracy versus rounding threshold

In this section, we explore the effect of the rounding threshold εround on both the residual accuracy and the computational efficiency of the TT-FEAST algorithm. To illustrate the theoretical arguments
<figure><img src="image_5.png" /><figcaption>Figure 2: 4-dimensional Laplacian – The errors and maximum residual norm of eigenpairs with iterations for different Gauss points $p$, with $n = [10, 12, 9, 11]$, $m_0 = m = 9$, and $\varepsilon_{\text{round}} = \varepsilon_{\text{Amen}} = 10^{-8}$.</figcaption></figure>

described in Section 3.2, we carried out numerical tests on a discretized Laplacian operator using the TT-FEAST algorithm in matrix format. We considered a mesh of size $n = [4, 6, 3, 7]$, seeking the $m = 8$ eigenvalues in the interval $[355.15, 393.27]$ so that $\lambda_u = 393.27$. At each iteration, when we compute $Aq_i$ to form $AQ$ we add a random perturbation vector of norm $\delta \lambda_u$, with $\delta = 10^{-6}$. We conduce two sets of experiments. With a fixed perturbation magnitude $\delta$, we first investigate the behavior of the residual norm and its bound given in Proposition 1 that describes the impact of the smallest singular value of the search space basis on the attainable accuracy. A perfectly (up-to machine precision) orthonormal basis of the projection space is first computed using an Householder QR factorization of $W$. Based on it, a new basis $\tilde{Q}$ of the space spanned by $Q$ is then computed with a worth condition number but still normalized vectors. Figure 3 shows the convergence history of the largest residual-norm eigenpair for different moderate values of $\kappa(\tilde{Q})$. For large condition numbers, our analysis no longer applies, as the computed residual of the generalized eigenproblem (9) begins to significantly deviate from zero, introducing an additional source of error. The horizontal dashed lines indicate the value of the upper bound in Proposition 1. It is worth noting that the bound captures the general behavior, but is not necessarily tight.

We now illustrate the effect of the perturbation norm $\delta$ on the attainable accuracy in the case $\kappa(Q) = 1$, i.e., when $Q$ is retained as computed by the Householder transformation. We display in Figure 4, the convergence history of the largest residual norm along the iterations for different values of $\delta$. We observe that the residual standards steadily decrease to below the threshold of $\sqrt{m\delta}\lambda_u$, at which point they cease to improve.

In order to illustrate that the analogy between perturbation in matrix computation and TT-computation is relevant, we report in Figure 5 the convergence history of the largest residual norms for a 6d Hamiltonian problem ($H_2CO$ molecule) discretized with $n = [4, 4, 3, 3, 2, 2]$ for two intervals of sought eigenvalues,
<figure><img src="image_6.png" /><figcaption>Figure 3: 4-dimensional Laplacian with $n = [4, 6, 3, 7]$, $m_0 = m = 8$ – The residual norms versus iterations for different values of $\kappa(\tilde{Q})$ and $\delta = 10^{-6}$, dashed lines represent the upper bound value in Proposition 1.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4: 4-dimensional Laplacian with $n = [4, 6, 3, 7]$, $m_0 = m = 8$ – The residual norms versus iterations for different values of $\delta$ in matrix format with normwise random perturbation on $Aq$ to compute $A_Q$ and $\kappa(Q) = 1$.</figcaption></figure>

namely $I_1 = [8620, 8800]$ and $I_2 = [6900, 75300]$ that both contain $m = 4$ eigenvalues. We recall that, in our matrix computation perturbation experiment, only the matrix-vector product in the Rayleigh-Ritz extraction is perturbed, while the orthonormal basis computed via the Householder transformation is preserved (i.e., $\kappa(Q) = 1$ up to machine precision). It is observed that, in terms of attainable accuracy, TT-FEAST and randomly perturbed FEAST in matrix computations yield comparable residual norms. This tends to indicate that the main source of inaccuracy in the TT-computation is the Rayleigh-Ritz extraction, with the rounding in computing $A_Q$ having the most influence on the final accuracy, while other rounding operations seem to have no dominant impact.

We present simulation results for H₂CO in Figure 6, where we examine the impact of varying the rounding threshold $\varepsilon_{round}$. As previously observed in the convergence history, the attainable accuracy reaches the expected order $\sqrt{m\lambda_u\varepsilon_{round}}$ for all values of $\varepsilon_{round}$, except for $\varepsilon_{round} = 10^{-8}$, where the theoretical bound appears overly pessimistic. This lack of tightness can be understood by analyzing the maximum rank of the approximated eigenvectors (see Figure 6b). Specifically, for $\varepsilon_{round} = 10^{-8}$, rounding does not alter the rank, implying that the rounded vector remains identical to the original one. As a result, no perturbation is introduced by the rounding process, whereas the perturbation theory analysis assumes a systematic relative perturbation of all vectors with an amplitude of $\varepsilon_{round} = 10^{-8}$.

## 4.4 Auto-tuning of the projection space dimension

In this section, we show how our adaptive strategy explained in the Section 3.1 allows us to automatically determine the correct number of eigenvalues in the interval while limiting the memory footprint.
<figure><img src="image_8.png" /><figcaption>Figure 5: H₂CO vibrational Hamiltonian – Comparison of TT-FEAST with FEAST in matrix format with perturbation with $n = [4, 4, 3, 3, 2]$, $m_0 = 10$, $m = 4$, $I_1 = [8620, 8800]$ or $I_2 = [6900, 75300]$ and $\epsilon_{round} = \delta$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 6: H₂CO vibrational Hamiltonian – The residual norms and ranks versus iterations for different $\epsilon_{round}$, within the interval $I = [8620, 8800]$, with $n = [5, 5, 4, 4, 3, 3]$, $m_0 = 5$, $m = 4$, and $\epsilon_{Amen} = 10^{-8}$.</figcaption></figure>

As before, we consider the four-dimensional Laplacian in the interval [180, 194], which is known to contain $m = 9$ eigenvalues. The grid discretization is $n = [10, 12, 9, 11]$, both the AMEN solver and tensor rounding thresholds are set to $\epsilon_{Amen} = \epsilon_{round} = 10^{-8}$. Notably, we choose $m_0 = 8$, which is strictly less than the actual number of eigenpairs $m = 9$, to illustrate how the dimension of the subspace is increased during iterations.

Figure 7a illustrates the ratio $r_{i,i}/r_{1,1}$ of the TT-vectors $\mathcal{W}$ during the first two iterations of TT-FEAST when $m_0$ is insufficiently large. Ideally, this ratio should remain close to one if $\mathcal{Q}$ is close to an orthonormal basis, which is observed at the second iteration. However, since no decrease in this ratio is detected, the search space is expanded by adding four random TT-vectors for the third iteration. As shown in Figure 7b, during the third iteration, all ratios—except for the first—drop below one, with a rapid decline after the eighth vector. In the fourth iteration, the first nine ratios return to one, while a significant decay is observed at the tenth vector. Consequently, the dimension detection algorithm reduces the search space to nine—the number of eigenvalues within the interval—at the fifth iteration. At this stage, $\mathcal{Q}$ forms an almost orthonormal basis, as all ratios are close to one.

Figure 7c presents the residual norms associated with the eigenpairs computed at each iteration, using the same color to represent each iteration. The vertical light gray dashed lines indicate the locations of the analytical eigenvalues within the interval, which is bounded by vertical red lines. In the first two iterations (black and red), only 8 out of the 9 eigenvalues are approximated, and those at the extremes of the interval begin to converge. However, the corresponding eigenvectors remain poorly estimated, as indicated by the large residual norms. At the third iteration (blue), when the search space is expanded,
<figure><img src="image_10.png" /><figcaption>Figure 7: 4-dimensional Laplacian – Illustration of the adaptive strategy for $m_0 < m$, with $n=[10, 12, 9, 11]$, $m_0 = 8$, $\varepsilon_{\text{round}} = \varepsilon_{\text{Amen}} = 10^{-8}$, and $I = [180, 194]$.</figcaption></figure>

all eigenvalues within the interval are more accurately approximated, though the residual norms remain relatively high. Since the ratio $r_{i,i}/r_{1,1}$ had not yet strictly deviated from one, 12 approximated eigenpairs were computed at this stage. By the fourth iteration (green), a clear decay of $r_{i,i}/r_{1,1}$ is observed beyond the ninth vector, leading to a reduction in the search space dimension to nine, with only nine eigenpairs being approximated. While the eigenvalues are now well captured, the eigenvectors still exhibit residual norms exceeding the theoretical bound established in Section 3.2. Finally, at the last iteration (yellow), the accuracy of the nine eigenpairs is further improved, and all satisfy the stopping criterion.

Overall, these results illustrate the advantage of the orthogonalization procedure (Algorithm 3), which provides a reliable means to determine the $\varepsilon$-rank of the projector and consequently the number of eigenpairs inside the interval. It enables us to adapt the dimension of the search space. Note that the number of additional vectors is free to be chosen according to the size of the initial subspace.

## 4.5 Performance and advantages of $\varepsilon_{\text{round}}$ the continuation procedure

In high-dimensional settings, a key challenge is the substantial storage required for TT-vectors, whose ranks can be large before the convergence occurs. This happens especially during the first iterations when the projectors are constructed from random vectors. The ranks of the projectors can be very large, especially when the TT-matrix has high ranks. To mitigate this, we introduce a continuation procedure with a sequence of decreasing convergence threshold values (consequently TT-rounding threshold values), allowing the computed eigenvectors to gradually reach the desired accuracy. This strategy prevents excessive memory consumption by avoiding the need to store high-accuracy TT-vectors in the early stages of convergence.
To measure the memory impact, we define the compression ratio of a tensor as the ratio of the storage requirements of its TT-vectors to its storage in dense form. For a TT-vector $\mathcal{A} \in \mathbb{R}^{n_1 \times \dots \times n_d}$ of TT-rank $r_i$, the compression ratio is given by $\sum_{i=1}^d r_{i-1} n_i r_i / \prod_{j=1}^d n_j$. For a set of TT-vectors $[\mathcal{A}_1, \dots, \mathcal{A}_m]$, we consider the mean compression ratio defined as $1/m \sum_{i=1}^m \text{ratio}(\mathcal{A}_i)$.

<figure><img src="image_11.png" /><figcaption>Figure 8: 4-dimensional Laplacian – Overview of the continuation procedure starting with $\varepsilon_1$ as first convergence and rounding threshold, with $n = [32, 34, 31, 33]$, $m_0 = 6$, and $I = [180, 194]$.</figcaption></figure>

We first consider the 4-dimensional Laplacian problem with a target accuracy corresponding to $\varepsilon_{round} = 10^{-8}$. Figure 8a presents the convergence history for three cases: a fixed $\varepsilon_{round} = 10^{-8}$ (blue) and two continuation strategies with decreasing threshold sequences, $\varepsilon_{round} = [10^{-1}, 10^{-3}, 10^{-5}, 10^{-7}, 10^{-8}]$ (green) and $\varepsilon_{round} = [10^{-2}, 10^{-4}, 10^{-6}, 10^{-8}]$ (yellow). In the latter cases, we alternate between solid and dashed lines to indicate transitions in $\varepsilon_{round}$, occurring when convergence is achieved for the previous value. As illustrated in Figure 8a, the first two iterations are performed with the threshold $\varepsilon_1$, after which the continuation threshold is adaptively modified at each iteration. Although the continuation approach requires more iterations to converge, it is important to note that the intermediate iterations (those with $\varepsilon_{round}$ larger than the final target accuracy) are significantly computationally cheaper. This is due to smaller TT-ranks and a larger $\varepsilon_{Amen}$, leading to reduced computational cost and lower peak memory usage. This memory advantage is clear in Figure 8b, where the fixed $\varepsilon_{round} = 10^{-8}$ case is several orders of magnitude more expensive compared to the continuation-based strategies. We deliberately do not provide a quantitative evaluation of time-to-solution savings as our implementation is not optimized for performance in this regard. Nevertheless we can report that the saving is significant compared to constant $\varepsilon_{round}$. To better demonstrate the advantages of the continuation procedure, we examine a

<figure><img src="image_12.png" /><figcaption>Figure 9: CH$_3$CN vibrational Hamiltonian – The eigenpair residuals and compression ratios of the set of TT-vectors over iterations of the continuation procedure with $n = [18, 18, 17, 13, 12, 12, 10, 7, 7, 5, 5, 5]$, $m_0 = 4$, $\varepsilon_1 = 10^{-1}$, and $I = [12437.42, 12457.42]$.</figcaption></figure>

12-dimensional problem related to the CH$_3$CN molecule. We present results only for experiments using the continuation setting $\varepsilon_{round} = [10^{-1}, 10^{-2}, 10^{-3}]$, as a test with a fixed $\varepsilon_{round} = 10^{-3}$ ran out of memory. The results, shown in Figure 9, indicate that the continuation approach effectively manages the
memory storage. The final increase in memory usage is due to the eigenvectors not having a particularly
low rank, unlike in the case of the Laplacian. As a result, representing them with higher accuracy incurs
a larger memory footprint.

# 5 Concluding remarks

In this work, we have demonstrated how the FEAST algorithm can be extended to tackle high-
dimensional eigenvalue problems, up to twelve dimensions, within a low-rank tensor framework. A
central difficulty in this context lies in the absence of an estimator for the number of eigenvalues inside
a given interval, unlike in the standard matrix-based formulation. To overcome this limitation, we
developed an adaptive procedure that dynamically adjusts the dimension of the projection subspace by
incorporating a rank-revealing extension of the Modified Gram-Schmidt algorithm with pivoting for TT-
vectors. Furthermore, we established an explicit bound on the attainable residual norm as a function of
the compression threshold and the number of targeted eigenpairs, leading to a robust stopping criterion
for the TT-FEAST iterations. A continuation strategy was also introduced to progressively refine the
tensor representation, effectively mitigating the initial memory overhead caused by the high TT-rank of
randomly generated projectors. Numerical experiments on problems up to twelve dimensions confirm
that TT-FEAST is a possible alternative to existing high-dimensional eigensolvers, particularly when
only a limited number of eigenvalues in a specific region of the spectrum are sought.

In conclusion, the primary objective of this paper was to demonstrate the adaptability of FEAST
to the tensor setting and to address the specific challenges that arise in this context. A comparative
study with alternative solvers, aimed at identifying the regimes in which each method performs best,
would be a valuable direction for future research. Additionally, all numerical experiments were carried
out in Python on a standard computing platform. It would be interesting to investigate scalability on
very high-dimensional and large-scale problems to better understand the performance and limitations of
TT-FEAST, which may require the development of a parallel implementation and could open further
research opportunities.

# Acknowledgment

The authors are grateful to the referees for their thoughtful and supportive feedback, which con-
tributed to strengthening the manuscript. Experiments presented in this paper were carried out using
the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bor-
deaux, Bordeaux INP and Conseil Régional d'Aquitaine (see https://www.plafrim.fr).

# References

[1] Ravi P. Agarwal, Kanishka Perera, and Sandra Pinelas. *An Introduction to Complex Analysis*. Springer, Boston, MA, 2011.

[2] Gustavo Avila and Tucker Carrington. "Using nonproduct quadrature grids to solve the vibrational Schrödinger equation in 12D". In: *Journal of Chemical Physics* 134.5 (2011). ISSN: 00219606. DOI: 10.1063/1.3549817.

[3] Tucker Carrington. "Iterative Methods for Computing Vibrational Spectra". In: *Mathematics* 6.1 (2018). ISSN: 2227-7390. DOI: 10.3390/math6010013.

[4] Edoardo Di Napoli, Eric Polizzi, and Yousef Saad. "Efficient estimation of eigenvalue counts in an interval". In: *Numerical Linear Algebra with Applications* 23.4 (Aug. 2016), pp. 674–692. ISSN: 10705325. DOI: 10.1002/nla.2048.

[5] Sergey V. Dolgov, Boris N. Khoromskij, Ivan V. Oseledets, and D.V. V. Savostyanov. "Computation of extreme eigenvalues in higher dimensions using block tensor train format". In: *Computer Physics Communications* 185.4 (Apr. 2014), pp. 1207–1216. ISSN: 00104655. DOI: 10.1016/j.cpc.2013.12.017.
[6] Sergey V. Dolgov and Dmitry V. Savostyanov. "Alternating minimal energy methods for linear systems in higher dimensions". In: *SIAM Journal on Scientific Computing* 36.5 (2014), A2248–A2271. ISSN: 10957200. DOI: 10.1137/140953289.

[7] Romain Garnier, Marc Odunlami, Vincent Le Bris, Didier Bégué, Isabelle Baraille, and Olivier Coulaud. "Adaptive vibrational configuration interaction (A-VCI): A posteriori error estimation to efficiently compute anharmonic IR spectra". In: *The Journal of Chemical Physics* 144.20 (May 2016), p. 204123. ISSN: 0021-9606. DOI: 10.1063/1.4952414.

[8] Gene H. Golub and Charles F. van Loan. *Matrix Computations*. Fourth. JHU Press, 2013. ISBN: 1421407949 9781421407944.

[9] Vladimir A. Kazeev and Boris N. Khoromskij. "Low-Rank Explicit QTT Representation of the Laplace Operator and Its Inverse". In: *SIAM Journal on Matrix Analysis and Applications* 33.3 (Jan. 2012), pp. 742–758. DOI: 10.1137/100820479.

[10] James Kestyn, Vasileios Kalantzis, Eric Polizzi, and Yousef Saad. "PFEAST: A High Performance Sparse Eigenvalue Solver Using Distributed-Memory Linear Solvers". In: *SC16: International Conference for High Performance Computing, Networking, Storage and Analysis*. IEEE, Nov. 2016, 178–189. DOI: 10.1109/sc.2016.15.

[11] Vincent Le Bris, Marc Odunlami, Didier Bégué, Isabelle Baraille, and Olivier Coulaud. "Using computed infrared intensities for the reduction of vibrational configuration interaction bases". In: *Phys. Chem. Chem. Phys.* 22 (13 2020), pp. 7021–7030. DOI: 10.1039/D0CP00593B.

[12] Arnaud Leclerc and Tucker Carrington. "Calculating vibrational spectra with sum of product basis functions without storing full-dimensional vectors or matrices." In: *The Journal of chemical physics* 140.17 (2014), p. 174111. ISSN: 1089-7690. DOI: 10.1063/1.4871981.

[13] Steven J. Leon, Åke Björck, and Walter Gander. "Gram-Schmidt orthogonalization: 100 years and more". In: *Numerical Linear Algebra with Applications* 20.3 (June 2012), pp. 492–532. DOI: 10.1002/nla.1839.

[14] Marc Odunlami, Vincent Le Bris, Didier Bégué, Isabelle Baraille, and Olivier Coulaud. "A-VCI: A flexible method to efficiently compute vibrational spectra". In: *Journal of Chemical Physics* 146.21 (June 2017). DOI: 10.1063/1.4984266.

[15] I. V. Oseledets. *ttpy: Python implementation of the TT-Toolbox*. https://github.com/oseledets/ttpy. Feb. 2019.

[16] I. V. Oseledets and S. V. Dolgov. "Solution of Linear Systems and Matrix Inversion in the TT-Format". In: *SIAM Journal on Scientific Computing* 34.5 (2012), A2718–A2739. DOI: 10.1137/110833142.

[17] Ivan V. Oseledets. "Tensor-Train Decomposition". In: *SIAM Journal on Scientific Computing* 33.5 (2011), pp. 2295–2317. DOI: 10.1137/090752286.

[18] Eric Polizzi. "Density-matrix-based algorithm for solving eigenvalue problems". In: *Physical Review B* 79.11 (Mar. 2009), p. 115112. ISSN: 1098-0121. DOI: 10.1103/PhysRevB.79.115112.

[19] Maxim Rakhuba, Alexander Novikov, and Ivan Oseledets. "Low-rank Riemannian eigensolver for high-dimensional Hamiltonians". In: *Journal of Computational Physics* 396 (2019). DOI: 10.1016/j.jcp.2019.07.003.

[20] Maxim Rakhuba and Ivan V. Oseledets. "Calculating vibrational spectra of molecules using tensor train decomposition". In: *Journal of Chemical Physics* 145.12 (2016). ISSN: 00219606. DOI: 10.1063/1.4962420.

[21] Yousef Saad. *Numerical Methods for Large Eigenvalue Problems*. Classics in Applied Mathematics. Society for Industrial and Applied Mathematics, Jan. 2011. DOI: 10.1137/1.9781611970739.

[22] Ping Tak, Peter Tang, and Eric Polizzi. "FEAST as a subspace iteration eigensolver accelerated by approximate spectral projection". In: *SIAM Journal on Matrix Analysis and Applications* 35.2 (2014), pp. 354–390. DOI: 10.1137/13090866X.
Inria

RESEARCH CENTRE
Centre Inria de l'Université de Bordeaux

200 avenue de la Vieille Tour
33405 Talence Cedex

Publisher
Inria
Domaine de Voluceau - Rocquencourt
BP 105 - 78153 Le Chesnay Cedex
inria.fr

ISSN 0249-6399