# Beyond Generalized Eigenvalues in Lattice Quantum Field Theory

George T. Fleming<EMAIL_ADDRESS>Theory Division, Fermi National
Accelerator Laboratory, Batavia, IL

###### Abstract

Two analysis techniques, the generalized eigenvalue method (GEM) or Prony’s
(or related) method (PM), are commonly used to analyze statistical estimates
of correlation functions produced in lattice quantum field theory
calculations. GEM takes full advantage of the matrix structure of correlation
functions but only considers individual pairs of time separations when much
more data exists. PM can be applied to many time separations and many
individual matrix elements simultaneously but does not fully exploit the
matrix structure of the correlation function. We combine both these methods
into a single framework based on matrix polynomials. As these algebraic
methods are well known for producing extensive spectral information about
statistically-noisy data, the method should be paired with some information
criteria, like the recently proposed Bayesean model averaging.

††preprint: FERMILAB-PUB-23-495-T

##  I Introduction

We briefly review the construction of lattice correlation functions, which has
been the subject of many reviews, e.g. [1], and even textbooks. Our review is
not comprehensive and serves only to fix the notation.

A standard computational task of lattice quantum field theory is to generate
statistical samples of Euclidean time correlation functions, which are often
written as vacuum-to-vacuum transition amplitudes

$\displaystyle C_{ab}(\vec{p},|t-t_{0}|)=$ (1)
$\displaystyle\sum_{\vec{x}}e^{i\vec{p}\cdot(\vec{x}-\vec{x}_{0})}\left\langle
0\left|\mathcal{O}_{a}(\vec{x},t)\mathcal{O}_{b}^{\dagger}(\vec{x}_{0},t_{0})\right|0\right\rangle$

where we have made explicit use of translation invariance of the ensemble
average so that it is clear the correlation function is independent of the
spatial source location $\vec{x}_{0}$ and only depends on the relative time
separation of the source and sink $|t-t_{0}|$.
$\mathcal{O}_{b}^{\dagger}(\vec{x}_{0},t_{0})$ is a creation operator that
creates a quantum state with a given set of quantum numbers associated with a
particular irreducible representation of the lattice Hamiltonian and
$\mathcal{O}_{a}(\vec{x},t)$ is an annihilation operator that annihilates a
quantum state of a given set of quantum numbers. If the two sets of quantum
numbers are not commensurate then the correlation is zero by symmetry. There
are many possible operators that can create and annihilate a given set of
quantum numbers, so we give them labels $a,b\in[1,N]$ and the correlation
function can be a $\mathbb{C}^{N\times N}$ matrix in this operator space.

A standard exercise is to expand this correlation function in a complete set
of eigenstates of the lattice Hamiltonian. Associated with each such state is
an energy eigenvalue $E_{k}$ and eigenvector $\left|k\right\rangle$. For
simplicity, we have suppressed labels associated with other conserved quantum
numbers. The correlation function is now written

$C_{ab}(\vec{p},|t-t_{0}|)=\sum_{k}z_{ak}(\vec{p})\ z_{bk}^{*}(\vec{p})\
e^{-E_{k}(\vec{p})\ t}$ (2)

but from here on we will suppress the dependence on momentum $\vec{p}$ and
freely assume time translation invariance. We will assume for each creation
operator $\mathcal{O}_{a}^{\dagger}$ a conjugate annihilation operator can be
defined $\mathcal{O}_{a}$ such that $z_{ak}=\left\langle
0\left|\mathcal{O}_{a}\right|k\right\rangle$ and $z_{ak}^{*}=\left\langle
k\left|\mathcal{O}_{a}^{\dagger}\right|0\right\rangle$. We will also assume
that any operators that are not linearly independent have been removed from
the operator basis and so $z_{ak}$ can be viewed as elements of a matrix with
full rank. In matrix notation, the correlation function can be written
$\bm{C}(t)=\bm{Z}\bm{\Lambda}^{t}\bm{Z}^{\dagger}$. Formally, the total number
of eigenstates should be considered very much larger than the dimension of
$\bm{C}(t)$ that can be computed in a standard lattice calculation. Hence, the
typical data analysis problem is to extract as much reliable information as
possible about the eigenvalue spectrum $E_{k}$ and the matrix elements
$z_{ak}$ given a finite number of samples of the matrix function $\bm{C}(t)$.

Throughout the paper we use the following notation:

$C_{ab}(t)$

Scalar correlation function $\in\mathbb{C}^{1\times 1}$.

$t$

Integer-spaced Euclidean time separations.

$\bm{C}(t)$, $\bm{C}_{t}$

Matrix correlation function $\in\mathbb{C}^{N\times N}$.

$N$

Dimension of operator basis.

$a,b,n$

Operator indices $\in[1,N]$.

$L$

Exact polynomial order (could be $\infty$).

$K$

Assumed polynomial order.

$\ell$

Index $\in[0,L-1]$ or $[0,K-1]$.

$k$

Index of states contributing to $\bm{C}(t)$.

$\lambda_{k}$

$\exp(-E_{k})$ energy contributing to $\bm{C}(t)$.

$\bm{\Lambda}$

Diagonal matrix of energies $\in\mathbb{R}^{LN\times LN}$ or
$\mathbb{R}^{KN\times KN}$.

$z_{ak}$

$\left\langle 0\left|\mathcal{O}_{a}\right|k\right\rangle$ amplitude
contributing to $\bm{C}(t)$.

$\bm{Z}$

Matrix of amplitudes $\in\mathbb{C}^{N\times LN}$ or $\mathbb{C}^{N\times KN}$

$\bm{\mathcal{C}}$

Companion matrix $\in\mathbb{C}^{LN\times LN}$ or $\mathbb{C}^{KN\times KN}$.

$\bm{P}_{\ell}$

Linear prediction matrices $\in\mathbb{C}^{N\times N}$.

$p(M|D)$

Bayesian model probability [2].

$q$

Number of model parameters, usually $LN(N+1)$ or $KN(N+1)$.

###  I.1 Generalized Eigenvalue Method

Early in the development of lattice quantum field theory, Wilson recognized
[3, 4] that correlation function data could be analyzed by truncating the
expansion in Eq. (2) to a relatively small number of states and that the
energies computed in such a fashion would be variational estimates. A simple
method to analyze a scalar correlation function is called effective energy

$\lambda_{\text{eff}}(t,t_{0})^{t-t_{0}}=e^{-E_{\text{eff}}(t,t_{0})\
(t-t_{0})}=\frac{C(t)}{C(t_{0})}$ (3)

Asymptotically,
$E_{\text{eff}}(t,t_{0})=E_{1}+\mathcal{O}(e^{-(E_{2}-E_{1})t})$ as
$t\to\infty$, where $E_{1},E_{2}$ are the lowest and next-to-lowest energy
eigenvalues, assuming the states created by the operators used in the
construction have at least some overlap with the true ground state,
$\left\langle 0\left|\mathcal{O}_{a}\right|1\right\rangle\neq 0$.

The method of using generalized eigenvalues to extract information about the
spectrum of the Hamiltonian is well-studied in lattice QFT [5, 6, 7]. For
simplicity of discussion, we will assume that only $N$ states contribute to
the Hamiltonian and that we have a complete basis of $N$ operators so the
matrix $\bm{Z}$ is of full rank and invertable. Generalized eigenvalues are
computed by solving the following problem

$\lambda_{\text{eff},k}(t,t_{0})^{t-t_{0}}\ \bm{C}(t_{0})\
\bm{v}_{k}(t,t_{0})=\bm{C}(t)\ \bm{v}_{k}(t,t_{0})$ (4)

In our simple case, if $\bm{v}_{k}$ is the $k$-th column vector of
$[\bm{Z}^{\dagger}]^{-1}$ then $\lambda(t,t_{0})=e^{-E_{k}}$ and there will be
$N$ such solutions. Since the general case is that the number of Hamiltonian
eigenstates is much larger than the dimension of the operator basis, the
generalized eigenvalues $\lambda(t,t_{0})$ will only be variational estimates
and will depend on the choice $(t,t_{0})$. Thus, the generalized eigenvalue
method can be viewed as an extension of the effective energy method to matrix-
valued correlation functions.

One approach to bounding the generalized eigenvalues is the eigenvalue
interlacing theorem. Simply stated, if $\bm{A}\in\mathbb{C}^{K\times K}$ is a
Hermitian matrix with eigenvalues
$\alpha_{1}\geq\alpha_{2}\geq\cdots\geq\alpha_{K}$ and
$\bm{B}\in\mathbb{C}^{N\times N},N<K$ is also a Hermitian matrix with
eigenvalues $\beta_{1}\geq\cdots\geq\beta_{N}$ and an orthogonal projection
matrix $\bm{P}\in\mathbb{C}^{N\times K}$ exists such that
$\bm{P}\bm{A}\bm{P}^{\dagger}=\bm{B}$ then the eigenvalues interlace:
$\alpha_{k}\geq\beta_{k}\geq\alpha_{k+K-N}$. To apply these bounds to
generalized eigenvalues requires $\bm{C}(t_{0})$ to be invertable so the
problem can be transformed to a standard eigenvalue problem. Then,
$\lambda_{k}\geq\lambda_{\text{eff},k}(t,t_{0})\geq\lambda_{k+K-N}$ which is
true for any $(t,t_{0})$. Of course, if $K\gg N$ then the right hand bound is
not very restrictive but the more important bound is to the left which
indicates the generalized eigenvalues are always an over-estimate of the true
eigenvalue.

Asymptotic bounds for $\lambda_{\text{eff},k}(t,t_{0})$ as $t\to\infty$ are
also known [7]. For fixed $t_{0}$ as $t\to\infty$ [6]

$\displaystyle E_{\text{eff},k}(t,t_{0})$ $\displaystyle=$ $\displaystyle
E_{k}+\mathcal{O}\left(e^{-\Delta E_{k}t}\right)$ $\displaystyle\Delta E_{k}$
$\displaystyle=$ $\displaystyle\min_{m\neq k}\left|E_{m}-E_{k}\right|.$ (5)

If $t_{0}$ is increased such that $t/2\leq t_{0}$ as $t\to\infty$ then

$E_{\text{eff},k}(t,t_{0})=E_{k}+\mathcal{O}\left(e^{-(E_{N+1}-E_{k})t}\right).$
(6)

It would be desirable to find a method that extends the generalized eigenvalue
method to simultaneously use the values of matrix-valued correlation functions
computed for more than two time separations. The approach we would like to
consider is viewing the generalized eigenvalue method as a matrix polynomial
[8] of linear degree and ask whether we can construct higher degree matrix
polynomials using more time separations. An monic $N\times N$ matrix
polynomial of degree $L$ is

$\bm{\mathcal{P}}_{L}(\lambda)=\bm{I}\lambda^{L}+\sum_{\ell=0}^{L-1}\bm{P}_{\ell}\lambda^{\ell},\
\bm{P}_{\ell}\in\mathbb{C}^{N\times N}$ (7)

Generally, $\det\bm{\mathcal{P}}_{L}(\lambda_{k})=0$ admits $NL$ solutions
$\lambda_{k}$. If the matrix coefficient of $\lambda^{L}$ is not the identity,
as long as it’s invertable, then the monic form can be constructed. For
example, a first-order monic matrix polynomial can be derived from Eq. (4)

$\bm{\mathcal{P}}_{1}(\lambda)=\bm{I}\lambda^{t-t_{0}}-{\bm{Q}^{\dagger}}^{-1}\
\bm{C}(t)\ \bm{Q}^{-1}$ (8)

where $\bm{C}(t_{0})=\bm{Q}^{\dagger}\bm{Q}$ and
$\det\bm{\mathcal{P}}_{1}(\lambda_{k})=0$ has the same $N$ solutions as Eq.
(4).

However, it is not sufficient to construct just any matrix polynomial. For
example, we could add Eq. (4) twice to form
$2\bm{C}(t_{0}+2)\bm{v}=\lambda\bm{C}(t_{0}+1)\bm{v}+\lambda^{2}\bm{C}(t_{0})\bm{v}$,
which would be equivalent to a second-order matrix polynomial. While this
equation yields two solutions in the $N=1$ case, we leave as an exercise for
the reader to see that the solutions are not exact when only two states
contribute to the correlation function [9]. So, our goal is to construct
degree $L$ matrix polynomials from the correlation function data that admit
$NL$ solutions and those solutions are the exact solutions when only $NL$
states contribute. In particular, the matrix polynomials should be equivalent
to the polynomials produced by Prony’s method when $N=1$, since those are
known to produce exact solutions.

###  I.2 Prony’s (and related) Method

Another extension of the effective energy method would be Prony’s method,
first discovered in 1795 [10], and since rediscovered many times in several
related guises. The method extracts multiple effective energies by using a
single correlation function at more than two time separations. The basic idea
is to assume that the correlation function over a range of time separations
can be well-modeled by a spectrum of $L$ effective energies
$\lambda_{k}=\exp(-E_{k})$. This spectrum is used to construct a
characteristic polynomial

$\mathcal{P}_{L}(\lambda)=\prod_{k=1}^{L}(\lambda-\lambda_{k})=\lambda^{L}+\sum_{\ell=0}^{L-1}p_{\ell}\lambda^{\ell}$
(9)

whose linear prediction coefficients $p_{\ell}$ can be computed as functions
of the input spectrum $\lambda_{k}$. $\mathcal{P}_{L}(\lambda_{k})=0$ can be
solved

$\lambda_{k}^{t}=\lambda_{k}^{t-L}\lambda_{k}^{L}=-\lambda_{k}^{t-L}\sum_{\ell=0}^{L-1}p_{\ell}\lambda_{k}^{\ell},\
t\geq L$ (10)

and used in combination with Eq. (2) for a scalar function with a compact
notation, $C_{t}=C(t)$ to write the linear system

$\left[\begin{array}[]{c}C_{L}\\\ C_{L+1}\\\ \vdots\\\
C_{2L-1}\end{array}\right]=-\left[\begin{array}[]{ccc}C_{0}&\cdots&C_{L-1}\\\
C_{1}&\cdots&C_{L}\\\ \vdots&\ddots&\vdots\\\
C_{L-1}&\cdots&C_{2L-2}\end{array}\right]\left[\begin{array}[]{c}p_{0}\\\
p_{1}\\\ \vdots\\\ p_{L-1}\end{array}\right]$ (11)

The structured matrix shown here is commonly called a Hankel matrix even
though Hankel’s work in this area occurred nearly a century later. Prony’s
solution is to take the computed scalar correlation function data $C_{t}$ and
solve the linear system to find the unknown $p_{\ell}$. Then, given the
$p_{\ell}$, the polynomial equation $\mathcal{P}_{L}(\lambda)=0$ can be solved
to find the $\lambda_{k}$.

Once the solution for the $\lambda_{k}$ has been computed, the coefficients
can be determined by solving another structured linear system

$\left[\begin{array}[]{c}C_{0}\\\ C_{1}\\\ C_{2}\\\ \vdots\\\ C_{2L-1}\\\
\end{array}\right]=\left[\begin{array}[]{cccc}1&1&\cdots&1\\\
\lambda_{1}&\lambda_{2}&\cdots&\lambda_{L}\\\
\lambda_{1}^{2}&\lambda_{2}^{2}&\cdots&\lambda_{L}^{2}\\\
\vdots&\vdots&\ddots&\vdots\\\
\lambda_{1}^{2L-1}&\lambda_{2}^{2L-1}&\cdots&\lambda_{L}^{2L-1}\\\
\end{array}\right]\left[\begin{array}[]{c}a_{1}\\\ a_{2}\\\ \vdots\\\
a_{L}\end{array}\right]$ (12)

using singular value decomposition (SVD). The structured matrix here is called
a Vandermonde matrix. Vandermonde, at least, was a contemporary of Prony but
it’s not clear the connection between their work was recognized at the time.
There have been numerous approaches to combining data from multiple matrix
elements but we have yet to see an approach that fully exploits the structure
of the matrices. In particular, as the scalar Prony method involves solving a
polynomial equation to find the spectrum, we anticipate that a fully
structured approach to the matrix problem will involve solving a matrix
polynomial.

In preparation for constructing a block Prony method, we see how a naive
application of linear prediction applied element-by-element to correlator
matrices fails in a case where only two states contribute to the correlation
function

$C_{ab}(t)=\lambda_{1}^{t}z_{a1}z_{b1}^{*}+\lambda_{2}^{t}z_{a2}z_{b2}^{*}$
(13)

Now, use $L=2$ linear prediction

$\displaystyle C_{ab}(2)$ $\displaystyle=$
$\displaystyle-(p_{1}\lambda_{1}+p_{0})z_{a1}z_{b1}^{*}-(p_{1}\lambda_{2}+p_{0})z_{a2}z_{b2}^{*}$
(14) $\displaystyle=$ $\displaystyle-p_{1}C_{ab}(1)-p_{0}C_{ab}(0)$
$\displaystyle\bm{C}_{t}$ $\displaystyle=$
$\displaystyle-\bm{C}_{t-1}(p_{1}\bm{I})-\bm{C}_{t-2}(p_{0}\bm{I})$ (15)

Substituting Eq. (15) into Eq. (11) reveals the problem: solving the block
linear system does not guarantee the resulting linear prediction matrices
$\bm{P}_{0}$, $\bm{P}_{1}$ are proportional to the identity matrix or even
have degenerate eigenvalues.

##  II Block Prony Method

As previously mentioned, there are a number of previous attempts to
incorporate matrix-valued correlation functions into Prony’s method, or other
algebraically-related methods [11, 12, 13, 14]. We will name this method the
Block Prony Method (BPM) as we are not aware of any other method called by
this name. To demonstrate the construction, we will assume we are given $2L$
Hermitian positive definite correlation matrices
$\bm{C}_{t}\in\mathbb{C}^{N\times N}$ computed at equally-spaced time
separations $t\in[0,2L-1]$. We will also assume that precisely $NL$ energies
contribute to the correlation matrices, $\lambda_{k}=\exp(-E_{k})$ and
$k\in[1,NL]$, according to Eq. (2). A matrix
$\bm{\mathcal{C}}\in\mathbb{C}^{NL\times NL}$, called the second companion
matrix, yields the correct characteristic polynomial
$\det[\lambda\bm{I}-\bm{\mathcal{C}}]\propto\prod_{k=1}^{NL}(\lambda-\lambda_{k})$.
In terms of $\mathbb{C}^{N\times N}$ blocks

$\bm{\mathcal{C}}=\left[\begin{array}[]{ccccc}0&0&\cdots&0&-\bm{P}_{0}\\\
\bm{I}&0&\cdots&0&-\bm{P}_{1}\\\ 0&\bm{I}&\ddots&\vdots&-\bm{P}_{2}\\\
\vdots&\ddots&\ddots&0&\vdots\\\
0&\cdots&0&\bm{I}&-\bm{P}_{L-1}\end{array}\right]$ (16)

and the matrix polynomial $\bm{\mathcal{P}}_{L}(\lambda)$ of Eq. (8) using
these blocks satisfies the same characteristic equation. To derive the block
equivalent of Eq. (10), we consider an eigenvector $\bm{v}_{k}$ of
$\bm{\mathcal{C}}$ partitioned into $L$ $\mathbb{C}^{1\times N}$ blocks

$\lambda_{k}\left[\begin{array}[]{c}\bm{v}_{k,0}\\\ \bm{v}_{k,1}\\\ \vdots\\\
\bm{v}_{k,L-1}\end{array}\right]=\left[\begin{array}[]{cccc}0&0&\cdots&-\bm{P}_{0}\\\
\bm{I}&\ddots&\cdots&-\bm{P}_{1}\\\ \vdots&\ddots&0&\vdots\\\
0&\cdots&\bm{I}&-\bm{P}_{L-1}\end{array}\right]\left[\begin{array}[]{c}\bm{v}_{k,0}\\\
\bm{v}_{k,1}\\\ \vdots\\\ \bm{v}_{k,L-1}\end{array}\right]$ (17)

which can also be written

$\displaystyle\lambda_{k}\bm{v}_{k,0}$ $\displaystyle=$
$\displaystyle-\bm{P}_{0}\bm{v}_{k,L-1}$ (18)
$\displaystyle\lambda_{k}\bm{v}_{k,\ell}$ $\displaystyle=$
$\displaystyle\bm{v}_{k,\ell-1}-\bm{P}_{\ell}\bm{v}_{k,L-1}.$ (19)

These equations can be reduced by substitution to

$\lambda_{k}^{L}\bm{v}_{k,L-1}=\left(-\sum_{\ell=0}^{L-1}\lambda_{k}^{\ell}\bm{P}_{\ell}\right)\bm{v}_{k,L-1}$
(20)

Given our assumptions, this can be repeated for any linear combination of
eigenvectors of $\bm{\mathcal{C}}$ which enables us to extend Eq. (11) to
block form

$\left[\begin{array}[]{c}\bm{C}_{L}\\\ \bm{C}_{L+1}\\\ \vdots\\\
\bm{C}_{2L-1}\end{array}\right]=-\left[\begin{array}[]{ccc}\bm{C}_{0}&\cdots&\bm{C}_{L-1}\\\
\bm{C}_{1}&\cdots&\bm{C}_{L}\\\ \vdots&\ddots&\vdots\\\
\bm{C}_{L-1}&\cdots&\bm{C}_{2L-2}\end{array}\right]\left[\begin{array}[]{c}\bm{P}_{0}\\\
\bm{P}_{1}\\\ \vdots\\\ \bm{P}_{L-1}\end{array}\right]$ (21)

Now, the algorithm is clear. Given the $\bm{C}_{t}$, solve Eq. (21) to find
the $\bm{P}_{\ell}$ and then construct the second companion matrix
$\bm{\mathcal{C}}_{L}$ and solve for the eigenvalues $\lambda_{k}$. As an
aside, if more than $2L$ equally-spaced $\bm{C}_{t}$ are available, additional
rows could be appended to the block Hankel in Eq. (21) while keeping the
number of columns fixed. This transforms the problem into an overconstrained
system could be analyzed by least squares using SVD.

##  III Bilinear System (BLS) of Equations

Once the $\lambda_{k}$ have been computed following the block Prony method of
the previous section, it is desirable to use the given $\bm{C}_{t}$ and the
previously computed $\lambda_{t}$ to compute the coefficients $z_{ak}$ of Eq.
(2). In the scalar case, this can be done by solving the Vandermonde system
Eq. (12) but in the matrix case, the coefficients appear in pairs, leading to
a bilinear system (BLS) of equations. A summary of current approaches for
solving BLS problems can be found in [15].

First, let’s count the number of data inputs and the number of free parameters
to see whether such a solution can exist. Given $2L$ Hermitian $\bm{C}(t)$,
there are $2LN(N+1)/2=LN(N+1)$ inputs, $LN$ eigenvalues $\lambda_{k}$ and
$LN^{2}$ matrix elements $z_{nk}$. So the equal number of inputs and free
parameters suggest a solution is possible. If we naively construct the block
Vandermonde system analogous to Eq. (12)

$\left[\begin{array}[]{c}\bm{C}_{0}\\\ \bm{C}_{1}\\\ \vdots\\\
\bm{C}_{2L-1}\\\
\end{array}\right]=\left[\begin{array}[]{ccc}\bm{I}&\cdots&\bm{I}\\\
\lambda_{1}\bm{I}&\cdots&\lambda_{NL}\bm{I}\\\ \vdots&\ddots&\vdots\\\
\lambda_{1}^{2L-1}\bm{I}&\cdots&\lambda_{LN}^{2L-1}\bm{I}\\\
\end{array}\right]\left[\begin{array}[]{c}\bm{A}_{1}\\\ \vdots\\\
\bm{A}_{LN}\end{array}\right]$ (22)

the input data is $\mathbb{C}^{2LN\times N}$, the solution is
$\mathbb{C}^{LN^{2}\times N}$. So, a unique solution exists by SVD if $N\leq
2$ and for $N=2$ the $\bm{A}_{k}$ should be rank-1. In our numerical tests we
have found that the $N=2$ solution is always rank-1, perhaps because our test
function always started from a rank-1 construction, but we don’t have a proof
that it should always be the case. For $N\geq 3$, the system is
underconstrained so an infinite number of general rank solutions exist but it
is still possible to find a unique rank-1 solution [15].

In our particular case, where the number of equations equals the number of
unknowns it should be possible to construct a solution where one unknown is
solved for each matrix element and the result used to eliminate that unknown
from the remaining system of equations by substitution. For example, start
with

$C_{12}(0)=\sum_{k=1}^{LN}z_{1k}z_{2k}^{*}$ (23)

and solve for $z_{11}$

$z_{11}=\frac{1}{z_{21}^{*}}\left[C_{12}(0)-\sum_{k=2}^{LN}z_{1k}z_{2k}^{*}\right]$
(24)

Next, write the equation for $C_{13}(0)$, eliminate $z_{11}$ by substitution

$C_{13}(0)=\frac{z_{31}^{*}}{z_{21}^{*}}\left[C_{12}(0)-\sum_{k=2}^{LN}z_{1k}z_{2k}^{*}\right]+\sum_{k^{\prime}=2}^{LN}z_{1k^{\prime}}z_{3k^{\prime}}$
(25)

and then solve for $z_{12}$, etc. The challenge of this approach is that the
reduced set of equations will become increasingly higher polynomial order in
the remaining unknowns. Still, with a sufficiently small number number of
unknowns, it should be possible for symbolic math program like Mathematica to
find a solution.

Another possibility would be to increase the amount of input data to $LN$
equally spaced $\bm{C}_{t}$ and form an over-constrained Hankel system Eq.
(21) to extract using SVD a least-squares estimate of the same $L$
$\bm{P}_{\ell}$ coefficient matrices. Now, a square Vandermonde system Eq.
(22) can be formed for any $N$ and solved uniquely. There is no guarantee that
the $\bm{A}_{k}$ matrices found this way are rank-1. Also, the amount of input
data exceeds the number of free parameters of a rank-1 solution, so even if a
rank-1 solution is constructed it should be considered, at best, a least-
squares solution. Another consideration is that in Euclidean lattice field
theory calculations it is often computationally less expensive to increase $N$
than it is to increase the number of available times $t$. It is even worse if
one considers the signal-to-noise of such correlation matrices typically
decreases exponentially with increasing time separations [16, 17].

In the absence of an exact method, the method we will advocate is given
$\bm{C}_{t}$ and $\lambda_{k}$ following the block Prony method, use non-
linear least squares minimization to find the solution for the $z_{nk}$. If
the numerical residual is sufficiently close to zero, we will take that as an
indication that we have found the unique rank-1 solution. Setting up the
minimization is not hard since the gradients and Hessian of the residual with
respect to $z_{nk}$ are easy to compute.

##  IV Error-Free Examples

First we will consider an $N=2$ example of $L=6$ polynomial order. We choose
the spectrum to be equally-spaced in the interval $0<\lambda_{k}<1$

$\lambda_{k}=\frac{2L+1-k}{2L+1}\in\left\\{\frac{12}{13},\frac{11}{13},\cdots,\frac{1}{13}\right\\}.$
(26)

The matrix elements $z_{nk}$ are drawn at random in Mathematica using
SeedRandom[11] and Transpose[RandomInteger[{-9,9},{$L$ $N$,$N$}]]

$\bm{Z}=\left[\begin{array}[]{rrrrrrrrrrrrrrrr}-3&-1&6&-3&6&4&4&-6&6&0&-6&6\\\
-3&2&-8&-5&6&5&-6&6&9&-1&-7&4\\\ \end{array}\right]$ (27)

$\det\bm{Z}\bm{Z}^{\dagger}=93894$ confirms the linear independence of the
rows. With these inputs

$\bm{C}_{0}=\left[\begin{array}[]{rr}267&90\\\ 90&382\end{array}\right]$ (28)
$\bm{C}_{1}=\left[\begin{array}[]{rr}\frac{2763}{17}&\frac{437}{17}\\\
\frac{437}{17}&\frac{4040}{17}\end{array}\right]$ (29)

and so forth. The reader can verify with these data that the block Prony
method reproduces the spectrum and the solution to Eq. (22) reproduces the
correct rank-1 solution $\bm{Z}$. As mentioned previously, we do not have a
proof for $N=2$ that the solution will be rank-1 but it is the case in this
example and all other random examples we tried for various sizes $L$.

Next, we consider an $N=3$ example of $L=4$ polynomial order so we can use the
same input spectrum as Eq. (26). The matrix elements $z_{nk}$ are drawn at
random in Mathematica using SeedRandom[13] and
Transpose[RandomInteger[{-9,9},{$L$ $N$,$N$}]]

$\bm{Z}=\left[\begin{array}[]{rrrrrrrrrrrr}5&1&8&4&0&2&-8&6&-4&8&5&-7\\\
6&5&1&-7&4&8&0&-6&1&-7&-2&-5\\\ 6&-1&0&2&8&-5&-6&-4&1&-7&-9&9\\\
\end{array}\right]$ (30)

$\det\bm{Z}\bm{Z}^{\dagger}=38,448,718$ confirms the linear independence of
the rows. With these inputs

$\bm{C}_{0}=\left[\begin{array}[]{rrr}364&-40&-117\\\ -40&306&56\\\
-117&56&394\\\ \end{array}\right]$ (31)
$\bm{C}_{1}=\left[\begin{array}[]{rrr}\frac{2042}{13}&\frac{6}{13}&14\\\
\frac{6}{13}&\frac{2098}{13}&\frac{489}{13}\\\
14&\frac{489}{13}&\frac{1856}{13}\\\ \end{array}\right]$ (32)

and so forth. Again, the reader can verify the block Prony method reproduces
the spectrum. Eq. (22) cannot be used to compute the amplitudes as it is
under-determined. If we define the components of the residual vector from
Eq.(2)

$r_{abt}=C_{ab}(t)-\sum_{k=1}^{LN}z_{zk}z_{bk}^{*}\lambda_{k}^{t}$ (33)

Then minimizing the squared norm of the residual vector

$r^{2}=\sum_{a=1}^{N}\sum_{b\geq a}\sum_{t=0}^{2L-1}|r_{abt}|^{2}$ (34)

using Mathematica’s NMinimize[] with 30 digits of precision finds $r\approx
0.00089$ and

$\bm{Z}\approx\left[\begin{array}[]{rrrrrrrrrrrr}5.0&-0.94&-8.0&-4.2&-0.096&1.8&7.6&-7.1&-2.3&-8.1&5.2&-7.0\\\
6.0&-5.0&-1.1&6.8&4.6&8.0&0.48&5.0&4.2&6.3&-2.2&-5.0\\\
6.0&0.95&0.12&-2.3&7.8&-5.3&6.1&3.0&3.2&6.5&-9.2&9.0\\\ \end{array}\right]$
(35)

whereas computing the residual using the computed eigenvalues and correct
amplitudes from Eq. (32) gives zero to machine precision. Clearly finding the
global minimum of the residual will be a challenging optimization problem.

##  V Variational Analysis

In practical applications, the number of eigenstates exceeds the number of
operators that is typically constructed so any eigenvalues computed using the
block Prony method (BPM) will, at best, be variational estimates. In Sec. I.1
we discussed various bounds and asymptotic limits of generalized eigenvalues
which confirm their variational nature. We would like to test the BPM
eigenvalues in a similar fashion. We continue to consider an error-free
correlation function constructed from 12 states as in Eq. (26) and generate
various combinations of $LN=12$. We then assume the polynomial order is
$K=L-1$ and compute $KN$ BPM eigenvalues for
$t_{\text{min}}=t_{\text{max}}-2K+1$, $t\in[t_{\text{min}},t_{\text{max}}]$.

For any time range, the eigenvalues are always bounded by the interlacing
theorem. They also approach the true values as $t_{\text{max}}\to\infty$. In
Figs. 1, 2, 3 and 4 we compute the exponential decay rate of the error of the
effective energies and compare with the asymptotic estimate in Eq. (6). In all
cases, the decay rate of the error approaches the estimate. In the figures we
plot a vertical line where points to the right satisfy
$t_{\text{min}}>t_{\text{max}}/2$. In [7], this was considered to be a
necessary condition for the asymptotic decay rate to dominate. In the figures,
it seems that this is not always sufficient, particularly if it occurs at
small time values. From this study, we conclude that the BPM eigenvalues are
natural extensions of the generalized eigenvalues to higher polynomial order
and are variational estimates of the true eigenvalues with well-defined
asymptotic behavior.

Figure 1:  Effective energies for $N=2,L=6$ example of Sec. IV, assuming $K=5$
and using $t_{\text{min}}=t_{\text{max}}-2K+1$,
$t\in[t_{\text{min}},t_{\text{max}}]$. Horizontal dashed lines are
$E_{NK+1}-E_{k}$ and points to the right of vertical dotted line satisfy
$t_{\text{min}}>t_{\text{max}}/2$ [7].  Figure 2:  Effective energies for
$N=3,L=4$ example of Sec. IV, assuming $K=3$ and using
$t_{\text{min}}=t_{\text{max}}-2K+1$, $t\in[t_{\text{min}},t_{\text{max}}]$.
Horizontal dashed lines are $E_{NK+1}-E_{k}$ and points to the right of
vertical dotted line satisfy $t_{\text{min}}>t_{\text{max}}/2$ [7].  Figure 3:
Effective energies for $N=4,K=2$ using same example spectrum of Sec. IV and
using $t_{\text{min}}=t_{\text{max}}-2K+1$,
$t\in[t_{\text{min}},t_{\text{max}}]$. Horizontal dashed lines are
$E_{NK+1}-E_{k}$ and points to the right of vertical dotted line satisfy
$t_{\text{min}}>t_{\text{max}}/2$ [7].  Figure 4:  Effective energies for
$N=6,K=1$ using same example spectrum of Sec. IV and using
$t_{\text{min}}=t_{\text{max}}-2K+1$, $t\in[t_{\text{min}},t_{\text{max}}]$.
Horizontal dashed lines are $E_{NK+1}-E_{k}$ and points to the right of
vertical dotted line satisfy $t_{\text{min}}>t_{\text{max}}/2$ [7].

##  VI Comments on Modeling Noisy Data

In practical applications, we must consider how best to apply the block Prony
method (BPM) to model noisy data. In previous studies of this problem [11,
13], because the algebraic method maps all the data into a complete set of
model parameters precisely describing the data, we focused on distinguishing
which model parameters encoded mostly the signal present in the data and which
encoded mostly the noise.

Recently, the method of Bayesian model averaging has been developed [2, 18,
19] for analysis of lattice field theory results and provides us with a tool
for weighing how much information is represented by a given subset of model
parameters. For this discussion, we will define the model probability

$p(M|D)\propto\exp(-\chi^{2}-2q)$ (36)

and the $\chi^{2}$ statistic is computed from residuals of Eq. (33) and some
estimate of the inverse of the data covariance matrix in the usual way. In the
sense described in [18], the full set of BPM parameters describe a “perfect”
model in that $\chi^{2}=0$. Even though the BPM model describes the data
perfectly, it is unlikely maximize $p(M|D)$. To see this, imagine a state $k$
encodes mostly noise. Removing parameters $\lambda_{k}$ and $z_{nk}$ will
reduce $2q$ by $2(N+1)$ but will not increase $\chi^{2}$ by anywhere near the
same amount due to the noise. Thus, the model with state $k$ removed will be
more probable that the full model. So, our heuristic in dealing with noisy
data will be to partition the BPM parameters by states which are mostly signal
and mostly noise by finding the subset of states $k$ which maximizes the model
probability. We will describe this process more completely in a future work.

##  VII Conclusion

In this work, we have noted that the generalized eigenvalue method (GEM) is an
extension of the usual effective energy method to Hermitian matrix-valued
correlation functions and can be reformulated as a first-order matrix
polynomial. We then extend GEM to higher-order matrix polynomials by finding
the generalization of Prony’s method to Hermitian matrix-valued correlation
functions, which we call the block Prony method (BPM). The companion matrix
plays a central role as it enables us to understand the space of solutions of
the matrix polynomial as the eigenspace of the companion matrix. We provided
some error-free numerical examples, including a demonstration that the BPM
eigenvalues are variational estimates with the expected asymptotic error
bounds [7]. We also provided some comments on our future plans for applying
the BPM to the analysis of data with noise.

## References

  * Cohen  _et al._ [2010] S. D. Cohen, G. T. Fleming, and H.-W. Lin, Exponential time series in lattice quantum field theory, in  _Exponential Data Fitting and Its Applications_, edited by V. Pereyra and G. Scherer (Bentham Science Publishers, 2010) Chap. 4, pp. 71–93.
  * Jay and Neil [2021] W. I. Jay and E. T. Neil, Bayesian model averaging for analysis of lattice field theory results, Phys. Rev. D 103, 114502 (2021), arXiv:2008.01069 [stat.ME] .
  * Wilson [1981] K. G. Wilson, Closing remarks (1981), Abingdon Meeting on Lattice Gauge Theories.
  * Berg [1982] B. Berg, Glueball calculations in lattice gauge theories, J. Phys. Colloq. 43, 272 (1982).
  * Michael and Teasdale [1983] C. Michael and I. Teasdale, Extracting Glueball Masses From Lattice QCD, Nucl. Phys. B 215, 433 (1983).
  * Luscher and Wolff [1990] M. Luscher and U. Wolff, How to Calculate the Elastic Scattering Matrix in Two-dimensional Quantum Field Theories by Numerical Simulation, Nucl. Phys. B 339, 222 (1990).
  * Blossier  _et al._ [2009] B. Blossier, M. Della Morte, G. von Hippel, T. Mendes, and R. Sommer, On the generalized eigenvalue method for energies and matrix elements in lattice field theory, JHEP 04, 094, arXiv:0902.1265 [hep-lat] .
  * Gohberg  _et al._ [2009] I. Gohberg, P. Lancaster, and L. Rodman,  _Matrix Polynomials_, Classics in Applied Mathematics (Society for Industrial and Applied Mathematics, Philadelphia, 2009).
  * Fleming [2005] G. T. Fleming, What can lattice QCD theorists learn from NMR spectroscopists?, in  _QCD and Numerical Analysis III_, Lecture Notes in Computational Science and Engineering No. 47, edited by A. Boriçi, A. Frommer, B. Joó, A. Kennedy, and B. Pendleton (Springer-Verlag, 2005) pp. 143–152, hep-lat/0403023 .
  * de Prony [1795] G. C. F. M. R. de Prony, Essai expérimental et analytique sur les lois de la dilatabilité et sur celles de la force expansive de la vapeur de l’eau et de la vapeur de l’alkool, à différentes températures, J. Ecole Poly. 1, 24 (1795), partial translation available in [20].
  * Fleming  _et al._ [2009] G. T. Fleming, S. D. Cohen, H.-W. Lin, and V. Pereyra, Excited-State Effective Masses in Lattice QCD, Phys. Rev. D 80, 074506 (2009), arXiv:0903.2314 [hep-lat] .
  * Beane  _et al._ [2009] S. R. Beane, W. Detmold, T. C. Luu, K. Orginos, A. Parreno, M. J. Savage, A. Torok, and A. Walker-Loud, High Statistics Analysis using Anisotropic Clover Lattices: (I) Single Hadron Correlation Functions, Phys. Rev. D 79, 114502 (2009), arXiv:0903.2990 [hep-lat] .
  * Cushman and Fleming [2020] K. K. Cushman and G. T. Fleming, Automated label flows for excited states of correlation functions in lattice gauge theory, Phys. Rev. E 102, 043303 (2020), arXiv:1912.08205 [hep-lat] .
  * Fischer  _et al._ [2020] M. Fischer, B. Kostrzewa, J. Ostmeyer, K. Ottnad, M. Ueding, and C. Urbach, On the generalised eigenvalue method and its relation to Prony and generalised pencil of function methods, Eur. Phys. J. A 56, 206 (2020), arXiv:2004.10472 [hep-lat] .
  * Johnson  _et al._ [2014] C. R. Johnson, H. Šmigoc, and D. Yang, Solution theory for systems of bilinear equations, Linear and Multilinear Algebra 62, 1553 (2014), arXiv:1303.4988 [math.RA] .
  * Parisi [1984] G. Parisi, The Strategy for Computing the Hadronic Mass Spectrum, Phys. Rept. 103, 203 (1984).
  * Lepage [1989] G. P. Lepage, The Analysis of Algorithms for Lattice Field Theory, in  _From Actions to Answers_ , Proceedings of the Theoretical Advanced Study Institute in Elementary Particle Physics (World Scientific, 1989) https://lib-extopc.kek.jp/preprints/PDF/1990/9003/9003479.pdf.
  * Neil and Sitison [2022] E. T. Neil and J. W. Sitison, Improved information criteria for Bayesian model averaging in lattice field theory, (2022), arXiv:2208.14983 [stat.ME] .
  * Neil and Sitison [2023] E. T. Neil and J. W. Sitison, Model averaging approaches to data subset selection, (2023), arXiv:2305.19417 [stat.ME] .
  * Vandevoorde [1996] D. Vandevoorde,  _A fast exponential decomposition algorithm and its applications to structured matrices_ , Ph.D. thesis, Rensselaer Polytechnic Institute (1996), http://wwwlib.umi.com/dissertations/fullcit/9806112.