# Variational quantum subspace construction via symmetry-preserving cost functions

Hamzat Akande, Alexandre Perrin, Bruno Senjean, Matthieu Saubanère

► To cite this version:

Hamzat Akande, Alexandre Perrin, Bruno Senjean, Matthieu Saubanère. Variational quantum subspace construction via symmetry-preserving cost functions. Physical Review A, 2025, 112 (3), pp.032623. 10.1103/6hpp-zl5h. hal-05296296

HAL Id: hal-05296296

https://hal.science/hal-05296296v1

Submitted on 3 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
Variational Quantum Subspace Construction via Symmetry-Preserving Cost Functions

Hamzat A. Akande,<sup>1,2,3</sup> Alexandre Perrin,<sup>1</sup> Bruno Senjean,<sup>3</sup> and Matthieu Saubanère<sup>1,*</sup>

<sup>1</sup>Université de Bordeaux, CNRS, LOMA, UMR 5798, F-33400 Talence, France

<sup>2</sup>ICTP-East African Institute for Fundamental Research, University of Rwanda, Kigali, Rwanda

<sup>3</sup>ICGM, Univ Montpellier, CNRS, ENSCM, Montpellier, France

Determining low-energy eigenstates in electronic many-body quantum systems is a key challenge in computational chemistry and condensed-matter physics. Hybrid quantum-classical approaches, such as the Variational Quantum Eigensolver and Quantum Subspace Methods, offer practical solutions but face limitations in circuit depth and measurement overhead. In this article, we propose a variational strategy based on symmetry-preserving cost functions to iteratively construct a reduced subspace for the extraction of low-lying energy states. We show that, under certain conditions, our approach leads to a tridiagonal representation similar to that obtained with the Lanczos algorithm. The iterative process allows control over the trade-off between circuit depth, the number of variational parameters, and the number of measurements required to achieve the desired accuracy, making it suitable for current quantum hardware. As a proof of concept, we test the proposed algorithms on $H_4$ chain and ring, targeting both the ground-state energy and the charge gap.

# I. INTRODUCTION

The many-body problem poses significant challenges in determining electronic low-lying eigenenergies using classical methods. Quantum computing offers a promising alternative to address the exponential complexity [1–4]. However, noisy intermediate-scale quantum (NISQ) devices are limited by short qubit coherence times and minimal error correction, making algorithms like phase estimation impractical [5, 6]. Consequently, hybrid approaches such as the Variational Quantum Eigensolver (VQE) and Quantum Subspace Methods (QSM) have emerged as feasible solutions for current quantum hardware [7–11].

On one hand, the VQE algorithm optimizes a variational trial state encoded in a quantum circuit using an ansatz, broadly classified as physically-motivated or hardware-efficient. Physically-motivated ansätze, such as the unitary coupled-cluster (UCC) ansatz [7, 12–21], are symmetry-preserving (for example number of electrons and spin, although the later might still be broken by using trotterized – also called disentangled – UCC [22]) and size-extensive but require deep circuits, which are impractical on current hardware [23–25]. The associated circuit depth can however be reduced using adaptive methods like ADAPT-VQE that selects only impactful excitation operators [26–38]. Despite these improvements, estimating ground-state energies for even small molecules remains challenging due to high gate counts and measurement overhead [26, 32, 39]. Shallower circuit can be constructed using Hardware-efficient ansätze (HEA) that employ layers of parametrized single-qubit rotations and entangling gates, enabling practical implementation on NISQ devices [40, 41]. However, HEA faces scalability challenges due to redundant parameters, lack of symmetry constraints, and barren plateaus, hinder-

ing optimization [42–47]. To address the aforementioned issues, symmetry-preserving ansätze [13, 23, 37, 48–52] and advanced optimization strategies [53–60] have been proposed. Additional strategies, such as adding auxiliary qubits [61], reducing the number of qubits using the quantum matrix product states algorithm [62], or using dissipative quantum algorithms [63] have also been developed recently to mitigate the barren plateau issue. We refer the reader to Refs. [45, 64] for recent reviews.

On the other hand, QSMs project the Schrödinger equation onto a reduced subspace, combining quantum measurements with classical eigenvalue solvers [9, 65]. The distinction among QSM approaches lies in how the reduced subspace is constructed. For instance, Quantum Subspace Expansion (QSE) constructs subspace vectors using excitation operators, inspired by classical methods like Multi-Reference Configuration Interaction Singles and Doubles [9, 65–70], but suffers from size-intensity for large systems [11]. Alternatives to the excitation operators include Pauli type operators [71–73] or the use of the equation-of-motion formalism for excited states [74, 75]. Quantum Krylov-based methods extend classical eigenvalue solvers to quantum systems, employing Hamiltonian powers to construct subspaces [11, 76]. Examples include quantum Lanczos [77–79], Chebyshev Krylov [80], Gaussian power [81], inverse power [82], Davidson [83] and Quantum Filter Diagonalization [10, 84–92]. Quantum Krylov based methods require significantly more quantum resources than QSE approaches, due to the construction of the block-encoded unitaries or the need of precise time evolution. However, just like their classical counterpart they excel in approximating precisely dominant eigenvalues with rapid convergence [89, 91]. For such approaches, synergy between quantum devices and classical solvers remains central, while challenges like noise and sampling overhead persist [72].

In this contribution, we propose a strategy for designing cost functions that enforce symmetry preservation, regardless of the ansatz used to define the unitary trans-

<footnote>* matthieu.saubanere@cnrs.fr</footnote>
formation. To obtain a low-energy eigenstate, we introduce an iterative procedure that allows control over the ratio between circuit depth, the number of variational parameters, and the number of measurements required to achieve the desired accuracy. In the spirit of QSMs, the iterative process is employed to construct variationally a reduced subspace from which low-lying energy states are extracted. Moreover, we show that the proposed strategy leads to a tridiagonal form analogous to that obtained with the Lanczos algorithm. As a proof of concept, we test the proposed algorithms on H₄ linear and square molecules, targeting both the ground-state energy and the charge gap. Convergence properties in terms of circuit depth and the number of iterations are also analyzed.

## II. DENSITY-MATRIX BASED VARIATIONALS PRINCIPLES

Let us consider the electronic Hamiltonian operator **H** of a system for which the ground state is not degenerate. The Hamiltonian is expanded in the many-body Fock-space basis set {$|\Phi_i\rangle$},

$$ \mathbf{H} = \sum_{ij} H_{ij} |\Phi_i\rangle\langle\Phi_j|, \qquad (1) $$

where {$|\Phi_i\rangle$} denote Slater determinants or configuration state functions, and we refer to $|\Phi_0\rangle$ as the Hartree-Fock (HF) state. **H** can be diagonalized using a unitary transformation matrix, $\bar{\mathbf{H}} = \mathbf{P}^\dagger\mathbf{H}\mathbf{P}$, where $\bar{\mathbf{H}}_{ij} = E_i\delta_{ij}$ and $E_i$ is the i-th eigenvalue sorted in ascending order. The eigenvectors are obtained as $|\Psi_i\rangle = \mathbf{P}|\Phi_i\rangle$, $|\Psi_0\rangle$ and $E_0$ thus correspond to the ground-state vector and energy, respectively. A standard approach to reach the lowest eigenvalue of **H** on a quantum computer relies on the variational principle,

$$ E_0 < E(\boldsymbol{\theta}^*) = \min_{\boldsymbol{\theta}} \langle \Phi_0 | \mathbf{P}^{\dagger}(\boldsymbol{\theta}) \mathbf{H} \mathbf{P}(\boldsymbol{\theta}) | \Phi_0 \rangle, \quad (2) $$

where the unitary transformation $\mathbf{P}(\boldsymbol{\theta})$ depends on parameters $\boldsymbol{\theta}$ to be optimized, and $\boldsymbol{\theta}^*$ denotes the optimal parameters. At the saddle point, $|\Psi_0(\boldsymbol{\theta}^*)\rangle = \mathbf{P}(\boldsymbol{\theta}^*)|\Phi_0\rangle$, and if $\mathbf{P}(\boldsymbol{\theta}^*) = \mathbf{P}$, then $|\Psi(\boldsymbol{\theta}^*)\rangle = |\Psi_0\rangle$ and $E(\boldsymbol{\theta}^*) = E_0$.

Alternatively, we rewrite the variational principle in Eq. (2) in terms of the many-body ground-state density matrices $\bar{\mathbf{\Gamma}}$ and $\mathbf{\Gamma}$ of which elements are defined as

$$ \bar{\mathbf{\Gamma}}_{ij} = \langle \Psi_0 | \Psi_i \rangle \langle \Psi_j | \Psi_0 \rangle, \quad \mathbf{\Gamma}_{ij} = \langle \Psi_0 | \Phi_i \rangle \langle \Phi_j | \Psi_0 \rangle. \quad (3) $$

As depicted in Fig. 1(a), $\bar{\mathbf{\Gamma}}_{ij} = \delta_{i0}\delta_{j0}$ for all i and j. In this framework, the ground-state energy $E_0$ is computed as the following convolution,

$$ E_0 = \mathrm{tr}(\mathbf{\Gamma}\mathbf{H}) = \mathrm{tr}(\bar{\mathbf{\Gamma}}\bar{\mathbf{H}}) = \bar{H}_{00}, \quad (4) $$

and the variational principle in Eq. (2) corresponds to

$$ E_0 < E(\boldsymbol{\theta}^*) = \min_{\boldsymbol{\theta}} \left[ \mathrm{tr}(\bar{\mathbf{\Gamma}}\mathbf{P}^\dagger(\boldsymbol{\theta})\mathbf{H}\mathbf{P}(\boldsymbol{\theta})) \right], \quad (5) $$

<figure><img src="image_2.png" /><figcaption>FIG. 1: (a) Schematic representation of the ground-state many-body density matrices $\mathbf{\Gamma}$, $\bar{\mathbf{\Gamma}}$, and $\tilde{\mathbf{\Gamma}}$ in the three different representations, namely the many-body basis set, eigenvector, and Householder representation. The unitary transformations **P**, **U**, and **R** linking the three representations are displayed. (b) Schematic representation of the two-level decomposition of the ground state. We propose to optimize the gain energy $E_G$ coming from the interaction between the good-guess $|\Phi_0\rangle$ and a variational vector $|\tilde{\Phi}_1(\boldsymbol{\theta})\rangle$.</figcaption></figure>

where tr denotes the trace.

We introduce a third intermediate representation in which the first vector $|\Phi_0\rangle$ remains unchanged, and the density matrix

$$ \tilde{\mathbf{\Gamma}} = \begin{bmatrix} \tilde{\mathbf{\Gamma}}^R & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{bmatrix} \qquad (6) $$

is block-diagonal. $\tilde{\mathbf{\Gamma}}^R$ is a 2 × 2 matrix that reads

$$ \tilde{\mathbf{\Gamma}}^R = \begin{bmatrix} \Gamma_{00} & \tilde{\Gamma}_{10} \\ \tilde{\Gamma}_{01} & \tilde{\Gamma}_{11} \end{bmatrix} = \begin{bmatrix} \omega^2 & \sqrt{\omega^2(1-\omega^2)} \\ \sqrt{\omega^2(1-\omega^2)} & 1-\omega^2 \end{bmatrix}, \quad (7) $$

where $\omega^2 = |\langle\Psi_0|\Phi_0\rangle|^2$ represents the weight of $|\Phi_0\rangle$ in the ground state. This representation is called the Householder representation since, among all unitary transformations **R** that provide the form of $\tilde{\mathbf{\Gamma}}$ in Eq. (6),

$$ \tilde{\mathbf{\Gamma}} = \mathbf{R}\mathbf{\Gamma}\mathbf{R}, \quad (8) $$

with $|\Phi_0\rangle = \mathbf{R}|\Phi_0\rangle$, the Householder reflection ($\mathbf{R} = \mathbf{R}^\dagger$) is uniquely defined from the first column of $\mathbf{\Gamma}$. Note that this transformation has been recently used in the context of embedding methods in electronic structure theory [94–97]. Using this transformation, the ground-state vector is decomposed as

$$ |\Psi_0\rangle = \omega |\Phi_0\rangle + \nu \sqrt{1 - \omega^2} \mathbf{R} |\Phi_1\rangle, \quad (9) $$

where $\nu = \pm 1$ is a relative phase factor, and the ground-state energy reads, in this new representation,

$$ E_0 = \mathrm{tr}(\tilde{\mathbf{\Gamma}}\tilde{\mathbf{H}}) = \mathrm{tr}(\tilde{\mathbf{\Gamma}}^R\tilde{\mathbf{H}}^R), \quad (10) $$
where $\tilde{\mathbf{H}} = \mathbf{RHR}$ and

$$ \tilde{\mathbf{H}}^R = \begin{bmatrix} H_{00} & \tilde{H}_{10} \\ \tilde{H}_{01} & \tilde{H}_{11} \end{bmatrix}. \qquad (11) $$

It follows that the variational principle in Eq. (5) can be rewritten in the Householder representation as

$$ E_0 < E(\boldsymbol{\theta}^*, \omega^*) = \min_{\omega, \boldsymbol{\theta}} \left[ \text{tr} \left( \tilde{\mathbf{\Gamma}}^R(\omega) \tilde{\mathbf{H}}^R(\boldsymbol{\theta}) \right) \right], \quad (12) $$

where $\boldsymbol{\theta}^*$ and $\omega^*$ are the minimizing parameters, and

$$ \tilde{\mathbf{H}}^R(\boldsymbol{\theta}) = \begin{bmatrix} H_{00} & \tilde{H}_{10}(\boldsymbol{\theta}) \\ \tilde{H}_{01}(\boldsymbol{\theta}) & \tilde{H}_{11}(\boldsymbol{\theta}) \end{bmatrix}, \quad (13) $$

with $\tilde{\mathbf{H}}(\boldsymbol{\theta}) = \mathbf{R}^\dagger(\boldsymbol{\theta})\mathbf{HR}(\boldsymbol{\theta})$. By introducing the Householder representation, we have split the variational process into two sub-processes: 1) the optimization of the reflection $\mathbf{R}(\boldsymbol{\theta})$ to block-diagonalize $\tilde{\mathbf{\Gamma}}$ into $\tilde{\mathbf{\Gamma}}$, and 2) the optimization of the unitary $\mathbf{U}(\omega)$ to diagonalize $\tilde{\mathbf{H}}^R(\boldsymbol{\theta})$, i.e. $\mathbf{P} = \mathbf{RU}$, see Fig. (1). Note that since $\tilde{\mathbf{\Gamma}}$ is idempotent, $\mathbf{U}(\omega)$ is only defined by a single parameter $\omega$. Optimizing $\mathbf{R}(\boldsymbol{\theta})$ consists in finding the complementary vector $|\tilde{\Phi}_1(\boldsymbol{\theta})\rangle = \mathbf{R}(\boldsymbol{\theta})|\Phi_1\rangle$ of $|\Phi_0\rangle$ that composes the ground state, as shown in Eq. (9).

### III. ITERATIVE AND VARIATIONAL SYMMETRY-PRESERVING QUANTUM EIGENSOLVER

#### A. Symmetry-preserving cost functions

Let us make the assumption that $|\Phi_0\rangle$ is a binding state i.e. $H_{00} = \langle\Phi_0|\mathbf{H}|\Phi_0\rangle \le 0$ and additionally that it is a good-guess, i.e. that $|\Phi_0\rangle$ as an overlap $\omega^2 = |\langle\Phi_0|\Psi_0\rangle|^2 \ge 0.5$ with the ground state. In that case, instead of using the energy as a cost function to be minimized, we consider the gain energy,

$$ E_G(\boldsymbol{\theta}) = \frac{\Delta(\boldsymbol{\theta})}{2} \left[ 1 - \sqrt{1 + \frac{4\tilde{H}_{01}(\boldsymbol{\theta})\tilde{H}_{10}(\boldsymbol{\theta})}{\Delta(\boldsymbol{\theta})^2}} \right], \quad (14) $$

where $\Delta(\boldsymbol{\theta}) = \tilde{H}_{11}(\boldsymbol{\theta}) - H_{00}$. $E_G(\boldsymbol{\theta})$ is negative if $\Delta(\boldsymbol{\theta}) > 0$ and that is computed in the two-level system composed of $|\tilde{\Phi}_1(\boldsymbol{\theta})\rangle$ and $|\Phi_0\rangle$, as shown in Fig. 1(b). It corresponds to the energy difference between the eigen-state of the two-level system $|\Psi_0\rangle$ and $|\Phi_0\rangle$, i.e. $E_G(\boldsymbol{\theta}) = \langle\Psi_0(\boldsymbol{\theta})|\mathbf{H}|\Psi_0(\boldsymbol{\theta})\rangle - \langle\Phi_0|\mathbf{H}|\Phi_0\rangle$. If $|\Phi_0\rangle$ is a good guess, the ground-state energy $E_0 = \langle\Psi_0|\mathbf{H}|\Psi_0\rangle$ is related with the Hartree-Fock energy $H_{00}$ and the gain energy as

$$ E_0(\boldsymbol{\theta}^*) = H_{00} + E_G(\boldsymbol{\theta}^*), \quad (15) $$

for the optimal parameter $\boldsymbol{\theta}^*$. From a numerical point of view, using $E_G(\boldsymbol{\theta})$ as a cost function for optimization may lead to instability when $\Delta(\boldsymbol{\theta}) \to 0$, i.e., when the

spectrum exhibits degeneracy or when the target state is not the lowest eigenvalue of the spectrum but the lowest-eigenvalue of a specific symmetry-restricted subspace. Alternatively, the interaction energy

$$ E_I(\boldsymbol{\theta}) = -\sqrt{\tilde{H}_{01}(\boldsymbol{\theta})\tilde{H}_{10}(\boldsymbol{\theta})} \quad (16) $$

can be also considered as a relevant cost function. In particular, for large $|\Delta(\boldsymbol{\theta})|$ compared to $|\tilde{H}_{10}(\boldsymbol{\theta})|$, both cost functions are expected to lead to the same saddle point as

$$ \frac{\partial E_G(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} \propto \frac{\partial E_I(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}}. \quad (17) $$

Moreover, we show in the following that minimizing $E_I(\boldsymbol{\theta})$ with the condition that $\langle\Phi_0|\mathbf{R}(\boldsymbol{\theta})\Phi_0\rangle = 1$ corresponds to emptying the first column/row of the Hamiltonian up to the diagonal term and the first off-diagonal element i.e., $\tilde{H}_{0j}(\boldsymbol{\theta}) = \sum_{kl} R_{0k}(\boldsymbol{\theta})H_{kl}R_{lj}(\boldsymbol{\theta}) = 0, \forall j > 1$, similarly as one obtains by using the corresponding Householder Reflection.

**Proof:** Let us consider the first diagonal element of the square of the Hamiltonian in the Householder representation $\tilde{\mathbf{H}}$ and in the variational representations $\tilde{\mathbf{H}}(\boldsymbol{\theta})$,

$$ \begin{align} [\tilde{\mathbf{H}}^2]_{00} &= H_{00}^2 + \tilde{H}_{01}\tilde{H}_{10}, \tag{18} \\ [\tilde{\mathbf{H}}(\boldsymbol{\theta})^2]_{00} &= \tilde{H}_{00}(\boldsymbol{\theta})^2 + \tilde{H}_{01}(\boldsymbol{\theta})\tilde{H}_{10}(\boldsymbol{\theta}) \nonumber \\ &\quad + \sum_{j>1} \tilde{H}_{0j}(\boldsymbol{\theta})\tilde{H}_{j0}(\boldsymbol{\theta}). \tag{19} \end{align} $$

Because the condition $\langle\Phi_0|\mathbf{R}(\boldsymbol{\theta})|\Phi_0\rangle = 1$ imposes that $\tilde{H}_{00} = H_{00}$ and $[\tilde{\mathbf{H}}^2]_{00} = [\tilde{\mathbf{H}}(\boldsymbol{\theta})^2]_{00}$, it follows that

$$ \tilde{H}_{01}\tilde{H}_{10} > \tilde{H}_{01}(\boldsymbol{\theta})\tilde{H}_{10}(\boldsymbol{\theta}), \quad (20) $$

except if for all $j > 1$, $\tilde{H}_{j0}(\boldsymbol{\theta}) = 0$. In that latter case, $\tilde{H}_{01}\tilde{H}_{10} = \tilde{H}_{01}(\boldsymbol{\theta})\tilde{H}_{10}(\boldsymbol{\theta})$. Consequently, the saddlepoint of minimizing $E_I(\boldsymbol{\theta})$ corresponds to the case where for all $j > 1$, $\tilde{H}_{j0}(\boldsymbol{\theta}^*) = \tilde{H}_{j0} = 0$ [93].

A first important remark deals with the fact that cost functions, $C(\boldsymbol{\theta}) = E_G(\boldsymbol{\theta})$ or $E_I(\boldsymbol{\theta})$, cancel out if $|\tilde{\Phi}_1(\boldsymbol{\theta})\rangle$ and $|\Phi_0\rangle$ do not belong to the same symmetry-defined subspace of the Hamiltonian, i.e., if

$$ \tilde{H}_{01}(\boldsymbol{\theta}) = \langle\Phi_0|\mathbf{H}|\tilde{\Phi}_1(\boldsymbol{\theta})\rangle = 0. $$

In this way, the desired symmetries are imposed by construction in the guess vector $|\Phi_0\rangle$. This property is expected to drastically simplify the energy landscape since only the states belonging to the same symmetry-subspace as $|\Phi_0\rangle$ can contribute to $C(\boldsymbol{\theta})$ while states that do not belong to the targeted symmetry subspace have a zero contribution. Consequently, it is naturally expected that for the optimal parameters $\boldsymbol{\theta}^*$ minimizing $C$, the state $|\tilde{\Phi}_1(\boldsymbol{\theta}^*)\rangle$ exhibits the same quantum numbers (e.g., $\langle N \rangle$, $\langle S^2 \rangle$, $\langle S_z \rangle$) as $|\Phi_0\rangle$, regardless of the ansatz used to construct $\mathbf{R}(\boldsymbol{\theta})$. The use of symmetry-preserving cost functions $C(\boldsymbol{\theta})$ is therefore highly compatible with HEA.
A second important remark relates to the fact that, in contrast to the VQE algorithm which relies on the minimization of energy, minimizing both $C(\theta)$ does not guarantee to reach the ground state. Indeed, $E_G(\theta)$ remains negative only for positive $\Delta(\theta)$. Consequently, minimizing $E_G(\theta)$ decreases the energy but constrains $\omega(\theta) = |\langle\Phi_0|\Psi_0(\theta)\rangle|^2$ to be larger than 0.5 to ensure that $\Delta(\theta)$ is positive. It follows that if $|\Phi_0\rangle$ is not a good guess, minimizing $E_G(\theta)$ does not lead to the ground state but rather to the lowest-energy state within the symmetry-subspace where $\omega(\theta) > 0.5$. Similarly, using $E_I(\theta)$ as a cost function does not guarantee convergence to the ground state but instead promotes maximal interaction between $|\Phi_0\rangle$ and $|\tilde{\Phi}_1(\theta)\rangle$ within the same symmetry-subspace. To overcome this limitation, an iterative scheme is proposed in the following.

## B. Iterative and variational quantum eigensolver (IVQE)

The good-guess condition is replaced by the less constraining requirement for $|\Phi_0\rangle$ to be a binding state that has a non-zero overlap with the targeted ground state. To obtain the ground state, we propose an iterative loop condition compatible with the cost function $C(\theta)$. More precisely, we set up an iterative process such that the ground state $|\Psi_0^{(n)}(\theta^{(n)*})\rangle$ of the two-level system (see Fig. 1(b)) obtained at the end of an iteration $n$ is used as the guess vector $|\Phi_0^{(n+1)}\rangle$ for the next iteration $n+1$. The algorithm is initialized at iteration $n=0$ with $|\Psi_0^{(n=0)}\rangle = |\Phi_0\rangle$. The process for $n > 0$ can be decomposed into the following steps:

1. Loop condition and guess vector definition: $|\Psi_0^{(n-1)}\rangle \rightarrow |\Phi_0^{(n)}\rangle$ becomes the new guess vector for iteration $n$.

2. Find optimal parameter $\theta^{(n)*}$ for the unitary $\mathbf{R}(\theta^{(n)})$. We optimize the cost function $C(\theta^{(n)})$. The classical optimization algorithm requires multiple evaluations of the cost function. Inner loop for cost function computation:

(a) Construct a quantum circuit (quantum computer) representing the variational trial state $|\tilde{\Phi}_1^{(n)}(\theta^{(n)})\rangle$.

(b) Measure (quantum computer) all the necessary elements of the overlap and Hamiltonian matrix and the Hamiltonian:

$$S_{pn} = \langle \tilde{\Phi}_{1}^{(p)}(\theta^{(p)*}) | \tilde{\Phi}_{1}^{(n)}(\theta^{(n)}) \rangle, \quad (21)$$

$$\mathcal{H}_{pn} = \langle \tilde{\Phi}_{1}^{(p)}(\theta^{(p)*}) | \mathbf{H} | \tilde{\Phi}_{1}^{(n)}(\theta^{(n)}) \rangle \quad (22)$$

for $0 \le p \le n$ and considering $|\tilde{\Phi}_1^{(0)}(\theta^{(0)})\rangle = |\Phi_0\rangle$.

(c) Orthogonalize (classical computer) the variational trial state $|\tilde{\Phi}_1^{(n)}\rangle$ with respect to the overlap matrix $\mathbf{S}$. Update all elements of $\mathcal{H}$ in the orthonormalized basis set.

(d) Compute the cost function (classical computer) by constructing the matrix $\tilde{\mathbf{H}}^R(\theta^{(n)})$ from the elements of $\mathcal{H}$ in Eq. (22). Then, use Eq. (14) or Eq. (16) to obtain the cost function $C(\theta^{(n)})$.

The minimizing parameters are denoted by $\theta^{(n)*}$.

3. Two-level system ground state. (classical computer) We diagonalize $\tilde{\mathbf{H}}^R(\theta^{(n)*})$ in the subspace spanned by $|\Phi_0^{(n)}\rangle (= |\Psi_0^{(n-1)}\rangle)$ and $|\tilde{\Phi}_1^{(n)}(\theta^{(n)*})\rangle$, obtaining the optimal $\omega^{(n)*}$ [see Eq. (9)]. This leads to:

$$|\Psi_0^{(n)}\rangle = \omega^{(n)*} |\Psi_0^{(n-1)}\rangle \pm \sqrt{1 - (\omega^{(n)*})^2} |\tilde{\Phi}_1^{(n)}(\theta^{(n)*})\rangle. \quad (23)$$

Ground state properties at iteration $n$ such as the energy can be evaluated by using the diagonal form of $\tilde{\mathbf{H}}^R(\theta^{(n)*})$

$$E_0^{(n)} = \langle \Psi_0^{(n)} | \tilde{\mathbf{H}}^R(\theta^{(n)*}) | \Psi_0^{(n)} \rangle.$$

4. Check convergence. (classical computer) If $(\omega^{(n)*})^2 > 1 - \epsilon_w$ and/or the minimized cost function $C(\theta^{(n)*}) > -\epsilon_e$, the process is considered converged. Here, $\epsilon_w$ and $\epsilon_e$ are positive convergence criterion parameters that must be user-defined. Otherwise, the algorithm restarts at step 1 with $n = n + 1$.

This procedure is depicted in Fig. 2. The convergence criteria are based on the fact that if $\omega^{(n)*} \neq 1$ or if the cost function is strictly negative $C(\theta^{(n)*}) \neq 0$, then $|\Psi_0^{(n-1)}\rangle$ is not an eigenstate. On the contrary, if $\omega^{(n)*} = 1$ or $C(\theta^{(n)*}) = 0$, then $|\Psi_0^{(n-1)}\rangle$ is an eigenstate in the same symmetry subspace as the initial guess $|\Phi_0\rangle$. Note that only $|\tilde{\Phi}_1^{(n)}(\theta^{(n)})\rangle$ are constructed on the QC, while $|\Psi_0^{(n)}\rangle$ are never constructed explicitly but retrieved using $\tilde{\mathbf{H}}^R(\theta^{(n)*})$.

This algorithm will be referred to as the Iterative Variational Quantum Eigensolver (IVQE) algorithm. It falls into the class of VQEs, as each iteration is a new variational process with an updated guess-vector, for which the cost function $C(\theta^{(n)})$ is minimized. Both the gain and interaction energy can be used as the cost function. The former is expected to speed up the convergence, while the latter is expected to provide more stable numerical optimization.
<figure><img src="image_3.png" /><figcaption>FIG. 2: Schematic depiction of the iterative algorithm (top panel). Only the states $\{|\tilde{\Phi}_1^{(n)}(\theta^{(n)})\rangle\}$ in blue are optimized on a quantum computer to minimize the cost function $C(\theta^{(n)})$.</figcaption></figure>

### C. Variational quantum subspace method (VQSM)

We can exploit the fact that at each iteration $n$, the set of optimal trial states $\{|\tilde{\Phi}_1^{(p)*}(\theta^{(p)*})\rangle\}$ optimized iteratively in $0 \le p \le n$ iterations, constitutes a reduced orthonormal basis in which the Hamiltonian is expanded. It follows that we can revise steps 2d and 3 of the algorithm presented in the previous section III B in the spirit of quantum subspace method. More precisely, instead of using the two-level system Hamiltonian $\tilde{\mathcal{H}}^{R(n)}(\theta^{(n)})$ at iteration $n$ to compute the cost function and define the iteration ground state, we could diagonalize the $(n+1) \times (n+1)$ Hamiltonian matrix $\mathcal{H}^{(n)}$ within the subspace $\{|\tilde{\Phi}_1^{(p)}(\theta^{(p)})\rangle\}$. In parallel, we also need to adapt the cost functions, Eq. (14) or Eq. (16), to the subspace expansion approach. The cost function Eq. (14) is not changed and consists of the energy gain computed in the two-level system constituted from the ground state $|\Psi_0^{(n-1)}\rangle$ of $\mathcal{H}^{(n-1)}$ of the previous iteration and the variational trial vector $|\tilde{\Phi}_1^{(n)}(\theta^{(n)})\rangle$. It takes the same form as Eq. (14) but $\tilde{H}_{01}(\theta) = \langle\tilde{\Phi}_1^{(n)}(\theta^{(n)})|\mathbf{H}|\Psi_0^{(n-1)}\rangle$ and $\tilde{H}_{11}(\theta) = \langle\tilde{\Phi}_1^{(n)}(\theta^{(n)})|\mathbf{H}|\tilde{\Phi}_1^{(n)}(\theta^{(n)})\rangle$. In the same line, Eq. (16) can also be used to maximize the interaction energy with the ground state of the previous iteration using $\tilde{H}_{01}(\theta) = \langle\tilde{\Phi}_1^{(n)}(\theta^{(n)})|\mathbf{H}|\Psi_0^{(n-1)}\rangle$.

Interestingly, following property (20), we can define another type of cost function,

$$E_I'(\theta^{(n)}) = -\sqrt{\mathcal{H}_{n(n-1)}^{(n)}(\theta)\mathcal{H}_{(n-1)n}^{(n)}(\theta)} \quad (24)$$

which upon minimization provides a tridiagonal form of $\mathcal{H}$. Ultimately, given a fixed starting vector, it leads to a tridiagonalization of the Hamiltonian equivalent to the Householder, Givens, and Lanczos tridiagonalization methods on classical computers [76, 98–100], since Pick and Tomasek showed that given a starting vector the tridiagonal form is unique [101]. Consequently, the convergence performances should also be comparable. Consequently, we can expect from this analysis that the variational iterative algorithm we propose will converge

rapidly to the dominant eigenvalue, similar to classical Lanczos/Householder tridiagonalization scheme.

Altogether, it follows that steps 2d and 3 of the IVQE algorithm can be replaced by:

2. ...

(d') Compute the cost function: (classical computer) First diagonalize $\mathcal{H}^{(n)}$, then use the generalization of Eq. (14) or Eq. (16) to obtain the cost function $C(\theta^{(n)})$.

3'. Subspace expansion ground state: (classical computer) Once the parameters $\theta^{(n)*}$ optimized, we diagonalize $\mathcal{H}^{(n)}$ expanded in the orthonormal subspace spanned by $\{|\tilde{\Phi}_1^{(p)}(\theta^{(p)*})\rangle\}$, $0 \le p \le n$. The ground state of $\mathcal{H}^{(n)}$ is denoted $|\Psi_0^{(n)}\rangle$. We compute $\omega^{(n)*} = \langle\Psi_0^{(n-1)}|\Psi_0^{(n)}\rangle$ as the fidelity of the obtained ground state with respect to the ground state obtained at the previous iteration. Also the ground-state properties at iteration $n$ such as the energy can be evaluated as

$$E_0^{(n)} = \langle \Psi_0^{(n)} | \mathcal{H}^{(n)} | \Psi_0^{(n)} \rangle.$$

This algorithm will be denoted as the Variational Quantum Subspace Method (VQSM) algorithm. It falls in the class of quantum subspace methods, as in each iteration the subspace is enlarged, optimizing gain or interaction energy with the ground state of the previous iteration. In terms of performance, VQSM is expected to converge faster to the dominant eigenvalue than IVQE, at the costs of an additional classical diagonalization of the subspace Hamiltonian. Regarding the classical computational overhead, we acknowledge that, unlike IVQE, VQSM requires a matrix diagonalization step at each iteration. However, this overhead remains moderate in practice, as the size of the matrix to be diagonalized is determined by the number of iterations, which typically remains small due to the geometric convergence behavior of such methods. Importantly, this subspace dimension usually grows only polynomially with system size. Thus, while VQSM introduces more classical post-processing, it compensates by offering significantly improved accuracy and convergence. In this sense, VQSM sits at the intersection of VQE and quantum subspace methods, combining the advantages of both approaches: it leverages very shallow quantum circuits enabled by symmetry-preserving cost functions typical of HEA-based VQE, while also inheriting the robustness, noise resilience, and fast convergence associated with QSM. For large-scale applications, a hybrid strategy combining IVQE in the early iterations and switching to VQSM in later stages could effectively balance classical overhead and algorithmic accuracy.
<figure><img src="image_4.png" /><figcaption>FIG. 3: Schematic representation of the HEA used in this work. The circuit inside the square brackets representing one layer $l$ is repeated $N_L$ times.</figcaption></figure>

## D. Practical implementation and measurements

**Construction of the Trial Variational State on the Quantum Computer.** Ideally, the trial state $|\tilde{\Phi}_1(\theta)\rangle$ is prepared by applying a reflection $\mathbf{R}(\theta)$ to an “easy-to-prepare” initial state, such as the Hartree–Fock (HF) state – similar to the procedure used in standard Householder tridiagonalization. Using $E_I$ is expected to mitigate optimization difficulties since it is linear in $\mathbf{R}$. Indeed, since $\mathbf{R}$ can be defined up to a phase factor, the quantity $\tilde{H}_{01}(\mathbf{R}) = \langle\Phi_0|\mathbf{HR}|\Phi_1\rangle$, where $|\Phi_0\rangle$ and $|\Phi_1\rangle$ are fixed vectors, can be arbitrarily chosen to be positive. By contrast, $E_G$ is generally non-convex (quadratic in $\mathbf{R}$), although it tends toward convexity close to the non-interacting limit, suggesting that severe optimization challenges may be avoided in that regime. However, despite recent progress, implementing a reflection on a quantum computer remains technically challenging, as it typically requires multi-controlled NOT gates [102, 103]. That said, Householder reflections are unique among many other unitary transformations capable of block-diagonalizing $\mathbf{H}$ while preserving $|\Phi_0\rangle$. Therefore, rather than constraining the transformation to be a reflection, we construct the trial state $|\tilde{\Phi}_1(\theta)\rangle$ using a variational unitary quantum circuit applied to a simple reference state. Optimally, the variational circuit would exactly reproduce the effect of the targeted reflection. While the cost functions $E_I(E_G)$ introduced above shows convex properties (or nearly convex properties), respectively, the discrepancy between the actual circuit and the ideal Householder reflection likely results in a highly non-convex function, thus complicating the optimization. As a consequence, replacing the explicit construction of $\mathbf{R}(\theta)$ with a variational unitary circuit may introduce the typical issues encountered in variational methods – such as barren plateaus and the proliferation of local minima – well documented in other ansatz-based approaches [45].

**Hardware efficient Ansatz.** The quantum circuit in Fig. 3 aims to construct $|\tilde{\Phi}_1(\theta)\rangle = \mathbf{U}_{qc}(\theta)|\Phi_I\rangle$ and is constituted by an initialization, i.e the ‘easy’ construction of $|\Phi_I\rangle$ followed by the unitary circuit $\mathbf{U}_{qc}(\theta)$. For the initialization step, one might choose an excited Slater determinant $|\Phi_I\rangle$ that has an overlap with the ground state, i.e. that belong to the same symmetry subspace than $|\Phi_0\rangle$. However to remain general, here we just initialize

the circuit with the zero state. To construct the unitary part $\mathbf{U}_{qc}(\theta)$, we exploit the fact that the cost functions $E_G(\theta)$ and $E_I(\theta)$ enforce $|\tilde{\Phi}_1(\theta)\rangle$ to be in the same symmetry-subspace as the initial guess vector $|\Phi_0\rangle$. This allows in principle the use of any ansatz, including HEA to reduce circuit depth. As a sake of simplicity and controllability, we construct $\mathbf{U}_{qc}(\theta) = \text{HEA}$ as a first layer of $R_y(\theta)$ gates on each qubit followed by series of $N_L$ layers, each containing two-qubit entangling gates and single-qubit rotational gates. Note if the circuit is initialized with an excited state, the first $R_y(\theta)$ layer becomes unnecessary. For instance, a simple case depicted in Fig. 3 consists in the 1D $R_y$CNOT HEA [40] that we propose to use in the following as a proof of concept. The ansatz consists in $N_L$ layers containing an entangling operator. The entangling operator consists of a linear chain of $N_{\text{ent}} = (N_{\text{qubit}} - 1)$ two-qubit CNOT gates, employing a linear entanglement strategy. Each layer also includes a set of $R_y(\theta)$ rotation gates applied to all qubits. The total number of variational parameters is $(N_L+1)N_{\text{qubit}}$, while the number of CNOT gates is $N_L(N_{\text{qubit}}-1)$. Since qubit connectivity affects the number of entangling gates per layer, $N_{\text{ent}}$ can be adjusted accordingly, becoming $N_{\text{ent}} = N_{\text{qubit}}(N_{\text{qubit}}-1)/2$ for a fully connected setup.

**Measurement.** Using $E_G(\theta)$, $E_I(\theta)$, or $E_I'(\theta)$ as the objective function requires additional measurements compared to standard VQE. Specifically, it necessitates measuring the elements of $\mathcal{H}$ and $\mathcal{S}$. The diagonal elements can be obtained through standard projective measurements. The measurement of off-diagonal matrix elements has received particular attention in the context of Quantum Subspace Methods [11], as these quantities are essential for constructing effective Hamiltonians in a reduced basis. A standard approach to estimating such matrix elements of the form $O_{\alpha\beta} = \langle v_0|\hat{U}_\alpha^\dagger\hat{O}\hat{U}_\beta|v_0\rangle$, where $\hat{U}_\alpha$, $\hat{U}_\beta$ are unitary operators and $\hat{O}$ is typically the Hamiltonian or the identity, is the Hadamard test. This method requires an ancilla qubit and the implementation of controlled versions of $\hat{U}_\alpha$ and $\hat{U}_\beta$, which significantly increases the circuit depth and the number of entangling gates. Specifically, each single-qubit gate in $\hat{U}$ must be replaced by its controlled counterpart (requiring two CNOT gates), and each CNOT gate must be promoted to a Toffoli gate, which typically decomposes into six single-qubit gates and four or six CNOT gates. To mitigate these costs, recent alternatives have been developed that avoid the Hadamard test. Notably, Cortés and Gray [86] proposed a measurement-based protocol that reconstructs the off-diagonal matrix element $O_{\alpha\beta}$ from fidelity measurements, using overlaps between specially prepared quantum states. Alternatively, simplifications are possible by exploiting physical symmetries. For instance, if spin conservation holds, one can avoid implementing costly controlled time-evolution operators by instead performing a controlled initialization of the reference state [92]. Consequently, while the use of the Hadamard test is expected to significantly increase the circuit depth and CNOT count, these alternative strategies
gies allow for substantial resource savings and make the method more suitable for near-term applications.

*Circuit depth.* To contextualize the performance of our approach, we compare the CNOT gate counts required to achieve chemical accuracy with those reported for state-of-the-art adaptive methods such as ADAPT-VQE. For example, the original fermionic ADAPT-VQE requires approximately 2208 CNOT gates for $H_4$ and 28632 for $H_6$, while optimized variants reduce these numbers significantly to 812 CNOT gates for $H_6$ [26, 28, 38]. In contrast, our one or two layers HEA circuit corresponds to only 108-208 CNOT gates for $H_4$ and 168-324 CNOT gates for $H_6$, depending on the number of HEA layers, respectively, and using the Hadamard test to measure of diagonal elements of $\mathcal{H}$ and $\mathcal{S}$.

## IV. RESULTS

In this section, we evaluate the performance of IVQE and VQSM algorithms using different proposed cost functions, with the $H_n$ linear chain as a test case. Note that for $H_2$, chemical accuracy is achieved after only one iteration, regardless of the cost function used or the bond length. We employ the STO-3G basis set and the Jordan-Wigner transformation, which requires $2n$ qubits for $H_n$, along with an additional ancilla qubit when using the Hadamard test.

### A. Convergence with respect to the HEA circuit depth

We first propose to focus on the results obtained at the first iteration, where there are no differences between the IVQE and VQSM approaches. Indeed, at this stage, the constructed subspace of VQSM matches the two-level system of IVQE. This allows us to carefully investigate the influence of the HEA circuit depth and to demonstrate that short circuits can be used. In parallel, we study the influence of choosing $E_G$ or $E_I$ in particular by showing that the topology of the energy landscape – and thus the ease of performing classical optimization – depends drastically on the choice of cost function.

In Fig. 4 we present the ground-state energy of $H_4$ chains as a function of the interatomic distance $d_{H-H}$ together with the error $\Delta E = E(\theta^*) - E^{\text{FCI}}$ with respect with the corresponding full configuration interaction (FCI) energy as a function of the number of layers $N_L$. Results for $E_G(\theta)$ (panels (a) and (c)) are compared with those of $E_I(\theta)$ (panels (b) and (d)). The optimizations are performed using the COBYLA method, either until convergence is reached or until the process is arbitrarily stopped after 150,000 cost function evaluations, a setup consistently used for all results presented throughout the manuscript. It shows that using the $C(\theta) = E_G(\theta)$, for small interatomic distances (i.e., $d_{H-H} \le 1.9$ Å), the energy closely matches the FCI energy

when $N_L = 5$ and that $\Delta E$ systematically decreases when increasing the circuit depth. More precisely, in that case the chemical accuracy can be achieved using only a few layers (e.g., six layers for $d_{H-H} = 1$, as seen in (b). However, when increasing $d_{H-H}$ beyond 1.9 Å, the ground-state energy appears to converge faster but saturates at a value significantly higher than the FCI energy. As discussed in Sec. III A, this behavior results from the violation of the *good-guess* assumption – that is, the HF weight $\omega^2$ in the ground state is lower than $1/2$. As shown in Appendix A, $\omega^2(\theta^*)$ quickly saturates at a value of $1/2$ after only a few layers for $d_{H-H} \ge 2$ Å. This drawback is overcome by starting the iterative scheme, as will be shown later on. When using $C(\theta) = E_I(\theta)$ an accurate description of the ground state is not expected in the first iteration since it optimizes the interaction energy with the Hartree-Fock state. It follows that in this case significant deviation with the FCI energy is expected. In particular, $\Delta E$ is always higher than 10 mHartree, regardless of the interatomic distance (see panel (d)). Interestingly, the errors on the energy converge rapidly and saturate within only a few HEA layers, suggesting that rather shallow circuits can be used. Note that as expected, using $E_G(\theta)$ or $E_I(\theta)$ leads to the same saddle point for large $d_{H-H}$, as suggested by Eq. (17) since at large $d_{H-H}$, $E_G \sim \sqrt{\tilde{H}_{01}\tilde{H}_{10}/\Delta} \propto E_I$.

Regarding the convergence of the cost function, we show in Fig. 5 the difference between the cost function $C(\theta^*)$ ($C = E_G$ or $E_I$) obtained within $N_L$ layers of HEA and its optimal value $C^{\text{min}}$ obtained classically by “brut-force” minimization of $C(\mathbf{R})$ where $\mathbf{R} = 1-2|\mathbf{v}\rangle\langle\mathbf{v}|$ is a variational reflection with $|\mathbf{v}\rangle$ a normalized vector containing the variational parameters, of length 256 for $H_4$. The minimization of $C(\mathbf{R})$ is numerically efficient and converges rapidly for both cost functions, even though $E_G(\mathbf{R})$ is not linear in $\mathbf{R}$, unlike $E_I(\mathbf{R})$. However, when $\mathbf{R}$ is replaced by a parametrized HEA quantum circuit, the optimization is more tedious. To analyze more deeply the optimization properties and the topology of the energy landscape, we conduct a stochastic analysis of the optimization process by running multiple optimizations with different random initial parameters. The results are shown in Fig. 5(b) and (d), where a box plot illustrates the difference between the optimized cost function $C(\theta^*)$ and its exact optimal value $C^{\text{min}}$ for $d_{H-H} = 2$ Å with respect to different number of layers $N_L$. The whiskers, empty boxes, and filled boxes highlight the regions containing the 75%, 50%, and 10% lowest optimized cost function values $C(\theta^*)$, respectively. Overall, the results indicate that the energy landscape is far from convex, making it challenging to reach the global minimum – particularly as the circuit depth increases. This suggests that global optimization approaches such as Basin-hopping [30] are required to reach the global minimum, rather than local optimizers such as gradient-based. Notably, the convergence rate appears to be significantly accelerated when using $E_I(\theta)$ instead of $E_G(\theta)$. Consequently, the mapping of $\mathbf{R}$ into
<figure><img src="image_5.png" /><figcaption>FIG. 4: *Ground-state energy at the first iteration.* Panels (a) and (b): ground-state energy (in Hartree) as a function of the interatomic distance $d_{H-H}$ (in Angström) for different circuit depths corresponding to $N_L = 1, 3$, or 5 layers in the HEA using $E_G(\theta)$ and $E_I(\theta)$ cost functions, respectively. Results are compared with values obtained using the Hartree-Fock approximation and FCI method. Panels (c) and (d): absolute error $\Delta E = E(\theta^*) - E^{\text{FCI}}$ (in Hartree) as a function of the number of layers in the HEA for selected interatomic distances of $d_{H-H} = 1.0, 2.0, 3.0$, and $4.0$ Å using $E_G(\theta)$ and $E_I(\theta)$ cost functions, respectively.</figcaption></figure>

a unitary parametrized quantum circuit, especially using HEA-type circuits, is itself non-convex and exhibits parameter redundancy. This makes the optimization of $E_I(\theta)$ still a complex task, particularly as the number of layers and variational parameters increase. Since $E_I(\theta)$ is expected to be numerically more stable and less prone to optimization issues than $E_G(\theta)$, we focus on it for the remainder of the manuscript. In the following, we address both the optimization challenges and the noise inherent in deep HEA circuits by trading circuit depth for a higher number of iterations. We show that increasing the number of iterations can yield better results even with shallow circuits composed of only one or two layers.

In summary, while optimization challenges are inherent to variational algorithms, the use of shallow circuits and a symmetry-preserving cost function is believed to enhances convergence and stability. Nonetheless, we acknowledge the limitations of hardware-efficient ansätze and anticipate that incorporating more expressive, problem-tailored circuits or adaptive strategies could further improve optimization performance in future work. Moreover, the iterative framework is expected to mitigate their impact by transforming a difficult global optimization into a sequence of smaller, more manageable subproblems.

## B. Trading number of layers for number of iterations

Fig. 6 presents the results of the IVQE and VQSM algorithms for computing the ground-state properties of the $H_4$ linear chain as a function of the interatomic distance $d_{H-H}$. The study uses a single-layer HEA and compares the results to the HF approximation and FCI solutions. The top panels (a) and (c) shows the ground-state energy (in Hartree) as a function of interatomic distance for different iteration numbers (1, 2, 3, 5, 10), using the $E_I(\theta)$ cost function. Initially, for small iteration counts, both IVQE and VQSM methods deviate significantly from the FCI reference, especially near the dissociation regime. However, as the number of iterations increases, the results get closer to the exact energy. From a number of iterations $n \ge 10$, the curve aligns well with the FCI solution, demonstrating that both IVQE and VQSM refine their estimates effectively. Panels (b) and (d) quantify the error on the energy $\Delta E$ of IVQE and VQSM, respectively, plotting the deviation from the FCI results. $\Delta E$ consistently decreases as the number of iterations $n$ increases, indicating that both methods converge monotonically with respect to the number of iterations. The improvement in accuracy is more significant at small and large bond lengths, with errors falling below chemical accuracy (1.6 mHartree) after a few iterations. This suggests that both methods efficiently cap-
<figure><img src="image_6.png" /><figcaption>FIG. 5: Convergence with respect to the circuit depth and optimization issues at first iteration. Panels (a) and (c): Difference between the cost function $C(\boldsymbol{\theta}^*)$ (in Hartree) with the optimal value $C^{\min}$, as a function of the number of layers in the HEA for selected interatomic distances of $d_{\text{H-H}} = 1.0, 2.0, 3.0$, and $4.0$ Å, using $E_G(\boldsymbol{\theta})$ and $E_I(\boldsymbol{\theta})$ cost functions, respectively. Panels (b) and (d): Difference between the cost function $C(\boldsymbol{\theta}^*)$ (in Hartree) with the optimal value $C^{\min}$, as a function of the number of layers in the HEA with box-plot representing 25% (colored box), 50% (empty box) and 75% (whiskers) of optimization outcomes for the cost function when starting with different random variational parameters. Results are given as a function of the number of layers in the HEA for interatomic distances of $d_{\text{H-H}} = 2.0$ Å, using $E_G(\boldsymbol{\theta})$ (b) and $E_I(\boldsymbol{\theta})$ (d) cost functions.</figcaption></figure>

ture electronic structure and correlation effects, even at large distances. However, residual errors larger than the chemical accuracy is persisting at intermediate regime $d_{\text{H-H}} \sim 2.0-3.0$ Å using IVQE, even after 15 iterations. A key difference between IVQE and VQSM emerges in their convergence behavior and efficiency. Both methods systematically improve with more iterations, but VQSM converges faster, especially for larger bond lengths. For instance, chemical accuracy is reached after $n = 6, 9, 15$, and $5$ iterations for $d_{\text{H-H}} = 1.0, 2.0, 3.0$ and $4.0$ Å, respectively, for IVQE; compared to $n = 5, 7, 10$, and $3$ iterations for VQSM. Therefore, the subspace expansion in VQSM provides a more efficient algorithm, accelerating convergence to the exact solution, at the expense of larger subspace Hamiltonian $\mathcal{H}^{(n)}$ to be diagonalized classically. Overall VQSM demonstrates superior convergence properties, highlighting its potential as a more efficient variational algorithm for electronic structure calculations. Accordingly, we focus on VQSM associated with $E_I(\boldsymbol{\theta})$ as a cost function in the rest of the manuscript.

### C. Convergence rate of VQSM versus the standard Lanczos algorithm

In Fig. 7, we analyze the convergence rate of the VQSM algorithm for the linear $H_4$ molecule and compare it with the standard (classical) Lanczos algorithm. As previously discussed in the manuscript, both VQSM using $E_I'(\boldsymbol{\theta})$ and the Lanczos algorithm are expected to produce the same tridiagonal matrix within the constructed subspace. Consequently, we also expect a geometric convergence for VQSM characterized by an exponential decay of the error: $\epsilon^{(n)} = \epsilon^0 r^n$, where $\epsilon^{(n)}$ denotes the error at iteration $n$, $\epsilon^0$ is a constant, and $r$ is the convergence rate. Panel (a) shows the energy error $\Delta E$ obtained using a single-layer HEA as a function of the number of iterations $n$, for different interatomic distances $d_{\text{H-H}}$. For small bond lengths ($d_{\text{H-H}} \le 2$ Å), $\Delta E$ decreases almost linearly (on a logarithmic scale) with the iteration number exhibiting the expected geometric convergence. The logarithm of the convergence rate $\ln(r)$ can be directly extracted as the slope of the error vs iteration curve when plotted on a logarithmic scale. For $d_{\text{H-H}} \le 2$ Å, we obtain $r \sim 0.5$ using a linear regression, meaning that at each iteration, the error is divided by a factor of two. However, for $d_{\text{H-H}} > 2$ Å, the error initially decreases rapidly in the first few iterations but then slows significantly, leading to
<figure><img src="image_7.png" /><figcaption>FIG. 6: Single layer HEA for linear H₄ molecule IVQE vs VQSM. Ground-state energy (a) and corresponding error ΔE (b), in Hartree, as a function of the interatomic distance d<sub>H-H</sub> (in Ångström), using the cost function E<sub>I</sub>(θ) and within the IVQE algorithm. In panel (a), results are shown at iterations n = 1, 2, 3, 5, and 10. In panel (b), results are shown for n = 1–15, with color shading ranging from dark (early iterations) to light (later iterations), as illustrated for selected iterations in the legend. Results are compared to the Hartree–Fock approximation and exact diagonalization (FCI). The horizontal green dot-dashed line in panel (b) indicates the chemical accuracy threshold (1.6 mHartree). Panels (c) and (d) show the same quantities as (a) and (b), respectively, but for the VQSM algorithm.</figcaption></figure>

a “plateau-like” behavior.

The efficiency or convergence rate of the VQSM algo-
rithm is expected to depend on several factors. On one
hand, as other Krylov-based algorithm, it is influenced by
physical parameters, such as the first eigengap, i.e. the
energy gap between the ground and first excited states
within the symmetry-constrained Hilbert space. On the
other hand, numerical parameters also play a role, in-
cluding the circuit depth (which determines the ability
to span the Fock space) and optimization success (i.e.,
the ability to reach the global minimum in a complex
energy landscape). To investigate the effect of circuit
depth and optimization issue on VQSM convergence, we
present in Fig. 7(b) and (c) the ground-state energy er-
ror ΔE as a function of the iteration number for HEA
circuits with different layer counts, compared to classi-
cal algorithms. We focus on two representative cases:
d<sub>H-H</sub> = 1.0 Å (b) and d<sub>H-H</sub> = 3.0 Å (c). In parallel, Fig. 8
illustrates the convergence of the ground-state fidelity
with respect to the FCI solution. Panels (a) and (b)
show fidelity versus iteration number for different HEA
depths, again for d<sub>H-H</sub> = 1.0 Å and d<sub>H-H</sub> = 3.0 Å, re-
spectively. The results in Figs. 7(b), 7(c), 8(a) and 8(b)
are compared with the classical Lanczos algorithm (black
curves). Dark gray curves are obtained by directly mini-
mizing the cost function E'<sub>I</sub>(R) in Eq. (24) with respect

to a variational reflection **R** = 1 − 2|**v**⟩⟨**v**|, where |**v**⟩ is
a normalized variational vector of the size of the many-
body basis set. Similarly, light gray curves correspond
to results obtained by minimizing *E*<sub>*I*</sub>(*R*). Interestingly,
switching from *E*'<sub>*I*</sub>(*θ*) to *E*<sub>*I*</sub>(*θ*) has only a minor impact,
meaning that both cost function can also lead to the same
triangular form of the Hamiltonian in the reduced sub-
space. As expected, minimizing *E*'<sub>*I*</sub>(*θ*) reproduces quan-
titatively, up to numerical setting, the Lanczos results.
First of all, we remark that the “plateau-like” behavior,
observed in Fig. 7(a) close to dissociation, is also present
when using classical methods such as Lanczos. This sug-
gests that the slower convergence near dissociation has
a physical origin rather than resulting from optimiza-
tion issues. When turning to quantum implementation,
replacing the reflection operator by a parameterized uni-
tary HEA circuit (with *N*<sub>*L*</sub> = 1, 2, 3 layers) to minimize
*E*<sub>*I*</sub>(*θ*) reduces the convergence rate, since the resulting
“Krylov vector” is only approximated. For instance, at
*d*<sub>H-H</sub> = 1.0 Å, the convergence rate for the energy error
drops from *r* ~ 0.50 for a single-layer HEA to *r* ~ 0.47
and *r* ~ 0.40 for two- and three-layer HEAs, respectively.
Even though the number of layers is increased, it remains
far from the convergence rate of the Lanczos algorithm
in this case, which shows *r* = 0.07. Note that the conver-
gence trend for fidelity closely follows that of the energy
<figure><img src="image_8.png" /><figcaption>FIG. 7: Convergence rate of VQSM for the linear H$_4$ molecule: (a) Energy error ΔE (in Hartree) as a function of the number of layers n, for various bond lengths d$_{H-H}$ = 1.0, 2.0, 3.0, and 4.0 Å (see legend), using the cost function E$_I$ (θ) and a single-layer HEA. (b) and (c): same as (a) for d$_{H-H}$ = 1.0 (b) and d$_{H-H}$ = 3.0 (c). In (b) and (c), results are shown for HEAs with 1, 2, and 3 layers (colored curves with symbols, see legend) and are compared with classical methods shown using shade of gray. In particular, results obtained using the standard Lanczos algorithm are shown in black, while results from classical minimization of E'$_I$ (R) and E$_I$ (R) are shown in dark gray and light gray, respectively.</figcaption></figure>

error.

Finally, Fig. 8(c) displays the fidelity as a function of d$_{H-H}$ at various iteration steps n, which can be compared to powers of the ratio (E$_1$/E$_0$)$^n$ in Fig. 8(d), where E$_0$ and E$_1$ are the ground- and first excited-state energies within the same symmetry subspace. Indeed, for the power method, r is known to be proportional to the ratio E$_1$/E$_0$. While the convergence rate of the Lanczos algorithm is more precisely governed by the ratio between the first eigengap E$_1$ − E$_0$ and the diameter of the remain-

ing spectrum [99], the ratio E$_1$/E$_0$ already provides a rough estimate of the convergence rate. In both cases, for nearly degenerate systems where E$_1$/E$_0$ → 1 (i.e., E$_1$ − E$_0$ → 0), as encountered near dissociation (see the inset of Fig. 8(d)), the convergence becomes slow and can exhibit plateau-like behavior. In contrast, for systems with a large energy gap (E$_1$/E$_0$ → 0 or E$_1$ − E$_0$ ≫ 1), i.e., for shorter d$_{H-H}$, convergence is expected to be fast. This analysis explains the plateau behavior observed in both VQSM and Lanczos algorithm in Figs. 7(c) and 8(b), and confirms that it arises from physical limitations rather than from optimization artifacts.

Altogether, increasing the number of HEA layers N$_L$ improves the convergence rate, although a significantly deep circuit would be required to fully match the performance of the Lanczos algorithm and probably also associated with optimization issues. Nonetheless, a good trade-off can be achieved between (i) the number of iterations, (ii) the circuit depth, and (iii) the occurrence of optimization difficulties to reach a desired level of accuracy. This balance is especially crucial in near-degenerate regimes, where improving circuit expressibility becomes necessary but increasingly costly. In practice, choosing the minimal depth that ensures sufficient convergence while maintaining feasible optimization and coherence constraints on quantum hardware is key for the efficient implementation of VQSM on near-term devices. Finally, as a sake of completeness we show in Fig. 9 results for H$_n$ chains 3 ≤ n ≤ 6. It shows that as the system size increases, the convergence rate slows down slightly, consistent with a decreasing energy gap and the limited ansatz expressivity. However, the geometric convergence is preserved, demonstrating the robustness of the VQSM as variational Krylov-based approach across different system sizes.

## D. Beyond Ground-State Properties

Beyond the accurate estimation of ground-state energies, the presented IVQE and VQSM algorithm enables access to low-lying excited-state properties. In this section, we demonstrate the capability of VQSM to compute charge and spin gaps, which are essential descriptors of electronic excitation spectra and chemical reactivity. In the following we focus on VQSM using cost function E$_I$ (θ) with a single HEA layer, since this has been shown to converge rapidly without optimization issues.

Fig. 10 shows the evaluation of the charge gap, defined as

$$ \Delta E_c = E(\mathbf{H}_4^+) + E(\mathbf{H}_4^-) - 2E(\mathbf{H}_4), \quad (25) $$

as a function of d$_{H-H}$. From a practical perspective, the energies of the charged states E(H$_4^+$) and E(H$_4^-$) are obtained independently from the ground-state energy of E(H$_4$) by constructing subspaces initialized from the corresponding charged Hartree-Fock states.

At low iteration numbers, the estimated charge gaps significantly deviate from the exact (FCI) results, espe-
<figure><img src="image_9.png" /><figcaption>FIG. 8: Convergence rate of the Fidelity using VQSM for the linear $H_4$ chain with $E_I(\theta)$ as the cost function.</figcaption></figure>

Panels (a) and (b): Fidelity as a function of the iteration number $n$ at bond lengths $d_{H-H} = 1.0$ Å and $3.0$ Å, respectively. The colored lines refer to VQSM results obtained for different HEA circuit depths having $N_L = 1, 2$ and $3$ layers, distinguished by the different symbols as specified in the legend. They are compared with values obtained using Lanczos method (black curve) or by brute-force minimization of $E_I'(\mathbf{R})$ (dark gray) and $E_I(\mathbf{R})$ (light gray) with respect to the reflection matrix $\mathbf{R}$. Panel (c): Fidelity as a function of the bond length $d_{H-H}$ (in Å) at different iterations ($n = 1-15$) with color shading ranging from dark (early iterations) to light (later iterations), as illustrated for selected iterations in the legend. Panel (d): Powers of the ratio between the ground-state energy $E_0$ and the first excited energy belonging to the same symmetry subspace $E_1$ as a function of $d_{H-H}$ for several value of the iteration $n$. The inset highlights the case for $n = 1$.

<figure><img src="image_10.png" /><figcaption>FIG. 9: Geometrical convergence of VQSM on hydrogen chains $H_n$ with $n = 3$ to 6 using a fixed one-layer HEA ansatz and $E_I(\theta)$ Cost Function, for $n = 3$ and $5$ we have considered $H_n^-$. For a sake of clarity the relative error is renormalized by the error at the first iteration.</figcaption></figure>

In the dissociation regime ($d_{H-H} > 2.5$ Å), where electron correlation becomes dominant. This underestimation reflects the limited expressiveness of the Krylov subspace when only a small number of vectors is used. However, as the number of iterations increases, the method systematically improves, and the computed $\Delta E_c$

<figure><img src="image_11.png" /><figcaption>FIG. 10: VQSM Band Gap with $E_I(\theta)$ Cost Function and a Single Layer HEA. Charge gap</figcaption></figure>

$\Delta E_c = E(H_4^+) + E(H_4^-) - 2E(H_4)$ as a function of the interatomic distance $d_{H-H}$ (in Å) at different iterations ($n = 1, 2, 3, 5$, and $10$). Results are compared with values obtained using the Hartree-Fock approximation and exact diagonalization (FCI).

values converge towards the FCI results. Already at iteration $n = 5$, the method captures the correct qualitative trend with respect to bond length, and by iteration
<figure><img src="image_12.png" /><figcaption>FIG. 11: VQSM Singlet-Triplet Gap with $E_I$ Cost Function and a Single Layer HEA. Singlet-Triplet gap $\Delta E_{ST} = E_T(H_4) + E_S(H_4)$ as a function of the interatomic distance $d_{H-H}$ (in Å) at different iterations $n$ for the open (a) and closed (b) H$_4$ molecules. Results are compared with values obtained using the Hartree-Fock approximation and FCI method.</figcaption></figure>

$n=10$, the charge gap is in good quantitative agreement with the reference data across the entire bond-length range. In contrast, the Hartree-Fock approximation fails to capture the correct asymptotic behavior at large distances, converging to zero in the dissociation limit.

Let us now explore the singlet-triplet gap,

$$ \Delta E_{ST} = E_T(H_4) - E_S(H_4), \quad (26) $$

for both linear-chain and square geometries of the neutral H$_4$ molecule. The gap is plotted in Fig. 11 as a function of $d_{H-H}$ for various iteration numbers ($n=1, 2, 3, 5$, and $10$ for the linear case, $n=7, 10$, and $15$ for the square H$_4$ molecule), and is compared with both FCI and Hartree-Fock results. In the linear geometry Fig. 11(a), the Hartree-Fock approximation incorrectly predicts a singlet-triplet gap that increases from negative to positive values as the bond length increases, suggesting a ground-state spin crossover. However, FCI data confirm that the singlet remains the ground state throughout, with the spin gap decreasing to zero near dissociation due to diminishing electron-electron repulsion. VQSM shows rapid convergence with respect to the number of iteration, capturing the correct trend of $\Delta E_{ST}$ already at iteration $n=3$, and being nearly indistinguishable from the FCI reference over the entire bond-length range at iteration $n=10$. For the square geometry Fig. 11(b), the singlet-triplet gap remains small (on the order of milli-

hartrees), and thus demands high numerical precision. Here too, VQSM exhibits excellent performance: at iteration $n=10$, the method reproduces both the curvature and the magnitude of the FCI reference. In this case, a triplet-to-singlet transition occurs near $d_{H-H} \sim 0.8$ Å. In the H$_4$ ring, all hydrogen atoms are closer together, enhancing the effects of Fock exchange, which stabilizes ferromagnetic interactions and thus a triplet configuration. The system thus undergoes a subtle competition between ferromagnetic exchange and electron delocalization across the ring, leading to a delicate spin gap highly sensitive to correlation effects.

## V. CONCLUSION

In this manuscript, we have proposed alternative symmetry-preserving cost functions, $E_G(\theta)$ and $E_I(\theta)$, for the VQE algorithm, aiming to reduce optimization challenges, particularly when employing HEA. Since the ground-state energy is not directly minimized, ground-state properties can still be obtained by activating an iterative procedure. Two iterative algorithms have thus been derived: IVQE and VQSM. In both approaches, at each iteration $n$, a variational trial vector $|\tilde{\Phi}_1^{(n)}(\theta)\rangle$ is optimized by minimizing the considered cost function. The main distinction between the two algorithms lies in the additional classical numerical tasks required by VQSM, where the set $\{|\tilde{\Phi}_1^{(p)}(\theta)\rangle\}$, for $0 \le p \le n$, forms a reduced Krylov-like subspace.

Both algorithms, along with the proposed cost functions $E_G(\theta)$ and $E_I(\theta)$, have been implemented and tested on the H$_4$ molecule, and compared against Full Configuration Interaction and Hartree-Fock reference results. Regarding the cost functions, both $E_G(\theta)$ and $E_I(\theta)$ perform similarly. However, $E_I(\theta)$ has been favored throughout the manuscript due to its expected greater numerical stability and reduced sensitivity to optimization issues, arising from its linear properties. In particular, it has been shown that using HEA as an ersatz reflection can lead to a highly non-convex energy landscape, especially as the number of layers increases.

In terms of algorithmic performance, VQSM outperforms IVQE at the cost of modest additional classical computation, which remains tractable. Furthermore, a formal connection between VQSM and the classical Lanczos algorithm has been established, as both methods ideally yield the same tridiagonal representation of the Hamiltonian in the reduced subspace. Remarkably, even with a shallow, single-layer ansatz, the VQSM algorithm exhibits geometric convergence. In practice, a significantly deeper circuit would be required to fully match the performance of the Lanczos algorithm, potentially at the expense of optimization stability. As VQSM can be seen as a variational construction of the Krylov subspace, we expect that a small number of iterations is needed to achieve a given accuracy. Moreover, the algorithm has
been shown to converge within a single iteration if the good-guess condition $\omega = 0.5$ is fulfilled.

Nonetheless, a good trade-off can be achieved between (i) the number of iterations, (ii) the circuit depth, and (iii) the complexity of the optimization process, to reach a desired level of accuracy. This balance is particularly crucial in near-degenerate regimes, where enhancing circuit expressibility becomes necessary but computationally costly. In practical applications, choosing the minimal circuit depth that ensures sufficient convergence while preserving feasible optimization and coherence constraints on quantum hardware is essential for the efficient deployment of VQSM on near-term quantum devices.

Finally, we have highlighted the strength of the VQSM approach in accurately and efficiently capturing both charge and spin excitations using a compact variational ansatz. The method demonstrates systematic and robust convergence with iteration number, achieving chemical accuracy for both neutral and charged states using a limited number of vectors in the reduced subspace. These results position VQSM as a versatile variational quantum algorithm, extending beyond ground-state energy estimation to address spectroscopic features and electronic correlations in quantum many-body systems.

# ACKNOWLEDGMENTS

The authors thank D. Lacroix from IJCLab in Saclay (France) for fruitful discussions. HAA is grateful to the Hubert Curien program (CAMPUS FRANCE) for the support as well as the Abdus Salam ICTP for the HPC facilities. Computer time for this study was also provided by the computing facilities of the MCIA (Mésocentre de Calcul Intensif Aquitain). BS and MS acknowledge the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-23-PETQ-0006.

# AUTHORS CONTRIBUTION

Algorithm implementation and data acquisition: H.A. and M.S.(equal), A. P. (support); Concept and methodology: M.S. (lead) and B.S. (support); Supervision and project management: M.S.; Writing original draft: M.S (lead) and B.S. (support). All authors edited and commented on the manuscript.

# Appendix A: Hartree-Fock weight

Fig. 12 shows that $\omega^2(\theta^*)$ saturates to a value of 1/2 after only a few layers for $d_{H-H} = 2, 3$, and $4$ Å. When $\omega^2 \le 1/2$, the initial state $|\Phi_0\rangle$ is no longer a good state, meaning it does not sufficiently overlap with the true ground state. As a result, the cost function $E_G(\theta)$ fails to converge within a single iteration, necessitating the use of an iterative procedure.

<footnote n="1">1 A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, *Science* **309**, 1704 (2005).</footnote>

<footnote n="2">2 B. Bauer, S. Bravyi, M. Motta, and G. K.-L. Chan, *Chem. Rev.* **120**, 12685 (2020).</footnote>

<footnote n="3">3 S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, *Rev. Mod. Phys.* **92**, 015003 (2020).</footnote>

<footnote n="4">4 Z. Zimborás, B. Koczor, Z. Holmes, E.-M. Borrelli, A. Gilyén, H.-Y. Huang, Z. Cai, A. Acín, L. Aolita, L. Banchi, F. G. S. L. Brandão, D. Cavalcanti, T. Cubitt, S. N. Filippov, G. García-Pérez, J. Goold, O. Kálmán, E. Kyoseva, M. A. C. Rossi, B. Sokolov, I. Tavernelli, and S. Maniscalco, "Myths around quantum computation before full fault tolerance: What no-go theorems rule out and what they don't," (2025), arxiv:2501.05694 [quant-ph].</footnote>

<footnote n="5">5 K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.-C. Kwek, and A. Aspuru-Guzik, *Rev. Mod. Phys.* **94**, 015004 (2022).</footnote>

<footnote n="6">6 P. J. J. O'Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, et al., *Phys. Rev. X* **6**, 031007 (2016).</footnote>

<footnote n="7">7 A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O'Brien, *Nat Commun* **5**, 4213 (2014).</footnote>

<footnote n="8">8 J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, *New J. Phys.* **18**, 023023 (2016).</footnote>

<footnote n="9">9 J. R. McClean, M. E. Kimchi-Schwartz, J. Carter, and W. A. de Jong, *Phys. Rev. A* **95**, 042308 (2017).</footnote>

<footnote n="10">10 N. H. Stair, R. Huang, and F. A. Evangelista, *J. Chem. Theory Comput.* **16**, 2236 (2020).</footnote>

<footnote n="11">11 M. Motta, W. Kirby, I. Liepuoniute, K. J. Sung, J. Cohn, A. Mezzacapo, K. Klymko, N. Nguyen, N. Yoshioka, and J. E. Rice, *Electron. Struct.* **6**, 013001 (2024).</footnote>

<footnote n="12">12 M.-H. Yung, J. Casanova, A. Mezzacapo, J. McClean, L. Lamata, A. Aspuru-Guzik, and E. Solano, *Scientific reports* **4**, 3589 (2014).</footnote>

<footnote n="13">13 P. K. Barkoutsos, J. F. Gonthier, I. Sokolov, N. Moll, G. Salis, A. Fuhrer, M. Ganzhorn, D. J. Egger, M. Troyer, A. Mezzacapo, S. Filipp, and I. Tavernelli, *Phys. Rev. A* **98**, 022322 (2018).</footnote>

<footnote n="14">14 J. Romero, R. Babbush, J. R. McClean, C. Hempel, P. J. Love, and A. Aspuru-Guzik, *Quantum Sci. Technol.* **4**, 014008 (2018).</footnote>

<footnote n="15">15 J. Lee, W. J. Huggins, M. Head-Gordon, and K. B. Whaley, *J. Chem. Theory Comput.* **15**, 311 (2019).</footnote>
<figure><img src="image_13.png" /><figcaption>FIG. 12: Weight $w^2$ at the first iteration and as a function of the H-H bond length using the $E_G$ cost function. Results are given considering different number of HEA layers $N_L$ and compared with the HF weight in the FCI ground state (black curve).</figcaption></figure>

16 F. A. Evangelista, G. K.-L. Chan, and G. E. Scuseria, *The Journal of Chemical Physics* **151**, 244112 (2019).

17 I. O. Sokolov, P. K. Barkoutsos, P. J. Ollitrault, D. Greenberg, J. Rice, M. Pistoia, and I. Tavernelli, *The Journal of Chemical Physics* **152**, 124107 (2020).

18 A. Khamoshi, G. P. Chen, F. A. Evangelista, and G. E. Scuseria, *Quantum Sci. Technol.* **8**, 015006 (2022).

19 S. E. Smart and D. A. Mazziotti, *Phys. Rev. Lett.* **126**, 070504 (2021).

20 Q. Marécat, B. Senjean, and M. Saubanère, *Phys. Rev. B* **107**, 155110 (2023).

21 D. Materia, L. Ratini, C. Angeli, and L. Guidoni, *J. Phys. Chem. A* (2024), 10.1021/acs.jpca.4c03756.

22 G. Greene-Diniz and D. Muñoz Ramo, *Int. J. Quantum Chem.* **121**, e26352 (2021).

23 B. O’Gorman, W. J. Huggins, E. G. Rieffel, and K. B. Whaley, “Generalized swap networks for near-term quantum computing,” (2019), arxiv:1905.05118 [physics, physics:quant-ph].

24 H. R. Grimsley, D. Claudino, S. E. Economou, E. Barnes, and N. J. Mayhall, *J. Chem. Theory Comput.* **16**, 1 (2019).

25 M. Motta, E. Ye, J. R. McClean, Z. Li, A. J. Minnich, R. Babbush, and G. K.-L. Chan, *npj Quantum Inf* **7**, 1 (2021).

26 H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, *Nat Commun* **10**, 3007 (2019).

27 J. Liu, L. Wan, Z. Li, and J. Yang, *J. Chem. Theory Comput.* **16**, 6904 (2020).

28 H. L. Tang, V. Shkolnikov, G. S. Barron, H. R. Grimsley, N. J. Mayhall, E. Barnes, and S. E. Economou, *PRX Quantum* **2**, 020310 (2021).

29 C. Feniou, M. Hassan, D. Traoré, E. Giner, Y. Maday, and J.-P. Piquemal, *Commun Phys* **6**, 1 (2023).

30 H. G. A. Burton, D. Marti-Dafcik, D. P. Tew, and D. J. Wales, *npj Quantum Inf* **9**, 1 (2023).

31 D. Mondal, S. Halder, D. Halder, and R. Maitra, “On-the-fly Tailoring towards a Rational Ansatz Design for Digital Quantum Simulations,” (2023), arxiv:2302.03405 [quant-ph].

32 N. Vaquero-Sabater, A. Carreras, R. Orús, N. J. Mayhall, and D. Casanova, *J. Chem. Theory Comput.* **20**, 5133 (2024).

33 Z. Sun, J. Liu, Z. Li, and J. Yang, “Circuit-Efficient Qubit-Excitation-based Variational Quantum Eigensolver,” (2024), arxiv:2406.11699 [quant-ph].

34 P. G. Anastasiou, Y. Chen, N. J. Mayhall, E. Barnes, and S. E. Economou, *Phys. Rev. Res.* **6**, 013254 (2024).

35 C. Feniou, O. Adjoua, B. Claudon, J. Zylberman, E. Giner, and J.-P. Piquemal, *J. Phys. Chem. Lett.* **15**, 3197 (2024).

36 J. Liu, Z. Li, and J. Yang, *Electron. Struct.* **6**, 015007 (2024).

37 H. G. A. Burton, *Phys. Rev. Res.* **6**, 023300 (2024).

38 M. Ramôa, P. G. Anastasiou, L. P. Santos, N. J. Mayhall, E. Barnes, and S. E. Economou, “Reducing the Resources Required by ADAPT-VQE Using Coupled Exchange Operators and Improved Subroutines,” (2024), arxiv:2407.08696.

39 M.-A. Filip, N. Fitzpatrick, D. Muñoz Ramo, and A. J. W. Thom, *Phys. Rev. Res.* **4**, 023243 (2022).

40 A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, *Nature* **549**, 242 (2017).

41 M. Ganzhorn, D. Egger, P. Barkoutsos, P. Ollitrault, G. Salis, N. Moll, M. Roth, A. Fuhrer, P. Mueller, S. Woberner, I. Tavernelli, and S. Filipp, *Phys. Rev. Appl.* **11**, 044092 (2019).

42 J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, *Nat Commun* **9**, 4812 (2018).

43 L. Bittel and M. Kliesch, *Phys. Rev. Lett.* **127**, 120502 (2021).

44 R. D'Cunha, T. D. Crawford, M. Motta, and J. E. Rice, *J. Phys. Chem. A* **127**, 3437 (2023).

45 M. Larocca, S. Thanasilp, S. Wang, K. Sharma, J. Bia monte, P. J. Coles, L. Cincio, J. R. McClean, Z. Holmes, and M. Cerezo, *Nat Rev Phys* **7**, 174 (2025).

46 L. Leone, S. F. E. Oliviero, L. Cincio, and M. Cerezo, *Quantum* **8**, 1395 (2024).

47 A. Pérez-Salinas, H. Wang, and X. Bonet-Monroig, *npj Quantum Inf.* **10**, 27 (2024).

48 B. T. Gard, L. Zhu, G. S. Barron, N. J. Mayhall, S. E. Economou, and E. Barnes, *npj Quantum Inf* **6**, 1 (2020).

49 K. Seki, T. Shirakawa, and S. Yunoki, *Phys. Rev. A* **101**, 052340 (2020).

50 G.-L. R. Anselmetti, D. Wierichs, C. Gogolin, and R. M. Parrish, *New J. Phys.* **23**, 113010 (2021).

51 D. Lacroix, E. A. Ruiz Guzman, and P. Siwach, *Eur. Phys. J. A* **59**, 3 (2023).

52 E. A. Ruiz Guzman and D. Lacroix, *Eur. Phys. J. A* **60**, 112 (2024).

53 I. G. Ryabinkin, T.-C. Yen, S. N. Genin, and A. F. Iz maylov, *J. Chem. Theory Comput.* **14**, 6317 (2018).

54 E. Grant, L. Wossnig, M. Ostaszewski, and M. Benedetti, *Quantum* **3**, 214 (2019).

55 K. M. Nakanishi, K. Fujii, and S. Todo, *Phys. Rev. Res.* **2**, 043158 (2020).

56 S. H. Sack, R. A. Medina, A. A. Michailidis, R. Kueng, and M. Serbyn, *PRX Quantum* **3**, 020365 (2022).

57 B. Choy and D. J. Wales, *J. Chem. Theory Comput.* **19**, 1197 (2023).

58 A. M. Alverlis, A. Khan, T. Iadecola, P. P. Orth, and N. Tubman, “Classical Benchmarks for Variational Quantum Eigensolver Simulations of the Hubbard Model,” (2024), arxiv:2408.00836 [cond-mat, physics:quant-ph].

59 L. Zambrano, A. D. Muñoz-Moller, M. Muñoz, L. Pereira, and A. Delgado, *Quantum Sci. Technol.* **9**, 025016 (2024).
60 V. Novák, I. Zelinka, and V. Snášel, arXiv:2506.01715 (2025).

61 Y. Yao and Y. Hasegawa, Physical Review A 111, 022426 (2025).

62 Y. Liu, G. Huang, Y. Wang, and J. Wu, arXiv:2506.08395 (2025).

63 E. Zapusek, I. Rojkov, and F. Reiter, arXiv:2507.02043 (2025).

64 J. Cunningham and J. Zhuang, *Quantum Information Processing* 24, 1 (2025).

65 N. Yoshioka, H. Hakushima, Y. Matsuzaki, Y. Tokunaga, Y. Suzuki, and S. Endo, Phys. Rev. Lett. 129, 020502 (2022).

66 T. Takeshita, N. C. Rubin, Z. Jiang, E. Lee, R. Babbush, and J. R. McClean, Phys. Rev. X 10, 011004 (2020).

67 W. J. Huggins, J. Lee, U. Baek, B. O'Gorman, and K. B. Whaley, New J. Phys. 22, 073009 (2020).

68 M. Urbanek, D. Camps, R. Van Beeumen, and W. A. de Jong, J. Chem. Theory Comput. 16, 5425 (2020).

69 U. Baek, D. Hait, J. Shee, O. Leimkuhler, W. J. Huggins, T. F. Stetina, M. Head-Gordon, and K. B. Whaley, PRX Quantum 4, 030307 (2023).

70 C. Umeano, F. Jamet, L. P. Lindoy, I. Rungger, and O. Kyriienko, Phys. Rev. Materials 9, 034401 (2025).

71 J. I. Colless, V. V. Ramasesh, D. Dahlen, M. S. Blok, M. E. Kimchi-Schwartz, J. R. McClean, J. Carter, W. A. de Jong, and I. Siddiqi, Phys. Rev. X 8, 011021 (2018).

72 K. Bharti and T. Haug, Phys. Rev. A 104, L050401 (2021).

73 J. C. Getelina, P. Sharma, T. Iadecola, P. P. Orth, and Y.-X. Yao, "Quantum subspace expansion in the presence of hardware noise," (2024), arxiv:2404.09132 [quant-ph].

74 P. J. Ollitrault, A. Kandala, C.-F. Chen, P. K. Barkoutsos, A. Mezzacapo, M. Pistoia, S. Sheldon, S. Woerner, J. M. Gambetta, and I. Tavernelli, Phys. Rev. Res. 2, 043140 (2020).

75 A. Asthana, A. Kumar, V. Abraham, H. Grimsley, Y. Zhang, L. Cincio, S. Tretiak, P. A. Dub, S. E. Economou, E. Barnes, and N. J. Mayhall, Chem. Sci. 14, 2405 (2023).

76 B. N. Parlett, *The Symmetric Eigenvalue Problem*, Classics in Applied Mathematics (Society for Industrial and Applied Mathematics, 1998).

77 S. McArdle, T. Jones, S. Endo, Y. Li, S. C. Benjamin, and X. Yuan, npj Quantum Inf 5, 1 (2019).

78 M. Motta, C. Sun, A. T. Tan, M. J. O'Rourke, E. Ye, A. J. Minnich, F. G. Brandão, and G. K.-L. Chan, Nat. Phys. 16, 205 (2020).

79 K. Yeter-Aydeniz, R. C. Pooser, and G. Siopsis, npj Quantum Inf 6, 1 (2020).

80 W. Kirby, M. Motta, and A. Mezzacapo, *Quantum* 7, 1018 (2023).

81 Z. Zhang, A. Wang, X. Xu, and Y. Li, *Quantum* 8, 1438 (2024).

82 O. Kyriienko, npj Quantum Inf 6, 1 (2020).

83 N. V. Tkachenko, L. Cincio, A. I. Boldyrev, S. Tretiak, P. A. Dub, and Y. Zhang, *Quantum Sci. Technol.* 9, 035012 (2024).

84 R. M. Parrish and P. L. McMahon, "Quantum Filter Diagonalization: Quantum Eigendecomposition without Full Quantum Phase Estimation," (2019), arxiv:1909.08925 [quant-ph].

85 J. Cohn, M. Motta, and R. M. Parrish, PRX Quantum 2, 040352 (2021).

86 C. L. Cortes and S. K. Gray, Phys. Rev. A 105, 022417 (2022).

87 E. N. Epperly, L. Lin, and Y. Nakatsukasa, SIAM J. Matrix Anal. Appl. 43, 1263 (2022).

88 K. Klymko, C. Mejuto-Zaera, S. J. Cotton, F. Wudarski, M. Urbanek, D. Hait, M. Head-Gordon, K. B. Whaley, J. Moussa, N. Wiebe, W. A. de Jong, and N. M. Tubman, PRX Quantum 3, 020323 (2022).

89 Y. Shen, K. Klymko, J. Sud, D. B. Williams-Young, W. A. de Jong, and N. M. Tubman, Quantum 7, 1066 (2023).

90 N. H. Stair, C. L. Cortes, R. M. Parrish, J. Cohn, and M. Motta, Phys. Rev. A 107, 032414 (2023).

91 W. Kirby, Quantum 8, 1457 (2024).

92 N. Yoshioka, M. Amico, W. Kirby, P. Jurcevic, A. Dutt, B. Fuller, S. Garion, H. Haas, I. Hamanura, A. Ivrii, R. Mujumbar, Z. Minev, M. Motta, B. Pokharel, P. Rivero, K. Sharma, C. J. Wood, A. Javadi-Abhari, and A. Mezzacapo Nat Commun 16, 5014 (2025).

93 An extension of that proof to non-Hermitian matrices is possible. Indeed, applying the Householder transformation $R$ to a non-Hermitian matrix $H$ leads to its lower Hessenberg form $\tilde{H} = RHR$. Moreover, since $R$ is unitary, we have $\tilde{H}^\dagger = RH^\dagger R$, which corresponds to the upper Hessenberg form of $H^\dagger$. Following the proof for Hermitian matrix, we can use that $[\tilde{H}^\dagger \tilde{H}]_{00} = [\tilde{H}^\dagger(\theta) \tilde{H}(\theta)]_{00}$ such that $[\tilde{H}^\dagger]_{01}[\tilde{H}]_{10} \geq [\tilde{H}^\dagger(\theta)]_{01}[\tilde{H}(\theta)]_{10}$.

94 S. Sekaran, M. Tsuchiizu, M. Saubanère, and E. Fromager, Phys. Rev. B 104, 035121 (2021), arxiv:2103.04194 [cond-mat].

95 S. Yalouz, S. Sekaran, E. Fromager, and M. Saubanère, *The Journal of Chemical Physics* 157, 214112 (2022).

96 Q. Marécat, B. Lasorne, E. Fromager, and M. Saubanère, Phys. Rev. B 108, 155119 (2023).

97 Q. Marécat and M. Saubanère, *Computation* 11, 203 (2023).

98 J. H. Wilkinson, *The Algebraic Eigenvalue Problem*, Monographs on Numerical Analysis (Clarendon press, Oxford, 1992).

99 G. H. Golub and C. F. Van Loan, *Matrix Computations*, 3rd ed., Johns Hopkins Studies in the Mathematical Sciences (Johns Hopkins University Press, Baltimore, 1996). G. H. Golub and H. A. van der Vorst, *Journal of Computational and Applied Mathematics Numerical Analysis* 2000. Vol. III: Linear Algebra, 123, 35 (2000).

100 Š. Pick and M. Tomášek, *Computer Physics Communications* 27, 101 (1982).

101 F. Orts, G. Ortega, and E. M. Garzón, in *Computational Science - ICCS 2022*, edited by D. Groen, C. de Mulatier, M. Paszynski, V. V. Krzhizhanovskaya, J. J. Dongarra, and P. M. A. Sloot (Springer International Publishing, Cham, 2022) pp. 122–128.

102 B. Claudon, J. Zylberman, C. Feniou, F. Debbasch, A. Peruzzo, and J.-P. Piquemal, Nat Commun 15, 5886 (2024).