# Reduced basis surrogates for quantum spin systems based on tensor networks

Paul Brehmer<EMAIL_ADDRESS>Institute for Theoretical Solid
State Physics, RWTH Aachen University, Otto-Blumenthal-Str. 26, 52074 Aachen,
Germany  Michael F. Herbst<EMAIL_ADDRESS>Mathematics for Materials
Modelling, Institute of Mathematics & Institute of Materials, EPFL, CH-1015
Lausanne  Stefan Wessel<EMAIL_ADDRESS>Institute for
Theoretical Solid State Physics, RWTH Aachen University, Otto-Blumenthal-Str.
26, 52074 Aachen, Germany  Matteo Rizzi<EMAIL_ADDRESS>Forschungszentrum Jülich GmbH, Institute of Quantum Control, Peter Grünberg
Institut (PGI-8), 52425 Jülich, Germany  Institute for Theoretical Physics,
University of Cologne, D-50937 Köln, Germany  Benjamin Stamm<EMAIL_ADDRESS>stuttgart.de Institute of Applied Analysis and Numerical Simulation,
University of Stuttgart, 70569 Stuttgart, Germany

###### Abstract

Within the reduced basis methods approach, an effective low-dimensional
subspace of a quantum many-body Hilbert space is constructed in order to
investigate, e.g., the ground-state phase diagram. The basis of this subspace
is built from solutions of snapshots, i.e., ground states corresponding to
particular and well-chosen parameter values. Here, we show how a greedy
strategy to assemble the reduced basis and thus to select the parameter points
can be implemented based on matrix-product-state (MPS) calculations. Once the
reduced basis has been obtained, observables required for the computation of
phase diagrams can be computed with a computational complexity independent of
the underlying Hilbert space for any parameter value. We illustrate the
efficiency and accuracy of this approach for different one-dimensional quantum
spin-1 models, including anisotropic as well as biquadratic exchange
interactions, leading to rich quantum phase diagrams.

##  I Introduction

A central topic in modern condensed matter theory is the exploration of
ground-state phase diagrams of quantum many-body systems. They uncover the
rich collective behavior of various physical models and often harbor
interesting states of matter that emerge from interaction effects. However,
only a small number of many-body Hamiltonians are solvable analytically, and
in most cases one needs to resort to computational methods to explore strongly
correlated quantum systems. Moreover, in order to relate microscopic models to
experimental findings, it is often required to extend beyond the most basic
model Hamiltonians that describe an emerging quantum many-body effect by
taking into account additional interaction terms or anisotropy effects. An
illustrative example from quantum magnetism is provided by its most basic
model system — the one-dimensional spin chain. In particular, it is by now
well established, dating back to the seminal works by Haldane [1, 2, 3], that
the spin-1 Heisenberg chain exhibits a gapped, quantum-disordered ground state
with symmetry-protected topological order [4, 5]. A quantitative description
of typical Haldane spin chain materials requires to account for both uniaxial
and rhombic-type single-ion anisotropies, yielding a rich overall ground-state
phase diagram [6, 7].

In the numerical study of complex phase diagrams of quantum many-body systems
we then face an exacerbated computational challenge: Not only is solving the
Hamiltonian already demanding due the curse of dimensionality, but the
solutions have to be obtained on many parameter points. Typically, this leads
either to the compromise of resorting to low-resolution scans of the phase
diagram or to the brute-force way of solving the problem in massively parallel
computations. When these approaches become impractical, we advocate here that
so-called reduced basis (RB) methods [8, 9] provide a powerful third option,
as we started exploring recently [10]. Harnessing the linear dependence of
eigenstates across the phase diagram, a low-dimensional surrogate model is
constructed based on a few select ground states at different parameter points.
Using the surrogate, one is then able to evaluate observables independent of
the Hilbert space dimension on any point in the phase diagram, providing
accuracy and efficiency at the same time.

The RB approach for parametrized eigenvalue problems originated from the
structural analysis of mechanical systems in the late nineties [11, 12]. It
has since then migrated to the applied mathematics community with the
pioneering work presented in Ref. [13]. While the RB method became popular for
studying parametrized partial differential equations (PDEs) [8, 9], with many
engineering applications and maturity for industrial applications, the RB
method for parametrized eigenvalue problems had not received much attention
thus far. Only years later, the RB methodology developed in Ref. [13] was
further extended and generalized to account for better error control and more
general settings [14, 15, 16, 17].

In the context of quantum many-body physics, RB methods for eigenvalue
problems found application only recently and are most prominent in the nuclear
physics community [18, 19, 20, 21, 22, 23], where their use emerged with the
eigenvector continuation (EC) method [24, 25, 26]. EC was retrospectively
identified as belonging to the family of RB methods, and more broadly, to the
field of model order reduction [27, 23]. In further physical applications, RB
surrogate models, referred to as emulators in EC language, proved viable in
determining phase diagrams of quantum spin systems [10]. Furthermore, EC
emulators were used as a subspace diagonalization method on quantum computers
[28], for emulating superconducting phenomena [29, 30], as well as in quantum
chemistry applications [31].

In this article, we expand upon our previous work [10] that utilizes exact
diagonalization (ED) techniques in combination with a greedy RB approach. As
suggested therein, we apply the density matrix renormalization group (DMRG)
[32, 33] to treat spin chain Hamiltonians and perform all necessary vector
operations using MPSs [34, 35, 36, 37], which belong to the family of tensor
network states [38, 39] that our approach could be naturally extended to (at
least for finite system sizes). While this approach has been very recently
adopted in [29], we focus on methodological aspects and develop the combined
RB-MPS method in application to quantum spin systems. In using compressed
MPSs, we gain access to larger many-body systems, allowing us to probe whether
RB methods remain a viable approach in the light of larger Hilbert spaces and
more complex parameter domains. We made a related open-source software package
available [40].

The remainder of this paper is organized as follows: In Sec. II we review the
greedy RB algorithm and introduce notation, before explaining the combined RB-
MPS approach. This allows us to apply the method to different spin-1 chain
Hamiltonians in Sec. III, namely the Haldane chain with uniaxial and rhombic
single-ion anisotropies as well as the bilinear-biquadratic model with a
uniaxial anisotropy, where we scan the ground-state phase diagrams for various
correlations. Section IV discusses the accuracy of the surrogate models in
different settings as well as the effects of MPS approximations on the RB
framework, after which we conclude in Sec. V.

##  II Method

###  II.1 The RB framework

#### Problem setting.

We begin by defining the physical problems that we aim to treat in the RB
approach. To that end, we consider a generic stationary quantum many-body
problem

$H(\bm{\mu})\ket{\bm{\Psi}(\bm{\mu})}=E(\bm{\mu})\ket{\bm{\Psi}(\bm{\mu})},$
(1)

that consists of finding the ground-state energy $E(\bm{\mu})$ and its
corresponding $m$ ground states
$\ket{\bm{\Psi}(\bm{\mu})}=(\ket{\Psi^{[1]}(\bm{\mu})},\dots,\ket{\Psi^{[m]}(\bm{\mu})})$.
The many-body Hilbert space $\mathcal{H}=\mathbb{C}^{\mathcal{N}}$ under
consideration is assumed finite but high-dimensional since $\mathcal{N}\gg 1$
grows exponentially with the number of physical constituents. We assume the
eigenvalue problem to be parametrized by a vector of physical model parameters
$\bm{\mu}\in\mathbb{P}$ that resides in the parameter space $\mathbb{P}$,
i.e., for each $\bm{\mu}$ we obtain a new Hamiltonian to solve. In order to
apply the RB framework, we consider a particular class of Hamiltonians that
can be expressed as affine decompositions

$H(\bm{\mu})=\sum_{q=1}^{Q}\theta_{q}(\bm{\mu})\,H_{q},$ (2)

with a number of terms $Q$ independent of the Hilbert space dimension
$\mathcal{N}$. Such linear combinations distinguish between parameter-
dependent coefficients $\theta_{q}:\mathbb{P}\to\mathbb{R}$ and parameter-
independent Hermitian matrices $H_{q}:\mathcal{H}\to\mathcal{H}$. Not only the
Hamiltonian must be affinely decomposable, but in fact all observables that we
may want to measure are assumed to be of the form

$O(\bm{\mu};p)=\sum_{r=1}^{R}\alpha_{r}(\bm{\mu};p)\,O_{r},$ (3)

in terms of operators $O_{r}:\mathcal{H}\to\mathcal{H}$, and where the
coefficients $\alpha_{r}(\bm{\mu};p)\in\mathbb{C}$ now include additional
parameters $p$ that are separate from $\mathbb{P}$. This generalization is
required, e.g., when studying observables in Fourier space. Note that the
individual operators $O_{r}$ need not necessarily be Hermitian. To then
extract physical information, one measures affine decompositions by taking the
expectation value

$\braket{O(\bm{\mu};p)}=\sum_{r=1}^{R}\alpha_{r}(\bm{\mu};p)\,\frac{1}{m}\sum_{i=1}^{m}\braket{\Psi^{[i]}(\bm{\mu})}{O_{r}}{\Psi^{[i]}(\bm{\mu})},$
(4)

where we average over the degenerate ground-state subspace, if needed.

In the particular setting we consider here, the goal is to scan domains of
$\mathbb{P}$ with fine resolution, i.e., to solve the many-body problem in Eq.
(1) on a large set of parameter points and then, based on the solutions,
measure observables on the same parameter domain. This is addressed by the RB
approach in two steps: In the first step, which is referred to as the offline
stage in RB parlance, a surrogate model is assembled based on a small number
of ground states across the parameter space, whereas in the second step, the
so-called online stage, observable measurements are efficiently obtained. This
approach is warranted by the insight that ground states at different parameter
points tend to show significant linear dependence. While being demonstrated in
various numerical applications (see e.g. [10, 23]), this can also be reasoned
in terms of analytical continuation in the context of eigenvector continuation
[24, 25]. One particular strategy to build RB surrogate models features a
greedy mechanism [10, 9, 8] to sample the parameter domain, which we want to
review next.

#### Offline stage.

The goal in the offline stage is to construct a low-dimensional reduced basis
space

$\mathbb{V}_{n}\coloneqq\operatorname{span}\big{\\{}\\!\ket{\bm{\Psi}(\bm{\mu}_{1})},\dots,\ket{\bm{\Psi}(\bm{\mu}_{n})}\\!\big{\\}}$
(5)

that is spanned by degenerate subspaces extracted at $n$ different parameter
points $\\{\bm{\mu}_{1},\dots,\bm{\mu}_{n}\\}$. Here we refer to
$\ket{\bm{\Psi}(\bm{\mu}_{j})}=(\ket{\Psi^{[1]}(\bm{\mu}_{j})},\dots,\ket{\Psi^{[m_{j}]}(\bm{\mu}_{j})})$
as a “snapshot”, such that one snapshot contains $m_{j}\geq 1$ states. Since
the snapshots might share linear dependent modes, the reduced basis dimension
$N\coloneqq\dim\mathbb{V}_{n}\leq M\coloneqq\sum_{j=1}^{n}m_{j}$ can in
general be smaller than the total number of obtained ground states $M$. For
the single ground-state solutions that span the RB space we use the shorthand
$\ket{\Psi_{j}}$ where $j=1,\dots,N$. Moreover, the solver used to obtain the
lowest eigenvalue and eigenvectors of $H(\bm{\mu})$ is called the “truth
solver”, which is to be understood in the sense that it can obtain the true
solution of the Hamiltonian at any point in the considered parameter domain,
up to high numerical accuracy. Beyond that, we do not need to further specify
the actual method yet, i.e., the RB approach is solver-agnostic. The idea is
then to represent ground states at any other parameter point
$\bm{\mu}\in\mathbb{P}$ as a linear combination of snapshots

$\ket{\Phi_{\mathrm{rb}}(\bm{\mu})}=\sum_{j=1}^{N}a_{j}(\bm{\mu})\ket{\Psi_{j}},$
(6)

with parameter-dependent coefficients $a_{j}(\bm{\mu})\in\mathbb{C}$. This
prompts the following two questions: First, how do we generate the RB space
$\mathbb{V}_{n}$ and in particular how are the parameter points chosen, and
second, how do we determine the coefficients $a_{j}(\bm{\mu})$?

To answer the first question, we turn to the greedy algorithm. We restrict the
snapshots to a subset of all parameter points, which we call the training grid
$\Xi_{\mathrm{train}}\subset\mathbb{P}$, that preselects the domain on which
the RB will be “trained”. Starting from a suitable first parameter point
$\bm{\mu}_{1}$ and its corresponding snapshot $\ket{\bm{\Psi}(\bm{\mu}_{1})}$,
the RB is constructed in an inductive manner. Namely at iteration $n$, using
the current RB space $\mathbb{V}_{n}$, we select the next parameter point
$\bm{\mu}_{n+1}$ where a truth solve will be performed. Based on the ground
states that are contained in the RB space

$\ket{\bm{\Phi}_{\mathrm{rb}}(\bm{\mu})}=\operatorname*{arg\,min}_{\ket{\Phi}\in\mathbb{V}_{n}}\frac{\braket{\Phi}{H(\bm{\mu})}{\Phi}}{\braket{\Phi}{\Phi}},\quad\forall\bm{\mu}\in\Xi_{\mathrm{train}},$
(7)

solved at all training points, we select the parameter point that meets the
greedy condition

$\bm{\mu}_{n+1}=\operatorname*{arg\,max}_{\bm{\mu}\in\Xi_{\mathrm{train}}}\,\mathrm{Res}_{n}(\bm{\mu}),$
(8)

i.e., which corresponds to the snapshot that maximizes the residual

$\mathrm{Res}_{n}(\bm{\mu})\coloneqq\sqrt{\sum_{i=1}^{m}\big{\lVert}H(\bm{\mu})\ket{\Phi_{\mathrm{rb}}^{[i]}(\bm{\mu})}-E_{\mathrm{rb}}(\bm{\mu})\ket{\Phi_{\mathrm{rb}}^{[i]}(\bm{\mu})}\big{\rVert}^{2}}.$
(9)

In this way, we sample first the parts of $\Xi_{\mathrm{train}}$ where the RB
performs worst, meaning that the surrogate ground states and energies least
accurately fulfill the eigenvalue problem of $H(\bm{\mu})$. Notice the
dependency of $\ket{\bm{\Phi}_{\mathrm{rb}}}$ and $E_{\mathrm{rb}}$ on the
number of greedy iterations since they were obtained from the $n$-snapshot RB
space $\mathbb{V}_{n}$. Following this, we perform a truth solve at
$\bm{\mu}_{n+1}$ and obtain a snapshot $\ket{\bm{\Psi}(\bm{\mu}_{n+1})}$ that
is appended to the basis, forming $\mathbb{V}_{n+1}$. This iteration is
repeated until an adequate exit condition, e.g., a target residual accuracy,
is reached.

In the current formulation, it seems difficult and expensive to solve Eq. (7)
in order to determine the greedy condition, i.e., to compute the residual on
all training points. To clarify this and to answer the question of determining
the linear coefficients of $\ket{\bm{\Phi}_{\mathrm{rb}}(\bm{\mu})}$, we now
map the problem onto the low-dimensional RB space via the reduced basis
$B\\!:\mathbb{V}_{n}\to\mathcal{H}$ by way of the Rayleigh-Ritz method. To
that end, we define the RB coefficients
$\bm{\varphi}_{\mathrm{rb}}(\bm{\mu})\in\mathbb{C}^{N\times m}$ by relating
them to

$\ket{\bm{\Phi}_{\mathrm{rb}}(\bm{\mu})}=B\bm{\varphi}_{\mathrm{rb}}(\bm{\mu})\in\mathbb{C}^{\mathcal{N}\times
m},$ (10)

corresponding to the Ritz vector, as well as the normalization matrix
$b=B^{\dagger}B$ and the reduced Hamiltonian

$h(\bm{\mu})=\sum_{q=1}^{Q}\theta_{q}(\bm{\mu})\,h_{q},\quad
h_{q}=B^{\dagger}H_{q}B.$ (11)

Then we can reformulate the variational problem of Eq. (7) as a generalized
eigenvalue problem

$h(\bm{\mu})\,\bm{\varphi}_{\mathrm{rb}}(\bm{\mu})=E_{\mathrm{rb}}(\bm{\mu})\,b\,\bm{\varphi}_{\mathrm{rb}}(\bm{\mu}),$
(12)

solving only for the lowest eigenvalue and corresponding eigenvectors. Indeed,
by expressing the residual in terms of reduced quantities and using the
eigenvalue Eq. (12), one finds the more efficient expression

$\displaystyle\mathrm{Res}_{n}^{2}(\bm{\mu})=$
$\displaystyle\sum_{q,q^{\prime}=1}^{Q}\theta_{q}(\bm{\mu})\,\theta_{q^{\prime}}(\bm{\mu})\,\sum_{i=1}^{m}\varphi_{\mathrm{rb}}^{[i]}(\bm{\mu})^{\dagger}\,h_{qq^{\prime}}\,\varphi_{\mathrm{rb}}^{[i]}(\bm{\mu})$
$\displaystyle-E^{2}_{\mathrm{rb}}(\bm{\mu})\sum_{i=1}^{m}\varphi_{\mathrm{rb}}^{[i]}(\bm{\mu})^{\dagger}\,b\,\varphi_{\mathrm{rb}}^{[i]}(\bm{\mu}),$
(13)

where $h_{qq^{\prime}}=B^{\dagger}H_{q}H_{q^{\prime}}B$ are the reduced
matrices of the $h^{2}(\bm{\mu})$ affine decomposition. We have thus obtained
the greedy assembly algorithm, where computing the greedy condition boils down
to the $N$-dimensional generalized eigenvalue problem of Eq. (12) and all
Hilbert space dimension dependent operations were isolated to the computation
of $b$, $h_{q}$, $h_{qq^{\prime}}$ and the truth solve, which are performed
only once per iteration.

Finally, to determine $B$ itself, we might make the simple ansatz of using the
snapshots as column vectors
$[\ket{\bm{\Psi}(\bm{\mu}_{1})}\cdots\ket{\bm{\Psi}(\bm{\mu}_{n})}]$. This,
however, leads to a poorly-conditioned normalization $b$, due to the
increasing linear dependence between the columns as we add snapshots, and
quickly renders the solution of Eq. (12) numerically unfeasible. For that
reason, one typically orthogonalizes the RB such that $b\simeq I$. The
orthogonalization procedure can be implemented on the level of the
coefficients $\bm{\varphi}_{\mathrm{rb}}$ or directly on the snapshots — this
depends on the truth solver and in particular the associated vector format,
which we leave unspecified for now.

#### Online stage.

Once the offline stage is finished, we are left with the low-dimensional basis
$B$ and the reduced quantities $b$, $h(\bm{\mu})$ and $h^{2}(\bm{\mu})$ as
byproducts of the greedy assembly algorithm. In order to compute expectation
values of affine decompositions $O(\bm{\mu};p)$, we compute the reduced
matrices $o_{r}=B^{\dagger}O_{r}B$, so that we now operate entirely in RB
space. Namely, by taking the expectation value with the emulated ground state
$\ket{\bm{\Phi}_{\mathrm{rb}}(\bm{\mu})}$, we obtain

$\braket{O(\bm{\mu};p)}_{\mathrm{rb}}=\sum_{r=1}^{R}\alpha_{r}(\bm{\mu};p)\,\frac{1}{m}\sum_{i=1}^{m}\varphi_{\mathrm{rb}}^{[i]}(\bm{\mu})^{\dagger}\,o_{r}\,\varphi^{[i]}_{\mathrm{rb}}(\bm{\mu}).$
(14)

It becomes clear that we can, as advertised, evaluate the above on any
parameter point $\bm{\mu}\in\mathbb{P}$, independent of $\mathcal{N}$ — thanks
to the affine decomposition, $o_{r}$ does not depend on $\bm{\mu}$ and we only
need to evaluate the coefficient functions and compute the RB coefficients
$\bm{\varphi}_{\mathrm{rb}}(\bm{\mu})$. This again amounts to solving the
$N$-dimensional generalized eigenvalue problem of Eq. (12). For further
details on the greedy RB approach, we refer to Refs. [10, 9, 8].

###  II.2 RB assembly using MPSs

We now aim to combine the RB framework with tensor network techniques from
quantum many-body physics. More specifically, we represent the snapshot many-
body ground states as MPSs

$\ket{\psi}=\\!\sum_{\alpha_{1}\cdots\alpha_{L}}\sum_{a_{1}\cdots
a_{L+1}}\\!\\!\\!M^{\alpha_{1}}_{a_{1}a_{2}}M^{\alpha_{2}}_{a_{2}a_{3}}\cdots
M^{\alpha_{L}}_{a_{L}a_{L+1}}\ket{\alpha_{1}\cdots\alpha_{L}},$ (15)

for systems of $L$ physical degrees of freedom $\alpha_{i}$. In doing so, the
treatment of large many-body Hamiltonians becomes feasible since efficient
algorithms for the computation of state overlaps, matrix elements and ground-
state MPSs exist, that scale polynomially in $L$. These rely on low-rank
approximations of the MPS tensors $M^{\alpha_{i}}_{a_{i}a_{i+1}}$ which can be
implemented by repeated singular value decompositions to reduce the matrix
rank $d$, known as the bond dimension in tensor network theory [39], by
dropping singular values according to a cutoff

$\mathtt{cut}_{\sigma}>\sqrt{\frac{\sum_{k\in\mathrm{trunc}}\sigma_{k}^{2}}{\sum_{k=1}^{d}\sigma_{k}^{2}}}.$
(16)

The sum $\sum_{k\in\mathrm{trunc}}$ is to be understood in the sense that we
remove the smallest singular values $\sigma_{k}$ once the singular value error
surpasses $\mathtt{cut}_{\sigma}$, which corresponds to a truncation of the
MPS tensors in the Frobenius norm. In order to perform truth solves in the MPS
format, we use DMRG, which is the most commonly used method for variational
ground MPS searches. For reviews on general MPS theory and the DMRG algorithms
we refer the reader to Refs. [38, 37, 39].

Carrying out the greedy offline stage in MPS formulation spawns two new
aspects: First, we need to use efficient MPS contractions for overlaps and
matrix elements when computing the reduced quantities $b$, $h$ and $h^{2}$
which generates additional inaccuracies on top of the truth solve, and second,
we need an orthogonalization scheme for $B$ that is compatible with MPS
snapshots. Here, we elaborate on these aspects by going through the greedy
algorithm once again, and postpone the discussion of MPS inaccuracies in the
RB framework to Sec. IV.2.

For each iteration of the greedy assembly, we obtain $m_{j}$ normalized
ground-state MPSs and append them to the matrix
$\Upsilon=[\ket{\Psi_{1}}\cdots\ket{\Psi_{M}}]\in\mathbb{C}^{\mathcal{N}\times
M}$ as column vectors. In practice, this matrix cannot be constructed
explicitly and only operations between its columns are allowed. To
orthogonalize this matrix, we make the ansatz $B=\Upsilon V$, where
$V\in\mathbb{C}^{M\times N}$ mixes the truth MPSs into orthogonal linear
combinations. However, these linear combinations are not computed explicitly
since this would entail inefficient MPS addition. Instead, we first compute
the overlap matrix

$S\coloneqq\Upsilon^{\dagger}\Upsilon,\quad
S_{ij}=\braket{\Psi_{i}}{\Psi_{j}},$ (17)

and matrix elements
$[\Upsilon^{\dagger}A\Upsilon]_{ij}=\braket{\Psi_{i}}{A}{\Psi_{j}}$ where
$i,j\in\\{1,\dots,M\\}$ and $A=H_{q},H_{q}H_{q^{\prime}},O_{r}$ might be any
operator of interest. Note that we assume the operator $A$ to be represented
in tensor format, e.g., as a matrix product operator (MPO) or a multi-site
operator. Then all reduced quantities are computed in the orthogonal basis by
transforming in RB space as

$b=V^{\dagger}SV,\quad a=V^{\dagger}\Upsilon^{\dagger}A\Upsilon V,$ (18)

where correspondingly $a=h_{q},h_{qq^{\prime}},o_{r}$. Note that by adding a
new MPS snapshot $\ket{\bm{\Psi}(\bm{\mu}_{n+1})}$ to $\Upsilon$, we need to
compute $M+m_{n+1}$ new overlaps and matrix elements per observable,
exploiting the hermiticity of $S$ and $\Upsilon^{\dagger}A\Upsilon$ (granted
that $A$ is Hermitian).

While there are multiple options to determine the orthogonalizing matrix $V$,
we here opt for the numerically efficient approach of Löwdin symmetric
orthogonalization [41]. To that end, we decompose the overlap matrix
$S=U\Lambda U^{-1}$ into its eigenvalues
$\Lambda=\mathrm{diag}(\lambda_{1},\dots,\lambda_{M})$, sorted in descending
order. Since $S$ is Hermitian, $U$ can be chosen to be a unitary matrix, so
that by demanding $b=V^{\dagger}SV\overset{!}{=}I$, we can immediately
identify

$V=U\Lambda^{-1/2}.$ (19)

Furthermore, the eigenvalue decomposition of $S$ provides a way to compress
the RB. Due to accumulating linear dependence, the eigenvalues $\lambda_{j}$
decrease as snapshots are appended until we cannot further resolve new
directions in RB space, given the truth solver’s accuracy. Hence we may
truncate $B$ in the Frobenius norm of $S$ according to an eigenvalue cutoff

$\mathtt{cut}_{\lambda}>\sqrt{\frac{\sum_{j\in\mathrm{trunc}}\lambda_{j}^{2}}{\sum_{j=1}^{M}\lambda_{j}^{2}}},$
(20)

i.e., when the smallest normalized sum of squared eigenvalues exceeds
$\mathtt{cut}_{\lambda}$ the corresponding snapshots are removed. This
approach of orthogonalizing a set of snapshots is akin to what is known as
proper orthogonal decomposition [9, 8] in RB theory, which provides an
(expensive) alternative to the greedy sampling algorithm.

On a slightly more technical note, we mention the possibility of using the RB
to produce initial guesses for the DMRG solver. By explicitly computing the
linear combination
$\ket{\Phi^{[i]}_{\mathrm{rb}}(\bm{\mu})}=B\varphi^{[i]}_{\mathrm{rb}}(\bm{\mu})=\sum_{j=1}^{N}[V\varphi^{[i]}_{\mathrm{rb}}(\bm{\mu})]_{j}\ket{\Psi_{j}}$
for $i=1,\dots,m$ at the selected parameter point and using it as the initial
MPS, one can speed up DMRG convergence. In order for this to work sufficiently
fast, it is necessary to heavily truncate the MPS while computing the linear
combination, since MPS addition leads to an additive increase of bond
dimensions. Nonetheless, this approach globally reduces the number of DMRG
sweeps and makes the RB assembly more stable as well as deterministic.

Before proceeding with the numerical results, we remark that the procedures we
covered here in the context of MPSs generalize to a larger class of vector
representations and truth solvers. In summary, one can identify the following
requirements for compatibility with the RB framework:

  1. 1.

Computation of state overlaps
$\langle\Psi(\bm{\mu}_{i})|\Psi(\bm{\mu}_{j})\rangle$.

  2. 2.

Computation of matrix elements
$\langle\Psi(\bm{\mu}_{i})|A|\Psi(\bm{\mu}_{j})\rangle$ for all relevant
operators $A=H_{q},H_{q}H_{q^{\prime}},O_{r}$.

  3. 3.

High truth solver accuracy: Large approximation errors on overlaps and matrix
elements prevent us from generating a meaningful surrogate. Since this turns
out to be a subtle point, we will further discuss this for MPSs in Sec. IV.2.

In consequence, we do not necessarily need access to the truth ground-state
vectors, only to contractions between them. This opens the door to further
tensor network architectures such as projected entangled pair states [42] and
tree tensor networks [43, 44, 45], or possibly, artificial neural network
based representations [46] as well as various truth solving methods, e.g.,
quantum Monte Carlo approaches [47].

##  III Results

We now turn to discuss various numerical results obtained using the RB-MPS
method. Alongside this work, we developed a code package implementing RB
methods for parametrized eigenvalue problems, in particular many-body
Hamiltonians, with the possibility for using DMRG and ED-based solvers as well
as custom truth solving methods. The code package is written in the Julia
programming language [48] and all MPS and DMRG procedures are performed using
the ITensor library [49, 50]. We made the code, including user instructions
and documentation, publicly available [40].

###  III.1 Haldane spin-1 chain with single-ion anisotropies

Figure 1: Ground-state phase diagram of the spin-1 chain in Eq. (21), where
the phase boundaries are taken from Ref. [6]. In the Néel phases the model
exhibits antiferromagnetic order in $x$, $y$ and $z$-direction, respectively.
The Haldane phase is gapped and an example of symmetry-protected topological
order and exhibits non-local string order. The uniaxial anisotropic coupling
leads to the large-$D$ phase that confines the spins to the $xy$-plane, since
the $D$-term penalizes configurations which are polarized in $z$-direction.
Similarly, due to the rhombic coupling the large-$E_{x}$ and $E_{y}$ phases,
which differ in the sign of $E$, favor configurations in the $yz$ and
$xz$-plane.

As a first application of the RB-MPS method, we consider the one-dimensional
Haldane spin-1 chain

$H_{\mathrm{HD}}\\!=\\!J\sum_{i=1}^{L-1}\bm{S}_{i}\cdot\bm{S}_{i+1}+D\sum_{i=1}^{L}(S_{i}^{z})^{2}+E\sum_{i=1}^{L}\big{[}(S_{i}^{x})^{2}-(S_{i}^{y})^{2}\big{]},$
(21)

with a uniaxial $D$ and a rhombic-type $E$ single-ion anisotropy, here using
open boundary conditions. With regard to the RB formalism we express
$H_{\mathrm{HD}}$ as a dimensionless affine decomposition using the parameter
vector $\bm{\mu}=(D/J,\,E/J)$ and the coefficient function
$\bm{\theta}(\bm{\mu})=(1,\mu_{1},\mu_{2})$, such that the matrices $H_{q}$
correspond to the summands of $H_{\mathrm{HD}}$. In order to benchmark the RB-
MPS method, we want to resolve the ground-state phase diagram of
$H_{\mathrm{HD}}$, which features: i) the symmetry-protected topological
Haldane phase [1, 2, 3, 4, 5] that is robust against small anisotropies, ii)
Néel-ordered phases, as well as so-called iii) large-$D$ and iv) large-$E$
phases. We summarize and sketch the full phase diagram in Fig. 1: cf. Refs.
[6, 7] for a detailed discussion and an overview of previous investigations.
DMRG studies of $H_{\mathrm{HD}}$ can profit from several parity quantum
numbers [51]. In particular, $H_{\mathrm{HD}}$ conserves the magnetization
parity $p_{m}=\sum_{i=1}^{L}S_{i}^{z}\\!\\!\mod 2\in\\{0,1\\}$ together with
the spatial parity $p_{s}=\pm 1$ and time reversal $t=\pm 1$. The ground-state
sector is described by $p_{m}=0$, whereas the corresponding $p_{s}$ and $t$
can be determined based on the valence-bond-solid picture [52] in dependence
of the boundary conditions [6]. In our RB-MPS approach we thus use DMRG with
Abelian quantum number conservation and fix $p_{m}=0$ to operate in the
ground-state sector. Note that $p_{s}$ and $t$ are not explicitly fixed here.

For the RB-MPS calculations we consider spin chains of length
$L\in\\{20,40,60,80\\}$ using MPSs with a singular value error cutoff 111Note
that in the ITensor convention, the singular value cutoffs are defined via the
squared Frobenius norm, such that the ITensor cutoffs correspond to
$\mathtt{cut}_{\sigma}^{2}$. of up to $\mathtt{cut}_{\sigma}=10^{-8}$ —
corresponding to bond dimensions of mostly $d\sim 100$ up to 1000, depending
on the phases — and DMRG sweeps which are converged up to an energy tolerance
of $10^{-9}$. Furthermore, we focus on the upper half-plane
$\mathbb{P}_{\mathrm{HD}}=[-2,2]\times[0,2]$ since the remaining phase diagram
is obtained by flipping the sign of $E$ which corresponds to swapping $x$ and
$y$ terms. On this domain, we use a regular training grid of $80\times 60$
parameter points and converge the surrogate up to an overlap eigenvalue
accuracy of at least $\mathtt{cut}_{\lambda}<10^{-8}$. Note that we here and
in all following examples do not target degenerate ground-states with DMRG
solves, such that the number of snapshots $n=N$ corresponds to the RB
dimension.

Figure 2: Maximal residual (upper panels) and the normalized minimal
eigenvalue
$\tilde{\lambda}^{\mathrm{min}}=\lVert\Lambda\rVert^{-1}_{F}\min(\Lambda)$ of
the overlap matrix $S$ (lower panels) in dependence of the RB dimension $N$.
Both quantities follow an exponential decay, although with different rates.
The right panels show the curves on rescaled $x$-axes which illustrates the
scaling behavior of $N$ with system size. More specifically, we find $N\sim
L^{\eta}$ with approximately $\eta=0.8$ for the residual and $\eta=0.65$ for
the minimal eigenvalues. Figure 3: Residual error estimate
$\mathrm{Res}(\bm{\mu})$ on $\mathbb{P}_{\mathrm{HD}}$. The dots indicate the
$N=173$ parameter points $\\{\bm{\mu}_{1},\dots,\bm{\mu}_{N}\\}$ corresponding
to the snapshots that span the RB space. The dotted white lines show the phase
boundaries of [6].

At these settings, the greedy algorithm takes $N=173$ snapshots to reach
convergence at the largest system size of $L=80$, with maximal residuals of
order $\max_{\bm{\mu}\in\Xi_{\mathrm{train}}}\mathrm{Res}(\bm{\mu})\sim
10^{-2}$. Only a fraction of the $80\times 60$ possible training points is
thus needed to converge the surrogate model. To illustrate the convergence, we
show the residual as well as the decay of the minimal eigenvalue of the
overlap matrix $S$ as a function of the RB dimension $N$ in Fig. 2. The
exponential decay of the overlap eigenvalues numerically demonstrates the
increasing linear dependence among ground states on $\mathbb{P}_{\mathrm{HD}}$
as snapshots are appended to the RB. Viewed from a physical angle, the fact
that the decrease of overlap eigenvalues coincides with decreasing residuals
over the parameter domain can be interpreted as larger eigenvalues being
associated with global ground-state behavior, while smaller eigenvalues
resolve more localized features in the phase diagram. Unsurprisingly, the
decay rates decrease with growing system sizes — in the thermodynamic limit
$L\to\infty$ we would in fact run into an orthogonality catastrophe, requiring
a diverging number of snapshots to assemble a RB. Note that this also hinders
the immediate use of infinite MPSs [54, 55, 56] in the RB framework, since all
overlaps and matrix elements between ground states at different parameter
points would vanish, effectively prohibiting any linear combinations between
snapshots. However, for $H_{\mathrm{HD}}$ on $\mathbb{P}_{\mathrm{HD}}$, the
RB dimension required to converge a surrogate up to a fixed residual increases
merely sublinearly in $L$, as opposed to the exponential Hilbert space growth,
which indicates that there exists a sweet spot for system sizes, where the
thermodynamic regime is approached while still being amenable to the RB-MPS
method.

It is also interesting to examine the residual on the entire parameter domain,
i.e., on a high-resolution online grid covering $\mathbb{P}_{\mathrm{HD}}$,
together with the snapshot parameter points in Fig. 3. In particular, many
snapshots are selected along phase boundaries and especially around the
transition point between the large-$E_{x}$ and Haldane phase, whereas deep in
the phases less snapshots are needed to resolve the ground-state behavior. We
also find that a noticeable number of snapshots are located along the boundary
of the parameter range $\mathbb{P}_{\mathrm{HD}}$, which is standard in greedy
RB algorithms. The clustering of snapshot parameter points hints at the
selection mechanism of the greedy algorithm that indicates domains in the
phase diagram where the ground states vary more rapidly. More specifically, a
higher density of sample points indicates a slower decay of the local
Kolmogorov $N$-width [57] and thus a higher linear independence of the
solutions under local parameter variation.

Figure 4: Different two-site correlation functions measured on
$\mathbb{P}_{\mathrm{HD}}$ for an $L=80$ chain at sites $r=20$ and
$r^{\prime}=60$. The dotted white lines show the phase boundaries obtained in
[6]. In the upper left and right panels, the Néel phases are clearly indicated
by the spin correlators in $z$ and $y$-direction, whereas the quadrupolar
correlator in the bottom left panel is smoothed along the phase boundaries.
The bottom right panel shows the modified string order parameter from Eq. (22)
where spin-spin correlations in $z$-direction were subtracted.

In order to reproduce the phase diagram of Fig. 1, we next measure various
correlation functions that distinguish the different types of ground-state
order [6]. The Néel-ordered phases are identified by measuring the spin-spin
correlation functions $\langle
S^{\alpha}_{r}S^{\alpha}_{r^{\prime}}\rangle_{\mathrm{rb}}$ in the $\alpha=z$
and $y$-direction, while the large-$D$ and large-$E_{x}$ phases are
characterized by quadrupolar correlations $\langle
Q^{\gamma}_{r}Q^{\gamma}_{r^{\prime}}\rangle_{\mathrm{rb}}$ where we focus on
the $\gamma=z^{2}$ order with
$Q^{z^{2}}_{r}=\big{[}3(S^{z}_{r})^{2}-2\big{]}/\sqrt{3}$. To detect the
Haldane phase, we consider the non-local string order operator [58]
$O_{rr^{\prime}}=-S^{z}_{r}\prod_{j=r+1}^{r^{\prime}-1}e^{i\pi
S^{z}_{j}}S^{z}_{r^{\prime}}$ from which we subtract the spin correlator in
$z$-direction to remove the trivial background signal in the $z$-Néel phase
and isolate the string order in the Haldane phase,

$\langle\tilde{O}_{rr^{\prime}}\rangle_{\mathrm{rb}}=-\big{\langle}S^{z}_{r}e^{i\pi\sum_{j=r+1}^{r^{\prime}-1}S^{z}_{j}}S^{z}_{r^{\prime}}\big{\rangle}_{\mathrm{rb}}-\langle
S^{z}_{r}S^{z}_{r^{\prime}}\rangle_{\mathrm{rb}}.$ (22)

We show the corresponding measurements obtained from the $L=80$ RB in Fig. 4.
The RB correctly reproduces the different ground-state phases, up to finite-
size effects, matching the results from Ref. [6]. It can be clearly seen how
the string order persists in the Haldane phase also for small values of the
anisotropies $D$ and $E$.

Overall, we find that the RB-MPS approach allows us to efficiently uncover the
parameter regions of the various ground-state phases at high resolution, based
on a comparably low number $N\sim 100$ of DMRG truth solves. We leave a
quantitative treatment of the accuracy of RB measurements to Sec. IV and move
to a more complex example application.

###  III.2 Bilinear-biquadratic spin-1 chain with uniaxial single-ion
anisotropy

Figure 5: Phase diagram of the bilinear-biquadratic spin-1 chain with uniaxial
coupling $D$. The phase boundaries are taken from Ref. [59], here shown
without error bars. The $\theta=0$ cut is equivalent to the $E=0$ cut of
$H_{\mathrm{HD}}$ and therefore already contained in Fig. 1. Besides the
Haldane, large-$D$ and $z$-Néel phase, further regimes are the fully polarized
ferromagnetic (FM) and in-plane XY-FM ordered phases, a dimerized phase with a
finite bond order parameter (24), and two critical phases, denoted critical A
and critical B. In the former, the central charge $c=1$, while for the latter,
two domains have been identified, where $c=1$ and $c=2$, respectively.

In the second example we stay in the realm of spin-1 chains. We now add a
biquadratic exchange interaction term, while removing the rhombic anisotropy,
resulting in the bilinear-biquadratic model with a uniaxial single-ion
anisotropy

$\displaystyle H_{\mathrm{BLBQ}}=\ $ $\displaystyle
J\sum_{i=1}^{L-1}\big{[}\cos(\theta)\bm{S}_{i}\cdot\bm{S}_{i+1}+\sin(\theta)(\bm{S}_{i}\cdot\bm{S}_{i+1})^{2}\big{]}$
$\displaystyle+D\sum_{i=1}^{L}(S^{z}_{i})^{2},$ (23)

on a chain with open boundaries. Written as a dimensionless affine
decomposition, we identify the parameter vector $\bm{\mu}=(\theta,\,D/J)$ and
coefficient function
$\bm{\theta}(\bm{\mu})=(\cos(\mu_{1}),\sin(\mu_{1}),\mu_{2})$. Various aspects
of the rather rich physics of this model, featuring multiple gapped as well as
critical phases, were previously studied using both analytical and
computational approaches, including extensive DMRG calculations. We refer in
particular to the overall ground-state phase diagram reported in Ref. [59],
which also provides an overview of previous studies. The ground-state phase
diagram of $H_{\mathrm{BLBQ}}$ is outlined and visualized in Fig. 5. Due to
the increased complexity in parameter space, the anisotropic bilinear-
biquadratic spin chain serves as an interesting stress test for the RB-MPS
approach.

To begin with, we consider the entire parameter domain of Fig. 5,
corresponding to $\mathbb{P}_{\mathrm{BLBQ}}=[-\pi,\pi]\times[-2,3]$, covering
all possible phases. Here, the training grid consists of $100\times 100$
uniformly spaced parameter points. Since this gives rise to a significant
linear independence between snapshots, we restrict ourselves to a small system
of $L=24$ spins. Again an energy convergence tolerance of $10^{-9}$ is used,
while the singular value cutoff is reduced to $\mathtt{cut}_{\sigma}=10^{-5}$,
due to increased computational complexity of DMRG in the critical phases. At
this relatively low precision, we face the problem that the induced MPS errors
prohibit the greedy algorithm to resolve the phase diagram up to the desired
accuracy — we merely reach a maximal residual of
$\max_{\bm{\mu}\in\Xi_{\mathrm{train}}}\mathrm{Res}(\bm{\mu})=0.2$ at $N=198$
basis snapshots, before terminating the assembly. Hence the results obtained
from the RB at these settings have to be interpreted with caution. We follow
up on the interplay of MPS accuracy and RB convergence in Sec. IV.2, and here
first discuss the numerical results. Despite the crude accuracy, we are able
to make several interesting observations.

Figure 6: Residual on $\mathbb{P}_{\mathrm{BLBQ}}$ as well as the $N=198$
selected parameter points for the $L=24$ bilinear-biquadratic model. The
critical A phase at central charge $c=2$ accumulates the bulk of the parameter
points, whereas the FM domain is emulated by only one snapshot, with the
surrogate obtaining a low residual across the FM phase.

We begin by illustrating the greedy parameter selection, where the snapshot
parameters and residual on $\mathbb{P}_{\mathrm{BLBQ}}$ are shown in Fig. 6.
It can be observed that the bulk of snapshots concentrates around the critical
A phase and the $c=2$ region features plateau-like structures in the residual.
Moreover, a less concentrated clustering can be seen around the
$\mathrm{SU}(3)$ point [60] at $\bm{\mu}=(-3/4\pi,0)$ and the entire
ferromagnetic phase is spanned by only one snapshot.

Figure 7: Correlation functions on the entire parameter domain
$\mathbb{P}_{\mathrm{BLBQ}}$. Since we consider only a small $L=24$ spin
chain, we observe finite-size effects such as shifted and blurred phase
boundaries, e.g., in the dimerized phase. The correlators are measured at
sites $r=6$ and $r^{\prime}=19$, and the dimer order parameter is computed in
the chain center. Note that the color bar scale is truncated at 1, which is
slightly surpassed by the dimer order parameter and quadrupolar correlator.

We proceed and investigate the gapped domains in parameter space. As in the
previous application, the $z$-Néel and large-$D$ phases are again indicated by
the $z$-spin and $z^{2}$-quadrupolar correlation functions, respectively, and
the Haldane phase is captured by the appropriate string order parameter
$\langle\tilde{O}_{rr^{\prime}}\rangle_{\mathrm{rb}}$ of Eq. (22). To detect
the dimerized phase, we consider the bond order parameter

$\mathcal{D}_{r}(\bm{\mu})=|H_{r}(\bm{\mu})-H_{r+1}(\bm{\mu})|,$ (24)

where
$H_{r}(\bm{\mu})=\cos(\theta)\bm{S}_{r}\cdot\bm{S}_{r+1}+\sin(\theta)(\bm{S}_{r}\cdot\bm{S}_{r+1})^{2}$,
i.e., the difference of next-neighbor interactions among three adjacent sites.
We present the RB measurements in Fig. 7. The most pronounced signal for
dimerization is obtained along the $D=0$ line. The RB results furthermore
indicate that the phase boundary between the dimerized and the large-$D$
regime is located at finite $D>0$ near the $D=0$ line. This is in accord with
earlier high-precision determinations of this phase boundary line [61, 62, 63,
64, 60, 59, 65, 66].

Next, we consider the critical phases. Previous investigations [60, 59]
pointed out that these exhibit dominant (algebraic) correlations of
quadrupolar spin-nematic operators. In order to systematically investigate
these correlations, we measured various quadrupolar structure factors

$\langle
Q^{\gamma}(-k)\,Q^{\gamma}(k)\rangle_{\mathrm{rb}}=\frac{1}{L}\sum_{r,r^{\prime}=1}^{L}\\!\exp[-ik(r-r^{\prime})]\,\langle
Q^{\gamma}_{r}Q^{\gamma}_{r^{\prime}}\rangle_{\mathrm{rb}},$ (25)

using the $L=24$ surrogate, with $\gamma\in\\{xy,x^{2}-y^{2},z^{2}\\}$,
probing thus both transverse and longitudinal quadrupolar correlations, where
$Q^{xy}_{r}=S^{x}_{r}S^{y}_{r}+S^{y}_{r}S^{x}_{r}$,
$Q^{x^{2}-y^{2}}_{r}=(S^{x}_{r})^{2}-(S^{y}_{r})^{2}$, and $Q^{z^{2}}_{r}$ was
given above. Note that the case $\gamma=x^{2}-y^{2}$ is in fact equivalent to
the $xy$ case, due to the residual $\mathrm{U}(1)_{xy}$ symmetry of the
Hamiltonian $H_{\mathrm{BLBQ}}$. Hence, we do not discuss this case separately
here.

Figure 8: Quadrupolar structure factor $\langle
Q^{xy}(-k)\,Q^{xy}(k)\rangle_{\mathrm{rb}}$, shown for $k=0$ in the critical B
domain (left) and at $k=\pi$ in the $c=1$ part of the critical A phase
(right).

The results for $\gamma=xy$ are shown in Fig. 8 in the vicinity of both the
critical A and B phases. From Fig. 8 we detect enhanced transverse
ferroquadrupolar correlations, i.e., at $k=0$, in the regime of the critical B
phase, as expected [60, 59]. Moreover, dominant transverse
antiferroquadrupolar correlations with $k=\pi$ are observed in the $c=1$
domain of critical A. This is in accord with the overall expectation for
enhanced nematic correlations in this regime [59]. The antiferroquadrupolar
character has however not been reported in Ref. [59]. Here, this additional
information is directly available from the RB-MPS approach.

For the $\gamma=z^{2}$ quadrupolar structure factor we observe an enhanced
signal in the $c=2$ domain of the critical A phase, as shown in Fig. 9. In
particular, along the $D=0$ line, the dominant wave vector equals $k=2\pi/3$,
as shown in the left panel of Fig. 9, in accord with earlier reports [64]. For
finite values of $D$ however, the signal at $k=2\pi/3$ reduces, and a plateau
structure is observed. The reduced signal could result from either an overall
suppression of the quadrupolar correlations or from a shift of the dominant
wave vector away from its $D=0$ value of $k=2\pi/3$. The emergence of the
plateau structure indeed already hints towards the latter scenario, as
modulations of the spin structure that are not commensurate with the chain
length of $L=24$ result in the pinning of the best-matching quadrupolar
structure over a finite parameter regime.

Figure 9: Quadrupolar structure factor $\langle
Q^{z^{2}}(-k)\,Q^{z^{2}}(k)\rangle_{\mathrm{rb}}$ on the critical A domain at
$k=2\pi/3$ (left) and on the parameter cut $\\{0.3\pi\\}\times[-2,3]$ (right),
which is indicated by the green dashed line. By varying the wave vector,
different plateaus are highlighted in the $c=2$ region of the critical phase.

In order to further investigate this effect, we narrow down the parameter
domain to the one-dimensional parameter cut $\\{0.3\pi\\}\times[-2,3]$, which
crosses the critical A phase, and on which we generate a second RB for $L=24$
spins. On this smaller parameter domain we reach a maximal residual of
$10^{-2}$ at $N=42$ snapshots, using a higher MPS accuracy of
$\mathtt{cut}_{\sigma}=10^{-6}$. Hence, the RB-MPS method this time requires
significantly fewer snapshots and simultaneously is able to rectify previous
inaccuracies. Again, we evaluate the $z^{2}$-quadrupolar structure factor,
which is shown for different wave vectors in the right panel of Fig. 9. It can
be observed that the dominant wave vector $k$ for the quadrupolar correlations
varies across the $c=2$ regime of the critical A phase, suggesting a
continuously varying $k$ to emerge in the thermodynamic limit. We are not
aware that these incommensurate quadrupolar correlations have been previously
reported for the Hamiltonian $H_{\mathrm{BLBQ}}$. On the other hand, the
presence of such incommensurate correlations provides a simple explanation for
the particular abundance of snapshots picked by the greedy algorithm
throughout the full $c=2$ domain of the critical A phase. Note that while a
relatively large number of snapshots is required to build the RB on this one-
dimensional parameter cut, a comparative scan of truth solves would be
significantly more expensive since a fine resolution is needed to detect the
sharp jumps between plateaus.

From this particular application it becomes apparent, that one of the
strengths of the RB-MPS method lies in being able to scan large domains and
cheaply vary additional parameters, such as the wave vector, for many
different observables with minimal overhead. This enables us to gain
additional insight during a post-processing stage by considering various order
parameters on the fly, which would be significantly more involved and
computationally expensive by way of truth solving only. Furthermore, it proves
viable to generate surrogates on smaller parameter domains, once a large-scale
scan has been performed, in order to improve accuracy and lower the required
RB dimension. Such a procedure of partitioning the parameter space and
respectively generating independent RBs could be systematically implemented
using the techniques from Refs. [67, 68, 69]. Partitioning approaches could
also allow treating larger system sizes for Hamiltonians, such as
$H_{\mathrm{BLBQ}}$, that exhibit complex phase diagrams.

##  IV Technical discussion

###  IV.1 Accuracy and convergence

We move on to assess the accuracy of the RB-MPS method. Similarly to Ref.
[10], we quantify the accuracy of RB surrogates by considering the differences
between truth and RB measurements and maximizing over the parameter domain
under consideration, to obtain the most conservative error estimates. In order
to sample the regions of the parameter domain that have not been selected by
the greedy algorithm and truth solved, we introduce a test grid
$\Xi_{\mathrm{test}}$ that has no mutual points with the training grid. In
practice, the test grid is obtained by shifting the training grid by half a
grid spacing in each parameter direction.

Taking these considerations into account, we estimate the eigenvalue error by

$\mathtt{err}_{\mathrm{val}}=\max_{\bm{\mu}\in\Xi_{\mathrm{test}}}\frac{|E(\bm{\mu})-E_{\mathrm{rb}}(\bm{\mu})|}{|E(\bm{\mu})|},$
(26)

i.e., the maximal relative difference between truth $E(\bm{\mu})$ and RB
energies $E_{\mathrm{rb}}(\bm{\mu})$. For errors of observables, we consider
absolute differences

$\mathtt{err}_{\mathrm{obs}}=\max_{\bm{\mu}\in\Xi_{\mathrm{test}}}|\langle
O(\bm{\mu})\rangle-\langle O(\bm{\mu})\rangle_{\mathrm{rb}}|,$ (27)

since for the observables that we will treat in the following — spin-spin
correlators and quadrupolar structure factors — the measured values become
exactly zero on certain subdomains. Note that the absolute errors in these
cases provide comparable error quantifiers alongside the relative errors
because the range of observable values concentrates around $\langle
O(\bm{\mu})\rangle_{\mathrm{rb}}\sim 0.1$ and $\sim 1$. The vector error
estimate is more subtle since it has to account for different global phases
between truth and RB solutions. We avoid such phase issues by computing the
difference of outer products

$\delta_{\mathrm{vec}}(\bm{\mu})=\frac{\lVert\bm{\Psi}(\bm{\mu})\bm{\Psi}^{\dagger}(\bm{\mu})-\bm{\Phi}_{\mathrm{rb}}(\bm{\mu})\bm{\Phi}^{\dagger}_{\mathrm{rb}}(\bm{\mu})\rVert_{F}}{\lVert\bm{\Psi}(\bm{\mu})\bm{\Psi}^{\dagger}(\bm{\mu})\rVert_{F}},$
(28)

where we denote
$\bm{\Psi}(\bm{\mu})\equiv(\ket{\Psi^{[1]}(\bm{\mu})},\dots,\ket{\Psi^{[m]}(\bm{\mu})})$
to declutter notation, and again maximize over the test grid,
$\mathtt{err}_{\mathrm{vec}}=\max_{\bm{\mu}\in\Xi_{\mathrm{test}}}\delta_{\mathrm{vec}}(\bm{\mu})$.
In this particular form, we would need to explicitly reconstruct Hilbert space
dimensional vectors from MPSs, which is exponentially hard in the system size.
This can be circumvented by computing the norm using the Frobenius inner
product and the cyclic property of the trace, producing the expression

$\delta_{\mathrm{vec}}(\bm{\mu})=\sqrt{1+\frac{\lVert\bm{\Phi}_{\mathrm{rb}}^{\dagger}(\bm{\mu})\bm{\Phi}_{\mathrm{rb}}(\bm{\mu})\rVert^{2}_{F}-2\lVert\bm{\Phi}_{\mathrm{rb}}^{\dagger}(\bm{\mu})\bm{\Psi}(\bm{\mu})\rVert^{2}_{F}}{\lVert\bm{\Psi}^{\dagger}(\bm{\mu})\bm{\Psi}(\bm{\mu})\rVert^{2}_{F}}},$
(29)

where
$\lVert\bm{\Phi}^{\dagger}_{\mathrm{rb}}(\bm{\mu})\bm{\Phi}_{\mathrm{rb}}(\bm{\mu})\rVert_{F}=\lVert\bm{\varphi}^{\dagger}_{\mathrm{rb}}(\bm{\mu})b\bm{\varphi}_{\mathrm{rb}}(\bm{\mu})\rVert_{F}$
and the norm of RB and truth overlaps is found to be

$\displaystyle\lVert\bm{\Phi}_{\mathrm{rb}}^{\dagger}(\bm{\mu})\bm{\Psi}(\bm{\mu})\rVert^{2}_{F}$
$\displaystyle=\sum_{i,j=1}^{m}\big{|}\\!\braket{\Phi_{\mathrm{rb}}^{[i]}(\bm{\mu})}{\Psi^{[j]}(\bm{\mu})}\\!\big{|}^{2}$
$\displaystyle=\sum_{i,j=1}^{m}\Big{|}\sum_{k=1}^{N}[V\varphi_{\mathrm{rb}}^{[i]}(\bm{\mu})]_{k}\braket{\Psi_{k}}{\Psi^{[j]}(\bm{\mu})}\\!\Big{|}^{2}.$
(30)

Thus, we are able to compute the eigenvector error using only RB coefficients
and MPS overlaps at the particular parameter point $\bm{\mu}$, without
resorting to $\mathcal{N}$-dependent operations.

Figure 10: Decay of maximal RB errors for the $L=40$ Haldane chain with
respect to the $\\{0\\}\times[0,2]$ parameter cut.

For the first example, we return to the Haldane chain of $L=40$ spins and
compute the above error quantities on the one-dimensional cut of parameters
$\bm{\mu}\in\\{0\\}\times[0,2]$. Note that we restrict ourselves to a medium
chain length and small parameter cut due to the high computational effort
associated with the required truth solves on $\Xi_{\mathrm{test}}$. At a
singular value cutoff of $\mathtt{cut}_{\sigma}=10^{-8}$ and $N=14$ snapshots,
we converge the surrogate to a maximal residual below $10^{-2}$. For the
observable error, we measure the spin-spin correlator $\langle
S^{y}_{r}S^{y}_{r^{\prime}}\rangle$. The resulting convergence of maximal
errors as the RB is generated is shown in Fig. 10. It is observed, that the
residual decays at a similar rate to the eigenvector error, while the
eigenvalue errors show a significantly faster convergence, as was previously
found for ED-based surrogate models [10]. Moreover, the data demonstrates that
the residual acts as an error surrogate by providing an upper bound to all
real errors — at a residual of $10^{-2}$, the considered error quantities have
decayed to the sub-percent range.

For the next example, we reexamine the $L=24$ bilinear-biquadratic chain on
the $\\{0.3\pi\\}\times[-2,3]$ parameter cut crossing the critical A phase. On
this parameter domain the surrogate has to resolve a large degree of ground-
state variation, and in particular discrete plateaus in the evaluated
observables, which serves as a good example to illustrate possible
difficulties in converging RBs. We adopt the offline settings from Sec. III.2
and again perform truth solves on a shifted test grid. As for the observable
error, we here consider the quadrupolar structure factor $\langle
Q^{z^{2}}(-k)\,Q^{z^{2}}(k)\rangle_{\mathrm{rb}}$ from before. In addition to
maximizing the error over the test grid, we also maximize with respect to all
possible wave vectors $k\in[0,2\pi]$. The maximal and median errors as
functions of the RB dimension are presented in Fig. 11. Most strikingly, we
observe that the maximal RB errors plateau above the maximal residual, which
still decays to $10^{-2}$. While this breaks the desired property of the
residual error estimator of providing an upper bound on all RB errors, we see
that the median errors do converge nicely below the residual line. This
disparity between maximal and median errors is owed to the discrete jumps
observed in the $c=2$ critical A phase as well as the sensitivity of the
Fourier transform in the structure factor to variations of $k$. As shown in
Fig. 12, at one particular point in $\Xi_{\mathrm{test}}$ a truth solve is
performed directly at the transition between two plateaus, where the RB
linearly interpolates between the plateau values instead of producing a sharp
gap, such that we obtain structure factor deviations of order $10^{-1}$.
Leaving this point aside, the DMRG and RB observable measurements are in
excellent agreement, whereas the decay rate of the vector errors is generally
slower.

Figure 11: RB errors for the $L=24$ bilinear-biquadratic model on
$\bm{\mu}\in\\{0.3\pi\\}\times[-2,3]$. The solid lines show the maximal errors
while dotted lines indicate the median errors with respect to
$\Xi_{\mathrm{test}}$. Figure 12: Structure factor $\langle
Q^{z^{2}}(-k)\,Q^{z^{2}}(k)\rangle_{\mathrm{rb}}$ on
$\bm{\mu}\in\\{0.3\pi\\}\times[-2,3]$ computed using the surrogate model
(solid lines) and DMRG truth solves on $\Xi_{\mathrm{test}}$ (markers). The
corresponding deviations between DMRG and RB ground states (vec), structure
factor measurements (obs) and energies (val) are shown below. Note that we
take the maximal structure factor deviations with respect to $k\in[0,2\pi]$.
The maximal observable error, indicated by the red dot, occurs at a
particularly error-sensitive transition between two plateaus in the $c=2$
domain, as shown in the inset, where the maximal deviation is produced by the
jump in the $k=0$ structure factor (encircled in red).

At the chosen MPS accuracies, we thus observe mostly well-behaved convergence
properties similar to those of ED-generated RBs, and the fact that the
surrogate operates based on approximate MPS snapshots and contractions thereof
does not hinder us from obtaining accurate observable measurements. This
prompts the question of when the MPS approximations do become significant and
thereby impede the generation of a meaningful surrogate, for which we turn to
the next section.

###  IV.2 Sensitivity to MPS approximations

One property of RBs which is independent of the specific truth solving method,
is that the surrogate accuracy is ultimately bounded by the truth accuracy.
More specifically, when assembling a RB using the greedy algorithm, we expect
the residual error estimator to stagnate at some value, since finite truth
solving accuracy implies that further snapshots cannot improve the RB
accuracy. In this regard, using DMRG in conjunction with RB methods does not
bring forth anything new since ED solvers are also numerically approximate,
albeit to higher accuracy. What is new, however, is that contractions of
snapshots with operators $\langle\Psi_{i}|A|\Psi_{j}\rangle$ are approximate
because applications of MPOs multiplicatively increase bond dimensions,
therefore requiring further MPS compressions.

To numerically probe the effects of MPS approximations on the greedy RB
assembly, we generate surrogates of the $L=20$ Haldane chain at different
singular value cutoffs on the $\mathbb{P}_{\mathrm{HD}}$ domain, while fixing
all other settings. In particular, we set the residual tolerance and
$\mathtt{cut}_{\lambda}$ to zero, such that the basis assembly could continue
indefinitely, and only terminate when reaching a predefined maximal number of
snapshots. We show the corresponding decay of maximal residuals and minimal
overlap eigenvalues in Fig. 13. Generally it is observed that lower bond
dimension cutoffs allow us to assemble surrogates with lower final residuals,
so to higher overall accuracy. Furthermore, for each greedy assembly we
encounter a point where the residuals suddenly spike and the minimal
eigenvalues drop to zero — which would be unexpected for ED-based surrogates,
where instead the residual is expected to plateau when reaching maximal
surrogate accuracy. The overlap eigenvalues vanish since the greedy algorithm
starts to select close-by or even equal parameter points over the span of the
last iterations, causing a strong increase in linear dependence. This is
accompanied by approximation errors in $h_{q}$, $h_{qq^{\prime}}$ and hence
$\bm{\varphi}_{\mathrm{rb}}$ that eventually prohibit the use of the Rayleigh-
Ritz method as in Eqs. (12, II.1, 14) and result in abrupt increases of the
residual. Note that, due to erroneous coefficients
$\bm{\varphi}_{\mathrm{rb}}$ that occur when, and sometimes right before,
spiking residuals are produced, RB expectation values of observables then may
exhibit unphysical artifacts.

Figure 13: Decay of maximal residuals and minimal normalized overlap
eigenvalues
$\tilde{\lambda}^{\mathrm{min}}=\lVert\Lambda\rVert^{-1}_{F}\min(\Lambda)$ for
different singular value error cutoffs $\mathtt{cut}_{\sigma}$. The red dots
indicate the final chosen dimension of the RB surrogates.

Fortunately, in numerical practice, one can devise simple heuristics to
terminate the greedy assembly at the appropriate time, i.e., right before
approaching a residual spike (see red dots in Fig. 13), in case sufficient MPS
accuracy cannot be guaranteed. Whenever a parameter point is selected twice,
which is forbidden in a correctly assembled greedy basis, the iteration is
stopped. Moreover, a suddenly ill-conditioned RB normalization
$b=B^{\dagger}B$ indicates the drop of minimal overlap eigenvalues, such that
the assembly is stopped when $\lVert I-b\rVert>\delta$ for some threshold
$\delta$ or equivalently when the condition number of $b$ becomes too large to
converge the eigenvalue problem of Eq. (12). We want to remark that even when
too many snapshots are included in the RB, it is possible, without further
computational expense, to roll back the surrogate model by removing the last
snapshots and the corresponding overlaps and matrix elements. For a more
practical discussion of these technical aspects, exemplified by further
physical models, we point the reader to our code documentation available via
the aforementioned repository [40].

##  V Conclusions

In the present work, we expanded the RB method for quantum spin systems with
greedy parameter selection, introduced in Ref. [10], to operate on MPSs
obtained from DMRG solves. MPS approximations allow us to access larger many-
body systems than with ED, and thereby open up the RB approach to more
physically interesting scenarios, while simultaneously providing sufficient
truth accuracy to generate accurate surrogate models. In order to combine the
RB framework with MPS snapshots, we adjusted the orthogonalization method for
the reduced basis $B$, and obtained the required overlaps and matrix elements
for the RB construction from efficient MPS procedures — beyond that the RB-MPS
approach treats the DMRG as a black-box solver (with finite accuracy), thus
profiting from the complementary structure of RB methods. The combined RB-MPS
method was then applied to one-dimensional spin models with rich phase
diagrams: the Haldane spin-1 chain with uniaxial and rhombic-type single-ion
anisotropies as well as the bilinear-biquadratic spin-1 chain with uniaxial
anisotropy. In both applications, the RB-MPS method proved to be a numerically
efficient and versatile tool for determining complex ground-state phase
diagrams.

In particular, the phase diagram of the anisotropic Haldane spin-1 chain was
efficiently obtained from RB surrogates, where the number of required
snapshots grows sublinearly in the system size $L$. Here, the degree of linear
independence across the parameter domain was monitored and controlled by the
eigenvalues of the overlap matrix. The greedy parameter selection mechanism
was found to pick up on physical features of the phase diagram: e.g., the
Haldane phase of $H_{\mathrm{HD}}$ was sampled more frequently, whereas the
ferromagnetic phase in the bilinear-biquadratic chain model
$H_{\mathrm{BLBQ}}$ was spanned by merely one snapshot. Therefore, we actually
gain an automatized indicator for regions of strong ground-state variation.
Furthermore, it was demonstrated that the RB-MPS approach functions on
comparatively large parameter domains containing multiple phases. Overall, RBs
of dimension $N\sim 100$ were sufficient to capture the whole variety of
correlations with fairly high resolution.

Especially in the application to $H_{\mathrm{BLBQ}}$, the RB-MPS approach
qualitatively captured a rather rich phase diagram from large-scale parameter
scans. Furthermore, based on a systematic exploration of quadrupolar structure
factors within the RB-MPS approach, we were able to identify dominant
antiferroquadrupolar correlations within the $c=1$ regime of the critical A
phase, and obtained characteristic signatures of incommensurate spin-nematic
correlations within the $c=2$ domain of the critical A phase. Both these
findings have apparently not been uncovered by previous investigations,
emphasizing the usefulness of a computational tool that efficiently scans
large parameter domains.

In the last part, we analyzed the convergence properties of the RB-MPS
approach. For this purpose, we quantified the maximal RB errors and found
well-behaved error convergence properties for $H_{\mathrm{HD}}$, i.e.,
exponentially decaying RB errors that enter the sub-percent range after $N\sim
10$ snapshots, where the residual constitutes a surrogate error estimator by
providing an upper bound on all RB errors. For the bilinear-biquadratic chain
on the critical A phase, the errors converge mostly well across the parameter
domain despite the rapid and discrete ground-state changes, except for sharp
transitions between structure factor plateaus that are sometimes inaccurately
interpolated by the surrogate models. Lastly, it was explored how the RB
accuracy ultimately depends on the chosen MPS singular value cutoff, where we
found that the greedy algorithm reaches a point of breaking down that can,
however, be avoided by simple heuristics and reversed without additional
computational effort.

The presence of MPS approximations does raise issues that depend on the
specific physics underlying the system on the sampled domain. On very
elaborate parameter domains with highly linearly independent ground states,
i.e., a slowly decaying Kolmogorov $N$-width, the number of required MPS
snapshots may become unwieldy — storing the entire RB in memory hence makes
the greedy assembly particularly memory-intensive. Additionally, we may
encounter, e.g., critical phases that necessitate higher bond dimensions and
thereby make the DMRG solves considerably more expensive. Yet, these issues
can be circumvented by first performing large-scale scans using a crude RB
surrogate and then, informed by the qualitative results, assembling RBs on
smaller parameter subdomains with higher accuracy. In this sense, the RB-MPS
approach offers a particularly convenient workflow for the exploration of
unknown phase diagrams: based on a given surrogate, refined RB models can be
devised in any relevant subregion of the parameter space upon adapting the
offline sampling domain. Moreover, such subregions can be identified based on
readily accessible observables during the online stage at low computational
costs. As one scales up the RB-MPS method to larger system sizes and possibly
higher dimensional parameter spaces, such strategies of incremental
exploration may become increasingly relevant.

Potential future directions of research on the RB-MPS approach include
(automated) grid refinements [67, 68, 69] and parallelization of the offline
stage (parallelizing the online stage is straightforward). Furthermore, we
anticipate generic ground-state probes such as wave-function overlaps [70] and
associated fidelity susceptibilities [71, 72], as well as entanglement
measures [73, 74] and fluctuations of conserved quantities [75, 76, 77, 78],
to be promising further diagnostic tools to combine with the RB-MPS approach
presented here. Finally, we could foresee other tensor network techniques [42,
43, 44, 45] (at least those for finite-size systems) to be naturally
integrated in the RB framework, thus extending the application realm to higher
spatial dimensions, too.

###### Acknowledgements.

We thank Norbert Schuch for helpful suggestions regarding the computation of
eigenvector errors using MPSs. Simulations were performed with computing
resources granted by RWTH Aachen University under project thes1253. M.R.
acknowledges support from the Deutsche Forschungsgemeinschaft (DFG), project
grant 277101999, within the CRC network TR 183 (sub- project B01). P.B. and
S.W. acknowledge support by DFG through RTG 1995.

## References

  * Haldane [1981] F. D. M. Haldane, Ground State Properties of Antiferromagnetic Chains with Unrestricted Spin: Integer Spin Chains as Realisations of the O(3) Non-Linear Sigma Model (1981), ILL preprint SP-81/95, arXiv:1612.00076 .
  * Haldane [1983a] F. D. M. Haldane, Nonlinear Field Theory of Large-Spin Heisenberg Antiferromagnets: Semiclassically Quantized Solitons of the One-Dimensional Easy-Axis Néel State, Phys. Rev. Lett. 50, 1153 (1983a).
  * Haldane [1983b] F. D. M. Haldane, Continuum dynamics of the 1-D Heisenberg antiferromagnet: Identification with the O(3) nonlinear sigma model, Phys. Lett. A 93, 464 (1983b).
  * Gu and Wen [2009] Z.-C. Gu and X.-G. Wen, Tensor-entanglement-filtering renormalization approach and symmetry-protected topological order, Phys. Rev. B 80, 155131 (2009).
  * Pollmann  _et al._ [2012] F. Pollmann, E. Berg, A. M. Turner, and M. Oshikawa, Symmetry protection of topological phases in one-dimensional quantum spin systems, Phys. Rev. B 85, 075125 (2012).
  * Tzeng  _et al._ [2017] Y.-C. Tzeng, H. Onishi, T. Okubo, and Y.-J. Kao, Quantum phase transitions driven by rhombic-type single-ion anisotropy in the $S=1$ Haldane chain, Phys. Rev. B 96, 060404 (2017).
  * Ren  _et al._ [2018] J. Ren, Y. Wang, and W.-L. You, Quantum phase transitions in spin-1 XXZ chains with rhombic single-ion anisotropy, Phys. Rev. A 97, 042318 (2018).
  * Hesthaven  _et al._ [2016] J. S. Hesthaven, G. Rozza, and B. Stamm,  _Certified Reduced Basis Methods for Parametrized Partial Differential Equations_, SpringerBriefs in Mathematics (Springer International Publishing, Cham, 2016).
  * Quarteroni  _et al._ [2016] A. Quarteroni, A. Manzoni, and F. Negri,  _Reduced Basis Methods for Partial Differential Equations_, UNITEXT, Vol. 92 (Springer International Publishing, Cham, 2016).
  * Herbst  _et al._ [2022] M. F. Herbst, B. Stamm, S. Wessel, and M. Rizzi, Surrogate models for quantum spin systems based on reduced-order modeling, Phys. Rev. E 105, 045303 (2022).
  * Aktas and Moses [1998] E. Aktas and F. Moses, Reduced Basis Eigenvalue Solutions for Damaged Structures, Mech. Struct. Mach. 26, 63 (1998).
  * Nair  _et al._ [1998] P. B. Nair, A. J. Keane, and R. S. Langley, Improved First-Order Approximation of Eigenvalues and Eigenvectors, AIAA J. 36, 1721 (1998).
  * Machiels  _et al._ [2000] L. Machiels, Y. Maday, I. B. Oliveira, A. T. Patera, and D. V. Rovas, Output bounds for reduced-basis approximations of symmetric positive definite eigenvalue problems, Comptes Rendus de l’Académie des Sciences - Series I - Mathematics 331, 153 (2000).
  * Fumagalli, Ivan  _et al._ [2016] Fumagalli, Ivan, Manzoni, Andrea, Parolini, Nicola, and Verani, Marco, Reduced basis approximation and a posteriori error estimates for parametrized elliptic eigenvalue problems, ESAIM: Math. Model. Numer. Anal. 50, 1857 (2016).
  * Horger, Thomas  _et al._ [2017] Horger, Thomas, Wohlmuth, Barbara, and Dickopf, Thomas, Simultaneous reduced basis approximation of parameterized elliptic eigenvalue problems, ESAIM: Math. Model. Numer. Anal. 51, 443 (2017).
  * Pla and Herrero [2019] F. Pla and H. Herrero, Reduced basis method applied to eigenvalue problems from convection, Int. J. Bifurc. Chaos Appl. Sci. Eng. 29, 1950028 (2019).
  * Pichi  _et al._ [2020] F. Pichi, A. Quaini, and G. Rozza, A Reduced Order Modeling Technique to Study Bifurcating Phenomena: Application to the Gross–Pitaevskii Equation, SIAM J. Sci. Comp. 42, B1115 (2020).
  * Ekström and Hagen [2019] A. Ekström and G. Hagen, Global Sensitivity Analysis of Bulk Properties of an Atomic Nucleus, Phys. Rev. Lett. 123, 252501 (2019).
  * Demol  _et al._ [2020] P. Demol, T. Duguet, A. Ekström, M. Frosini, K. Hebeler, S. König, D. Lee, A. Schwenk, V. Somà, and A. Tichai, Improved many-body expansions from eigenvector continuation, Phys. Rev. C 101, 041302 (2020).
  * König  _et al._ [2020] S. König, A. Ekström, K. Hebeler, D. Lee, and A. Schwenk, Eigenvector continuation as an efficient and accurate emulator for uncertainty quantification, Phys. Lett. B 810, 135814 (2020).
  * Furnstahl  _et al._ [2020] R. Furnstahl, A. Garcia, P. Millican, and X. Zhang, Efficient emulators for scattering using eigenvector continuation, Phys. Lett. B 809, 135719 (2020).
  * Drischler  _et al._ [2021] C. Drischler, M. Quinonez, P. Giuliani, A. Lovell, and F. Nunes, Toward emulating nuclear reactions using eigenvector continuation, Phys. Lett. B 823, 136777 (2021).
  * Bonilla  _et al._ [2022] E. Bonilla, P. Giuliani, K. Godbey, and D. Lee, Training and projecting: A reduced basis method emulator for many-body physics, Phys. Rev. C 106, 054322 (2022).
  * Frame  _et al._ [2018] D. Frame, R. He, I. Ipsen, D. Lee, D. Lee, and E. Rrapaj, Eigenvector Continuation with Subspace Learning, Phys. Rev. Lett. 121, 032501 (2018).
  * Sarkar and Lee [2021] A. Sarkar and D. Lee, Convergence of Eigenvector Continuation, Phys. Rev. Lett. 126, 032501 (2021).
  * Sarkar and Lee [2022] A. Sarkar and D. Lee, Self-learning emulators and eigenvector continuation, Phys. Rev. Res. 4, 023214 (2022).
  * Melendez  _et al._ [2022] J. A. Melendez, C. Drischler, R. J. Furnstahl, A. J. Garcia, and X. Zhang, Model reduction methods for nuclear emulators, J. Phys. G 49, 102001 (2022).
  * Francis  _et al._ [2022] A. Francis, A. A. Agrawal, J. H. Howard, E. Kökcü, and A. F. Kemper, Subspace diagonalization on quantum computers using eigenvector continuation (2022), arXiv:2209.10571 .
  * Baran and Nichita [2023] V. V. Baran and D. R. Nichita, Reduced basis emulation of pairing in finite systems, Phys. Rev. B 107, 144503 (2023).
  * Franzke  _et al._ [2023] M. C. Franzke, A. Tichai, K. Hebeler, and A. Schwenk, Eigenvector continuation for the pairing Hamiltonian (2023), arXiv:2302.08373 .
  * Mejuto-Zaera and Kemper [2023] C. Mejuto-Zaera and A. F. Kemper, Quantum Eigenvector Continuation for Chemistry Applications (2023), arXiv:2305.00060 .
  * White [1992] S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69, 2863 (1992).
  * White [1993] S. R. White, Density-matrix algorithms for quantum renormalization groups, Phys. Rev. B 48, 10345 (1993).
  * Fannes  _et al._ [1992] M. Fannes, B. Nachtergaele, and R. F. Werner, Finitely Correlated States on Quantum Spin Chains, Commun. Math. Phys. 144, 443 (1992).
  * Klümper  _et al._ [1993] A. Klümper, A. Schadschneider, and J. Zittartz, Matrix Product Ground States for One-Dimensional Spin-1 Quantum Antiferromagnets, EPL 24, 293 (1993).
  * Östlund and Rommer [1995] S. Östlund and S. Rommer, Thermodynamic Limit of Density Matrix Renormalization, Phys. Rev. Lett. 75, 3537 (1995).
  * Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Ann. Phys. 326, 96 (2011).
  * Verstraete  _et al._ [2008] F. Verstraete, V. Murg, and J. Cirac, Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems, Adv. Phys. 57, 143 (2008).
  * Orús [2014] R. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Ann. Phys. (N. Y.) 349, 117 (2014).
  * [40] https://github.com/mfherbst/ReducedBasis.jl.
  * Löwdin [1950] P. Löwdin, On the Non‐Orthogonality Problem Connected with the Use of Atomic Wave Functions in the Theory of Molecules and Crystals, J. Chem. Phys. 18, 365 (1950).
  * Verstraete and Cirac [2004] F. Verstraete and J. I. Cirac, Renormalization algorithms for Quantum-Many Body Systems in two and higher dimensions (2004), arXiv:cond-mat/0407066 .
  * Shi  _et al._ [2006] Y.-Y. Shi, L.-M. Duan, and G. Vidal, Classical simulation of quantum many-body systems with a tree tensor network, Phys. Rev. A 74, 022320 (2006).
  * Gerster  _et al._ [2014] M. Gerster, P. Silvi, M. Rizzi, R. Fazio, T. Calarco, and S. Montangero, Unconstrained tree tensor network: An adaptive gauge picture for enhanced performance, Phys. Rev. B 90, 125154 (2014).
  * Gerster  _et al._ [2017] M. Gerster, M. Rizzi, P. Silvi, M. Dalmonte, and S. Montangero, Fractional quantum hall effect in the interacting Hofstadter model via tensor networks, Phys. Rev. B 96, 195123 (2017).
  * Carleo and Troyer [2017] G. Carleo and M. Troyer, Solving the quantum many-body problem with artificial neural networks, Science 355, 602 (2017).
  * Gubernatis  _et al._ [2016] J. Gubernatis, N. Kawashima, and P. Werner,  _Quantum Monte Carlo Methods: Algorithms for Lattice Models_ (Cambridge University Press, 2016).
  * Bezanson  _et al._ [2017] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, Julia: A Fresh Approach to Numerical Computing, SIAM Rev. 59, 65 (2017).
  * Fishman  _et al._ [2022a] M. Fishman, S. R. White, and E. M. Stoudenmire, The ITensor Software Library for Tensor Network Calculations, SciPost Phys. Codebases , 4 (2022a).
  * Fishman  _et al._ [2022b] M. Fishman, S. R. White, and E. M. Stoudenmire, Codebase release 0.3 for ITensor, SciPost Phys. Codebases , 4 (2022b).
  * Tzeng [2012] Y.-C. Tzeng, Parity quantum numbers in the density matrix renormalization group, Phys. Rev. B 86, 024403 (2012).
  * Affleck  _et al._ [1987] I. Affleck, T. Kennedy, E. H. Lieb, and H. Tasaki, Rigorous results on valence-bond ground states in antiferromagnets, Phys. Rev. Lett. 59, 799 (1987).
  * Note [1] Note that in the ITensor convention, the singular value cutoffs are defined via the squared Frobenius norm, such that the ITensor cutoffs correspond to $\mathtt{cut}_{\sigma}^{2}$.
  * Vidal [2007] G. Vidal, Classical Simulation of Infinite-Size Quantum Lattice Systems in One Spatial Dimension, Phys. Rev. Lett. 98, 070201 (2007).
  * McCulloch [2008] I. P. McCulloch, Infinite size density matrix renormalization group, revisited (2008), arXiv:0804.2509 .
  * Zauner-Stauber  _et al._ [2018] V. Zauner-Stauber, L. Vanderstraeten, M. T. Fishman, F. Verstraete, and J. Haegeman, Variational optimization algorithms for uniform matrix product states, Phys. Rev. B 97, 045145 (2018).
  * Binev  _et al._ [2011] P. Binev, A. Cohen, W. Dahmen, R. DeVore, G. Petrova, and P. Wojtaszczyk, Convergence Rates for Greedy Algorithms in Reduced Basis Methods, SIAM J. Math. Anal. 43, 1457 (2011).
  * den Nijs and Rommelse [1989] M. den Nijs and K. Rommelse, Preroughening transitions in crystal surfaces and valence-bond phases in quantum spin chains, Phys. Rev. B 40, 4709 (1989).
  * De Chiara  _et al._ [2011] G. De Chiara, M. Lewenstein, and A. Sanpera, Bilinear-Biquadratic Spin-1 Chain Undergoing Quadratic Zeeman Effect, Phys. Rev. B 84, 054451 (2011).
  * Rodríguez  _et al._ [2011] K. Rodríguez, A. Argüelles, A. K. Kolezhuk, L. Santos, and T. Vekua, Field-Induced Phase Transitions of Repulsive Spin-1 Bosons in Optical Lattices, Phys. Rev. Lett. 106, 105302 (2011).
  * Fáth and Sólyom [1995] G. Fáth and J. Sólyom, Search for the nondimerized quantum nematic phase in the spin-1 chain, Phys. Rev. B 51, 3620 (1995).
  * Buchta  _et al._ [2005] K. Buchta, G. Fáth, O. Legeza, and J. Sólyom, Probable absence of a quadrupolar spin-nematic phase in the bilinear-biquadratic spin-1 chain, Phys. Rev. B 72, 054433 (2005).
  * Rizzi  _et al._ [2005] M. Rizzi, D. Rossini, G. De Chiara, S. Montangero, and R. Fazio, Phase Diagram of Spin-1 Bosons on One-Dimensional Lattices, Phys. Rev. Lett. 95, 240404 (2005).
  * Läuchli  _et al._ [2006] A. Läuchli, G. Schmid, and S. Trebst, Spin nematics correlations in bilinear-biquadratic ${S}=1$ spin chains, Phys. Rev. B 74, 144426 (2006).
  * Hu  _et al._ [2014] S. Hu, A. M. Turner, K. Penc, and F. Pollmann, Berry-Phase-Induced Dimerization in One-Dimensional Quadrupolar Systems, Phys. Rev. Lett. 113, 027202 (2014).
  * Weyrauch and Rakov [2018] M. Weyrauch and M. V. Rakov, Dimerization in the anisotropic bilinear-biquadratic Heisenberg model, J. Phys. Conf. Ser. 1041, 012004 (2018).
  * Eftang  _et al._ [2010] J. L. Eftang, A. T. Patera, and E. M. Rønquist, An ”$hp$” Certified Reduced Basis Method for Parametrized Elliptic Partial Differential Equations, SIAM J. Sci. Comp. 32, 3170 (2010).
  * Eftang  _et al._ [2011] J. L. Eftang, D. J. Knezevic, and A. T. Patera, An $hp$ certified reduced basis method for parametrized parabolic partial differential equations, Math. Comp. Model. Dyn. Syst. 17, 395 (2011).
  * Eftang and Stamm [2012] J. L. Eftang and B. Stamm, Parameter multi-domain ‘hp’ empirical interpolation, Int. J. Numer. Methods Eng. 90, 412 (2012).
  * Zanardi and Paunković [2006] P. Zanardi and N. Paunković, Ground state overlap and quantum phase transitions, Phys. Rev. E 74, 031123 (2006).
  * You  _et al._ [2007] W.-L. You, Y.-W. Li, and S.-J. Gu, Fidelity, dynamic structure factor, and susceptibility in critical phenomena, Phys. Rev. E 76, 022101 (2007).
  * Gu [2010] S.-J. Gu, Fidelity approach to quantum phase transitions, Int. J. Mod. Phys. B 24, 4371 (2010).
  * Eisert  _et al._ [2010] J. Eisert, M. Cramer, and M. B. Plenio, Colloquium: Area laws for the entanglement entropy, Rev. Mod. Phys. 82, 277 (2010).
  * Zeng  _et al._ [2019] B. Zeng, X. Chen, D.-L. Zhou, and X.-G. Wen,  _Quantum Information Meets Quantum Matter_, Quantum Science and Technology (Springer, New York, NY, 2019).
  * Song  _et al._ [2010] H. F. Song, S. Rachel, and K. Le Hur, General relation between entanglement and fluctuations in one dimension, Phys. Rev. B 82, 012405 (2010).
  * Song  _et al._ [2012] H. F. Song, S. Rachel, C. Flindt, I. Klich, N. Laflorencie, and K. Le Hur, Bipartite fluctuations as a probe of many-body entanglement, Phys. Rev. B 85, 035409 (2012).
  * Rachel  _et al._ [2012] S. Rachel, N. Laflorencie, H. F. Song, and K. Le Hur, Detecting Quantum Critical Points Using Bipartite Fluctuations, Phys. Rev. Lett. 108, 116401 (2012).
  * Contessi  _et al._ [2023] D. Contessi, A. Recati, and M. Rizzi, Phase diagram detection via gaussian fitting of number probability distribution, Phys. Rev. B 107, L121403 (2023).