soft open fences

# Quantum Maximum Entropy Inference and Hamiltonian Learning

Minbo Gao  Institute of Software, Chinese Academy of Sciences, Beijing, China
University of Chinese Academy of Sciences, Beijing, China  Zhengfeng Ji
Corresponding author<EMAIL_ADDRESS>Department of Computer
Science and Technology,
Tsinghua University, Beijing, China  Fuchao Wei  Department of Computer
Science and Technology,
Tsinghua University, Beijing, China

###### Abstract

Maximum entropy inference and learning of graphical models are pivotal tasks
in learning theory and optimization. This work extends algorithms for these
problems, including generalized iterative scaling (GIS) and gradient descent
(GD), to the quantum realm. While the generalization, known as quantum
iterative scaling (QIS), is straightforward, the key challenge lies in the
non-commutative nature of quantum problem instances, rendering the convergence
rate analysis significantly more challenging than the classical case. Our
principal technical contribution centers on a rigorous analysis of the
convergence rates, involving the establishment of both lower and upper bounds
on the spectral radius of the Jacobian matrix for each iteration of these
algorithms. Furthermore, we explore quasi-Newton methods to enhance the
performance of QIS and GD. Specifically, we propose using Anderson mixing and
the L-BFGS method for QIS and GD, respectively. These quasi-Newton techniques
exhibit remarkable efficiency gains, resulting in orders of magnitude
improvements in performance. As an application, our algorithms provide a
viable approach to designing Hamiltonian learning algorithms.

##  1 Introduction

Maximum entropy inference is a widely used method in machine learning,
particularly in the context of graphical models (McCallum et al., 2000;
Kindermann & Snell, 1980; Ackley et al., 1985; Bresler, 2015; Hamilton et al.,
2017) and natural language processing (Berger et al., 1996). In graphical
models, it is known as the _backward mapping_ , the problem of computing the
model parameters from the marginal information (Wainwright & Jordan, 2007).
The inverse problem of estimating marginal parameters from the model
parameters is called the _forward mapping_.

Maximum entropy inference is also a core concept in statistical physics
(Jaynes, 1957) known as the Jaynes’ principle which links statistical
mechanics and information theory. The Hammersley-Clifford theorem establishes
that, in the classical case, any positive probability distribution satisfying
the local Markov property can be represented as a Gibbs distribution (Lafferty
et al., 2001). Hence, maximum entropy inference is equivalent to learning the
classical Hamiltonian defining the Gibbs distribution (Li & Korb, 2020). In
recent years, sample-efficient and time-efficient algorithms with demonstrated
effectiveness for specific classical graphical models have emerged (Ravikumar
et al., 2010; Bresler, 2015; Hamilton et al., 2017; Klivans & Meka, 2017;
Vuffray et al., 2016). However, it is technically non-trivial to extend these
classical techniques into the quantum domain (Wiebe et al., 2014a, b; Evans et
al., 2019; Wang et al., 2017; Anshu et al., 2021; Amin et al., 2018; Bairey et
al., 2019, 2020; Haah et al., 2022; Bakshi et al., 2023) due to the non-
commutativity of quantum Hamiltonian terms.

In quantum physics, the maximum entropy inference problem naturally arise when
we only have partial local information about a quantum many-body system. Let
$H_{j}$ for $j=1,2,\ldots,m$ be local observables of, say, an $n$-qubit
quantum system. What the state of the system would likely be given only the
local marginal information $\alpha_{j}=\operatorname{tr}(H_{j}\rho)$? The
Jaynes’ principle suggests that we take the state to be the one with the
maximum von Neumann entropy, which is always a Gibbs state of the form
$\xi_{H(\mu)}=\frac{1}{Z(\mu)}e^{-\beta H(\mu)}$ where
$H(\mu)=\sum_{j=1}^{m}\mu_{j}H_{j}$ is a quantum Hamiltonian parameterized by
a vector of real numbers $\mu={(\mu_{j})}_{j=1}^{m}$. That is, given as input
the average values $\alpha={(\alpha_{j})}_{j=1}^{m}$ encoding local
information, the parameter vector $\mu$ of the Hamiltonian is uniquely
determined by $\alpha$ and the computation of $\alpha$ from $\mu$ is an
example of the backward mapping problems.

Maximum entropy problem

$\begin{split}\text{maximize:}\quad&S(\rho)\\\ \text{subject
to:}\quad&\left\langle H_{j},\rho\right\rangle=\alpha_{j},\\\ &\rho\text{\ is
a density matrix},\end{split}$

Dual problem

$\begin{split}\text{minimize:}\quad&\ln\operatorname{tr}\exp\Bigl{(}\sum_{j}\lambda_{j}H_{j}\Bigr{)}-\lambda\cdot\alpha\\\
\text{subject to:}\quad&\lambda_{j}\in.\end{split}$

Figure 1: The maximum entropy problem and its dual problem. Here,
$S(\rho)=-\operatorname{tr}(\rho\ln\rho)$ is the von Neumann entropy of
$\rho$.

Algorithms for solving the maximum entropy inference problem include the
generalized iterative scaling method (Darroch & Ratcliff, 1972; Della Pietra
et al., 1997) which has a rich history of study in statistics and gradient
descent and variable metric methods originated from optimization research
(Malouf, 2002). In this paper, we study quantum generalizations of both the
Generalized Iterative Scaling (GIS) algorithm and quasi-Newton methods for
solving the Hamiltonian inference problem. Overall, the paper offers three key
contributions. Firstly, we present a convergence rate analysis for both
Quantum Iterative Scaling (QIS) and Gradient Descent (GD) algorithms. This is
achieved by representing the iteration’s Jacobian in a concise and explicit
formula, which extends a fundamental result from Liang et al. (2004).
Secondly, we establish bounds on the eigenvalues of the Jacobian, proving
polynomial convergence of both the QIS and GD algorithms. This contribution
represents the main technical novelty of the paper. The non-commutativity of
the Hamiltonian terms makes the convergence analysis highly non-trivial and
divergent from classical proof techniques. Moreover, the lack of a simple
closed-form formula for the derivative of the matrix exponential function, due
to non-commutativity, further complicates the analysis. To address this, we
introduce a newly proposed quantum belief propagation method (Lemma C.3) that
enables the bounding of eigenvalues of the Hessian of the log-partition
function, even in cases where an explicit formula for the Hessian is
unavailable. In classical settings, this portion of the proof would typically
be more straightforward. Lastly, we investigate two variants of quasi-Newton
methods that accelerate the iterative processes for QIS and GD, respectively,
resulting in significant performance improvements. In the following, we
provide a detailed discussion of these three contributions, highlighting their
technical significance and implications.

#### Explicit Jacobian Formula for Iterations.

The QIS algorithm we consider here is a natural quantum counterpart to the GIS
algorithm. When the matrices are all diagonal, the QIS algorithm reduces to
the GIS algorithm for classical maximum entropy inference. Previous work (Ji,
2022) has shown that the algorithm indeed converges to the correct solution of
the maximum entropy problem using the auxiliary function method and matrix
inequalities. However, while the auxiliary function method is versatile, it
falls short in providing a precise assessment of the algorithms’ convergence
rates. To solve this problem, we perform a more tailored analysis of the QIS
algorithm’s convergence rate in this paper. In the classical case, a
corresponding convergence analysis for the GIS algorithm is detailed in Liang
et al. (2004). This classical analysis leverages Ostrowski’s theorem, which
bounds the convergence rate of an iterative procedure by the spectral radius
of the Jacobian matrix associated with the iteration. Our contribution extends
the classical analysis to the quantum setting and addresses the difficulty
that the gradient of the exponential function for matrices is much more
involved than the classical counterpart. We provide a closed-form formula for
the Jacobian of the QIS iteration, generalizing the approach of Liang et al.
(2004). The formula has a concise form $\mathds{1}-P^{-1}L$ where $P$ is a
diagonal matrix whose entries are the mean values of operators related to
$H_{j}$ over the Gibbs state and $L$ is the Hessian of the log-partition
function. This connection to the Hessian matrix $L$ is crucial for the
convergence proof.

The maximum entropy problem has a dual problem, which is an _unconstrained_
optimization problem regarding the log-partition function in Figure 1. As a
comparison, we consider the gradient descent (GD) algorithm for the dual
problem. The Jacobian of the GD update process can be computed as
$\mathds{1}-\eta L$ where $L$ is the Hessian of the log-partition function and
$\eta$ is the step size of GD. In a sense, the Jacobian of the QIS update
rule, $\mathds{1}-P^{-1}L$, can be seen as a mechanism to _adaptively_ choose
the step size for different directions. This is the main advantage that QIS
has over the GD algorithm. Numerical simulations in Figure 2 of Section 6 also
show that QIS converges significantly faster than GD.

#### Upper and Lower Bounds for the Jacobian.

As our main technical contribution, we analyze the eigenvalues of the Jacobian
matrix by establishing both the lower and upper bounds for them.

First, we prove that all eigenvalues of the Jacobian are non-negative. This
result is established by proving an _upper bound_ on the Hessian of the log-
partition function $L\preceq P$. The main difficulty for proving such a bound
arises from the fact that there is no simple explicit formula for the
derivative of the matrix exponential function $\odv*{e^{H+sV}}{s}$ when $H$
and $V$ do not commute. Hasting’s quantum belief propagation (Hastings, 2007)
expresses the derivative $\odv*{e^{H+sV}}{s}$ as the anti-commutator
$\bigl{\\{}e^{H+sV},\Phi(V)\bigr{\\}}$ for some quantum channel $\Phi$
depending on $H+sV$ and is the main technical tool used in many previous works
for addressing this difficulty. However, this form of quantum belief
propagation is not applicable in our case to prove the inequality because the
anti-commutator form only guarantees the Hermitian property of the derivative,
while the inequality requires positivity. We propose a modified quantum belief
propagation operator (see Lemma C.3) to circumvent the problem. In the
modified quantum belief propagation, we express $\odv*{e^{H+sV}}{s}$ as
$e^{(H+sV)/2}\Psi(V)e^{(H+sV)/2}$ for some quantum channel $\Psi(\cdot)$.

Second, we show that the largest eigenvalue of the Jacobian is bounded away
from $1$ by making a connection to the strong positivity of the log-partition
function proven in Anshu et al. (2021), a _lower bound_ on the Hessian of the
log-partition function. The upper and lower bounds of the Jacobian together
complete the convergence rate analysis of the QIS algorithm by using
Ostrowski’s theorem (Theorem B.1). As a corollary, this shows that for local
Hamiltonians satisfying certain technical conditions, the backward mapping
problem can be reduced to the forward mapping problem efficiently. The forward
mapping problem is unfortunately a hard problem in general and one may have to
resort to approximate inference methods such as variational inference
(Wainwright & Jordan, 2007; Cranmer et al., 2019) or Markov entropy
decomposition (Globerson & Jaakkola, 2007; Poulin & Hastings, 2011)

#### Accelerations.

While the QIS and GD algorithms enjoy provable convergence analysis and are
expected to converge in a polynomial number of iterations for local
Hamiltonians at constant temperature, it can exhibit sluggish performance in
practical scenarios. For classical learning of model parameters, quasi-Newton
methods are recommended for solving the maximum entropy inference problems, as
suggested in a systematic comparison of classical maximum entropy inference
algorithms performed in Malouf (2002). Even though the convergence analysis is
usually less established, quasi-Newton methods are usually much faster in
practice than iterative scaling and gradient descent algorithms.

In light of this, we investigate two families of quasi-Newton accelerations of
the algorithms. The first family of heuristic acceleration is based on the
Anderson mixing method (Anderson, 1965). The Anderson mixing algorithm is a
heuristic method for accelerating slow fixed-point iterative algorithms. It
can be seamlessly integrated to work with the QIS algorithm as QIS is indeed a
fixed-point iteration. The Anderson mixing accelerated QIS algorithm (AM-QIS)
has _exactly the same_ computational requirement of the QIS algorithm in terms
of the oracle access and the type of measurements required on the quantum
system. The second family is based on the BFGS method (Nocedal & Wright, 2006;
Yuan, 2015) and in particular the limited memory variant, L-BFGS is applied to
the GD algorithm (L-BFGS-GD). In our numerical simulations, we observed that
AM-QIS and L-BFGS-GD have comparable performance, usually faster by orders of
magnitude than the standard QIS and GD algorithms.

We believe that applying such quasi-Newton heuristics is important for quantum
optimization algorithms like the QIS algorithm considered here. While quantum
computing offers a promising new paradigm with the potential for substantial
speedups in specific problems, the practical construction of large-scale
quantum computers is still in its early stages. Current quantum computing
technology has limitations in terms of scale and suffers from errors. Hence,
quantum computing power remains a scarce and valuable resource. Given this
situation, the careful optimization of resources required to solve problems on
quantum computers emerges as a critical task. The use of Anderson mixing and
BFGS for Hamiltonian inference algorithms and potentially for other fixed-
point iterative quantum algorithms represents an attempt to achieve such
resource optimization. Notably, this approach is not unique to quantum
computing. In fact, the quasi-Newton method, which developed into an important
optimization heuristics, was initially developed by W. Davidon while working
with early classical computers, which often crashed before producing correct
results. In response, he devised faster heuristics to expedite calculations
later known as the first quasi-Newton method! Quantum computers are currently
in its very early stages. They are unstable and prone to errors just like
classical computers in the early days; hence, such heuristic speedups may be
critical for numerical quantum algorithms.

##  2 Preliminary

In this section, we introduce some notations used in this paper. For two real
vectors $x,y\in^{m}$, we define $x\cdot y$ as $\sum_{i=1}^{m}x_{i}y_{i}$. We
sometimes extend this notation to the case when $y$ is a vector of matrices
and write, for example, $\lambda\cdot F$ to mean the summation
$\sum_{j}\lambda_{j}F_{j}$. For matrices $A,B$, define $\left\langle
A,B\right\rangle=\operatorname{tr}(A^{\dagger}B)$. We use $A\succeq B$ or
$B\preceq A$ to mean $A-B$ is a positive semidefinite matrix. A density matrix
$\rho$ is a positive semidefinite matrix of unit trace. The set of density
matrices on Hilbert space $\mathcal{X}$ is denoted $\mathrm{D}(\mathcal{X})$.

Suppose $\mathcal{X}$ is a finite-dimensional Hilbert space and $f$ is a real
convex function. We use $\Delta$ to denote the domain $\mathrm{dom}\,{f}$ of
$f$, the interval on which $f$ takes well-defined finite values. Then $f$
extends to all Hermitian operators in $\mathrm{Herm}_{\Delta}(\mathcal{X})$ as
$f(X)=\sum_{k}f(\lambda_{k})\Pi_{k}$ where $X=\sum_{k}\lambda_{k}\Pi_{k}$ is
the spectral decomposition of $X$. Denote the interior and boundary of
$\Delta$ as $\Delta_{\mathrm{int}}$ and
$\Delta_{\mathrm{bd}}=\Delta\setminus\Delta_{\mathrm{int}}$ respectively. It
is easy to see that the domain of matrix function $f$ is
$\mathrm{Herm}_{\Delta}(\mathcal{X})$, and the interior of the domain is
$\mathrm{Herm}_{\Delta_{\mathrm{int}}}(\mathcal{X})$. For a subset $S$ of
Hermitian matrices, we use $\mathrm{cl}\left(S\right)$ to denote the closure
of $S$.

Given convex function $f$ as above, the Bregman divergence for matrices is
$D_{f}(X,Y)=\operatorname{tr}(f(X)-f(Y)-f^{\prime}(Y)(X-Y))$, where
$X\in\mathrm{Herm}_{\Delta}(\mathcal{X})$ and
$Y\in\mathrm{Herm}_{\Delta_{\mathrm{int}}}(\mathcal{X})$. An important case we
focus on in this paper is $f(x)=x\ln x-x$. In this case, the matrix Bregman
divergence becomes the Kullback-Leibler divergence
$D(X,Y)=\operatorname{tr}(X\ln X-X\ln Y-X+Y)$ defined for non-normalized
matrices $X,Y$. When $X,Y$ are positive semidefinite matrices of trace $1$, it
recovers the Kullback-Leibler relative entropy $D(X,Y)=\operatorname{tr}(X\ln
X-X\ln Y)$. We will need the matrix Bregman-Legendre projection
$\mathcal{L}(Y,\Lambda)$ and Bregman-Legendre conjugate $\ell(Y,\Lambda)$ for
convex function $f(x)=x\ln x-x$ defined as $\mathcal{L}(Y,\Lambda)=\exp(\ln
Y+\Lambda)$, $\ell(Y,\Lambda)=\operatorname{tr}\exp(\ln
Y+\Lambda)-\operatorname{tr}Y$. For $Y\propto\mathds{1}$,
$\ell(Y,\Lambda)=\operatorname{tr}\exp(\Lambda)$ and we omit $Y$ and write it
as $\ell(\Lambda)$.

In this paper, we consider spin Hamiltonians only and write them as
$H=\sum_{j=1}^{m}H_{j}$ where $H_{j}$’s are local terms acting on at most
constant number of neighboring qubits according to certain interaction
geometry. For example, a $\sigma_{z}\otimes\sigma_{z}$ term acting on the
first two qubits is
$\sigma_{z}\otimes\sigma_{z}\otimes\bigl{(}\mathds{1}^{\otimes n-2}\bigr{)}$
for Pauli operator $\sigma_{z}=\begin{pmatrix}1&0\\\ 0&-1\end{pmatrix}$. We
often use $\xi_{H}=\frac{1}{Z}e^{-\beta H}$ to represent the Gibbs state of
the Hamiltonian $H$ for inverse temperature $\beta$ specified in the context
or $\beta=1$ otherwise. Here $Z=\operatorname{tr}e^{-\beta H}$ is the
partition function normalizing the state to have trace $1$ and plays an
important role in statistical physics and also in our work.

##  3 Quantum Iterative Scaling

This section presents a version of the Quantum Iterative Scaling (QIS)
algorithm and discusses its applications in the Hamiltonian inference problem.

We first introduce some notations used in the following discussions. For a
given list of Hermitian matrices $F={(F_{j})}_{j=1}^{m}$, define the linear
family of quantum states $\mathscr{L}(\rho_{0})$ as $\bigl{\\{}\rho\succeq
0\mid\left\langle F_{j},\rho\right\rangle=\left\langle
F_{j},\rho_{0}\right\rangle\bigr{\\}}$. Define the exponential family
$\mathscr{E}(\sigma_{0})$ as
$\bigl{\\{}\frac{1}{Z}\exp(\ln\sigma_{0}+\lambda\cdot F)\bigr{\\}}$. We
introduce the new notation $F_{j}$ playing the role of $H_{j}$ in the previous
discussion as we will need certain normalization conditions. In the end, we
will choose $F_{j}=\frac{H_{j}+\mathds{1}}{2m}$ so $F_{j}$ is a scaled linear
shift of $H_{j}$ such that $F_{j}\succeq 0$ and
$\sum_{j}F_{j}\preceq\mathds{1}$.

We note that, in Ji (2022), the algorithms are designed for non-normalized
matrices and, therefore, there is no need to explicitly normalize $Y^{(t)}$ in
the update. Here, we perform explicit normalization to work with normalized
quantum states and their von Neumann entropy. For
$Y^{(t)}=\exp(\ln\sigma_{0}+\lambda\cdot F)$, the normalization is equivalent
to a linear update in the summation of the exponential function
$Y^{(t)}=\exp(\ln\sigma_{0}+\lambda\cdot F-\ln Z)$ where
$Z=\operatorname{tr}Y^{(t)}$. Hence, we have the following two methods to
handle the normalization. The first is to let the algorithm to find the
normalization implicitly, and this would require that $\mathds{1}$ is in the
span of the $F_{j}$’s. This will be the case if the assumption on $F_{j}$ is
that $\sum_{j}F_{j}=\mathds{1}$. The second is to perform the normalization
explicitly as we did in the algorithm. This approach is advantageous as it
works for all $F_{j}$’s satisfying $\sum_{j}F_{j}\preceq\mathds{1}$ even if
$\mathds{1}$ is not in the span of $F_{j}$’s.

An important special case of the algorithm is when $\sigma_{0}=\mathds{1}/d$
and $D(\rho,\sigma_{0})=\ln(d)-S(\rho)$ where $d$ is the dimension. Then, the
minimization over the linear family is now exactly the maximum entropy problem
as in Figure 1 with $H_{j}=F_{j}$, $\alpha_{j}=\left\langle
F_{j},\rho_{0}\right\rangle$, and $\lambda_{j}=-\beta\mu_{j}$. When all the
operators $F_{j}$’s are diagonal in the computational basis, the QIS algorithm
recovers the GIS algorithm (see e.g. Theorem 5.2 of Csiszár & Shields (2004)).

Require: $\rho_{0},\sigma_{0}\in\mathrm{D}(\mathcal{X})$ such that
$D(\rho_{0},\sigma_{0})<\infty$.
Input: $F=(F_{1},F_{2},\ldots,F_{k})\in{\mathrm{Pos}(\mathcal{X})}^{k}$ and
$\sum_{j=1}^{k}F_{j}\preceq\mathds{1}$.
Output: $\lambda^{(1)},\lambda^{(2)},\cdots$ such that

$\lim_{t\to\infty}D\bigl{(}\rho_{0},\mathcal{L}(\sigma_{0},\lambda^{(t)}\cdot
F)\bigr{)}=\inf_{\lambda\in^{k}}D\bigl{(}\rho_{0},\mathcal{L}(\sigma_{0},\lambda\cdot
F)\bigr{)}.$

1: Initialize $\lambda^{(1)}=(0,0,\ldots,0)$.

2: for $t=1,2,\ldots,$ do

3: Compute $Y^{(t)}=\mathcal{L}(\sigma_{0},\lambda^{(t)}\cdot F)$.

4: for $j=1,2,\ldots,k$ do

5: $\delta^{(t)}_{j}=\ln\left\langle
F_{j},\rho_{0}\right\rangle-\ln\left\langle
F_{j},Y^{(t)}/\operatorname{tr}Y^{(t)}\right\rangle$.

6: end for

7: Update parameters $\lambda^{(t+1)}=\lambda^{(t)}+\delta^{(t)}$.

8: end for

Algorithm 1:  Quantum iterative scaling algorithm.

To give some intuition behind the QIS algorithm, we define
$\xi^{(t)}=Y^{(t)}/\operatorname{tr}Y^{(t)}$ and note that the update in the
QIS algorithm is simply $\delta_{j}=\ln\left\langle
F_{j},\rho_{0}\right\rangle-\ln\left\langle F_{j},\xi^{(t)}\right\rangle$,
which is zero when the linear family constraint $\left\langle
F_{j},\rho\right\rangle=\left\langle F_{j},\rho_{0}\right\rangle$ is satisfied
by $\rho=\xi^{(t)}$. In this case, the algorithm stops updating $\lambda$ in
the $j$-th direction as expected. Otherwise, if the difference between the
current mean value $\left\langle F_{j},\xi^{(t)}\right\rangle$ and the target
value $\left\langle F_{j},\rho_{0}\right\rangle$ is big, so will be the update
$\delta_{j}$. The algorithm is, in this sense, _adaptive_ when compared to
algorithms like multiplicative weight update algorithms (Arora et al., 2012).

The maximum entropy problem has dual program in Figure 1 which is an
unconstrained problem. Hence, it is also attractive to work with the dual
using the gradient descent method (or corresponding quasi-Newton methods
discussed later in the paper). The gradient of the dual objective function is
$\pdv*{\bigl{(}\ln\ell(\lambda\cdot
F)-\lambda\cdot\alpha\bigr{)}}{\lambda_{j}}=\left\langle
F_{j},\xi_{\lambda\cdot F}\right\rangle-\alpha_{j}$, where $\xi_{\lambda\cdot
F}$ is the Gibbs state for Hamiltonian $\lambda\cdot F$. Therefore, in
gradient descent, the update in each step is
$\eta\bigl{(}\alpha_{j}-\left\langle F_{j},\xi_{\lambda\cdot
F}\right\rangle\bigr{)}$, where $\eta$ is the learning rate. This leads to the
gradient descent algorithm in Algorithm 2.

Require: $\rho_{0},\sigma_{0}\in\mathrm{D}(\mathcal{X})$ such that
$D(\rho_{0},\sigma_{0})<\infty$.
Input: $F=(F_{1},F_{2},\ldots,F_{k})\in{\mathrm{Pos}(\mathcal{X})}^{k}$ and
$\sum_{j=1}^{k}F_{j}\preceq\mathds{1}$.
Output: $\lambda^{(1)},\lambda^{(2)},\cdots$ such that

$\lim_{t\to\infty}D\bigl{(}\rho_{0},\mathcal{L}(\sigma_{0},\lambda^{(t)}\cdot
F)\bigr{)}=\inf_{\lambda\in^{k}}D\bigl{(}\rho_{0},\mathcal{L}(\sigma_{0},\lambda\cdot
F)\bigr{)}.$

1: Initialize $\lambda^{(1)}=(0,0,\ldots,0)$.

2: for $t=1,2,\ldots,$ do

3: Compute $Y^{(t)}=\mathcal{L}(\sigma_{0},\lambda^{(t)}\cdot F)$.

4: for $j=1,2,\ldots,k$ do

5: $\delta^{(t)}_{j}=\eta\left\langle
F_{j},\rho_{0}\right\rangle-\eta\left\langle
F_{j},Y^{(t)}/\operatorname{tr}Y^{(t)}\right\rangle$.

6: end for

7: Update parameters $\lambda^{(t+1)}=\lambda^{(t)}+\delta^{(t)}$.

8: end for

Algorithm 2:  Gradient descent algorithm for Kullback-Leibler divergence
minimization.

The QIS algorithm generally outperforms the bare-bones GD algorithm and avoids
the problem of choosing the step size completely. Consider the update of the
QIS algorithm $\ln\alpha_{j}-\ln\left\langle F_{j},\xi_{\lambda\cdot
F}\right\rangle$, and the update of the GD algorithm
$\eta\bigl{(}\alpha_{j}-\left\langle F_{j},\xi_{\lambda\cdot
F}\right\rangle\bigr{)}$, For $\alpha_{j}$ and $\left\langle
F_{j},\xi_{\lambda\cdot F}\right\rangle$ in $\interval[openleft]{0}{1}$, the
QIS update is more aggressive than the dual gradient descent for learning rate
$\eta\leq 1$ while still guarantees the convergence. This effect is more
evident when the two numbers $\alpha_{j}$ and $\left\langle
F_{j},\xi_{\lambda\cdot F}\right\rangle$ are small, which holds in most
applications. We will later see that choosing an appropriate learning rate
will improve the performance of the GD algorithm considered in Section 4, but
it is still less efficient compared to QIS.

##  4 Convergence Rate

In this section, we analyze the geometric convergence rate for the QIS
algorithm in Algorithm 1 for the case when $\sigma_{0}=\mathds{1}/d$. As a
comparison, we will also analyze the convergence rate of the GD algorithm
(Algorithm 2).

We will come across several matrices which are defined here for later
references. For $\lambda$ and $F$, we define the corresponding Gibbs state as
$\xi=\dfrac{\exp(\lambda\cdot F)}{\operatorname{tr}\exp(\lambda\cdot F)}$, and
for an operator $O$, we use $\left\langle
O\right\rangle=\operatorname{tr}(O\xi)$ to mean the average value of $O$ with
respect to $\xi$. Define diagonal matrix $P=\sum_{j}\left\langle
F_{j}\right\rangle\lvert\mspace{0.5mu}j\mspace{0.5mu}\rangle\langle\mspace{0.5mu}j\mspace{0.5mu}\rvert$.
Finally, define $L$ to be the Hessian of the log-partition function
$\ln\operatorname{tr}\exp(\lambda\cdot F)$ with $\lambda$ as the variables.

We have the following two results regarding the QIS and GD algorithms
respectively, proved in Appendix B.

######  Theorem 4.1.

The Jacobian of the iterative update map $\lambda^{(t)}\mapsto\lambda^{(t+1)}$
of Algorithm 1 for $\sigma_{0}=\mathds{1}/d$ is given by $\mathds{1}-P^{-1}L$
for $P$ and $L$ defined above with $\lambda=\lambda^{(t)}$.

######  Theorem 4.2.

The Jacobian of the iterative update map $\lambda^{(t)}\mapsto\lambda^{(t+1)}$
of Algorithm 2 is given by $\mathds{1}-\eta L$ where $L$ is the Hessian of the
log-partition function defined above.

By the Ostrowski theorem stated in Theorem B.1, the geometric convergence rate
of the QIS algorithm is, therefore, governed by the spectral radius of the
Jacobian $\mathds{1}-P^{-1}L$. Hence, we need to prove bounds on the spectral
radius. In Anshu et al. (2021), a non-trivial lower bound on $L$ is proved
(see Theorem C.6), giving an upper bound of the spectral radius. To prove the
lower bound, will need the following result proved in Appendix C.

######  Theorem 4.3.

Let $P$ and $L$ be matrices defined above. We have $L\preceq P$.

We are now able to state the convergence rate of the QIS and GD algorithms.

######  Theorem 4.4.

For Hamiltonian $H=\sum_{j=1}^{m}\mu_{j}H_{j}$ where conditions of Theorem C.6
are satisfied and $H_{j}$ are traceless terms with norm $\left\lVert
H_{j}\right\rVert\leq 1$, the QIS algorithm with
$F_{j}=\frac{\mathds{1}+H_{j}}{2m}$ and the GD algorithm with the same choices
of $F_{j}$ and $\eta=m$ solve the Hamiltonian inference problem with geometric
convergence rate $1-\Omega\bigl{(}\frac{1}{m^{2}}\bigr{)}$.

###### Proof of Theorem 4.4.

We consider the QIS algorithm first, which is the more difficult case. By
Theorem 4.1, $J_{\mathrm{QIS}}=\mathds{1}-P^{-1}L$. Hence, we can bound the
spectral radius as

$\begin{split}r\bigl{(}J_{\mathrm{QIS}}\bigr{)}&=r\bigl{(}\mathds{1}-P^{-1}L\bigr{)}\\\
&=r\bigl{(}\mathds{1}-L^{1/2}P^{-1}L^{1/2}\bigr{)}\\\
&=\left\lVert\mathds{1}-L^{1/2}P^{-1}L^{1/2}\right\rVert\\\ \end{split}$

Here, the second line follows from the fact that $AB$ and $BA$ have the same
set of eigenvalues and the third step follows as $L^{1/2}P^{-1}L^{1/2}$ is
Hermitian. By definition, $P$ is a diagonal matrix whose $(j,j)$-th entry is
$\left\langle
F_{j},\xi\right\rangle=\left\langle\frac{\mathds{1}+H_{j}}{2m},\xi\right\rangle\leq
1/m$. Theorem 4.3 then implies that $\mathds{1}-P^{-1}L$ has eigenvalues in
$\interval{0}{1}$ and

$r\bigl{(}J_{\mathrm{QIS}}\bigr{)}=1-\lambda_{\min}(L^{1/2}P^{-1}L^{1/2})\leq
1-m\lambda_{\min}(L).$ (1)

The Hamiltonian is $H(\mu)=\sum_{j}\mu_{j}(2mF_{j}-\mathds{1})$. Define
$\lambda_{j}=2m\mu_{j}$, and $\tilde{H}(\lambda)=\sum_{j}\lambda_{j}F_{j}$. We
have $\tilde{H}(\lambda)=H(\mu)+\mu_{\Sigma}\mathds{1}$ where
$\mu_{\Sigma}=\sum_{j}\mu_{j}=\sum_{j}\lambda_{j}/(2m)$. We compute

$\begin{split}L_{j,k}&=\pdv{\ln\operatorname{tr}\exp(\lambda\cdot
F)}{\lambda_{j},\lambda_{k}}\\\ &=\pdv{\ln\operatorname{tr}\exp(\lambda\cdot
F)}{\mu_{j},\mu_{k}}\pdv{\mu_{j}}{\lambda_{j}}\pdv{\mu_{k}}{\lambda_{k}}\\\
&=\frac{1}{4m^{2}}\pdv{\bigl{(}\ln\operatorname{tr}\exp(H(\mu))+\sum_{j}\mu_{j}\bigr{)}}{\mu_{j},\mu_{k}}\\\
&=\frac{1}{4m^{2}}\nabla^{2}_{\mu}\ln\operatorname{tr}\exp(H(\mu)).\end{split}$

By Theorem C.6, we have
$\lambda_{\min}(L)\geq\Omega\Bigl{(}\frac{1}{m^{3}}\Bigr{)}$. Together with
Equation 1, this completes the proof using Theorem B.1.

By a similar calculation, we can prove the claim for the GD algorithm. ∎

The above analysis shows that the QIS algorithm has a better geometric
convergence rate even if we set $\eta=m$ in the GD algorithm. Numerical
simulations in Section 6 also confirm this observation. In some sense, the QIS
algorithm is an adaptive gradient descent that can automatically choose the
appropriate learning rate for different dimensions as $\left\langle
F_{j}\right\rangle$ may differ for each $j$.

##  5 Acceleration by Quasi-Newton Methods

The convergence analysis in Section 4 is of theoretical interest but
polynomial convergence proved there is usually not enough for practical
applications. In this section, we explore the application of quasi-Newton
methods, which can significantly improve the efficiency of the adaptive
learning algorithms considered in Section 3. In particular, we study two
families of methods, the Anderson mixing (Anderson, 1965) method and the BFGS
method (Nocedal & Wright, 2006).

Anderson mixing (abbreviated as AM in the following) is a widely used method
employed in numerical and computational mathematics to accelerate the
convergence of fixed-point iterations. It particularly excels in scenarios
where traditional iterative methods may converge slowly or struggle to find
solutions efficiently. The essence of the Anderson mixing algorithm lies in
its ability to dynamically combine and update a finite set of historical
iterates. It adaptively selects a linear combination of these historical
iterates, leveraging the past information to guide the algorithm toward
convergence more effectively. This technique finds applications in various
scientific and engineering domains, including quantum chemistry (Garza &
Scuseria, 2012), machine learning (Sun et al., 2021), and solving complex
systems of equations (Brezinski et al., 2022), where it often delivers
substantial acceleration in computational tasks.

Applying the Anderson mixing method to the Hamiltonian inference problem,
specifically to the Quantum Iterative Scaling (QIS) algorithm, is
straightforward due to the inherent nature of QIS as a fixed-point iterative
update algorithm. The Anderson-accelerated QIS (AM-QIS) algorithm combines
both the QIS iterative step and simple classical processing, so it has exactly
the same requirement as the standard QIS algorithm for the oracle access to
the Gibbs state or the average values. Since the fixed-point map $g(x)$ in QIS
iteration is a contraction, we can set the mixing parameter $\beta_{t}\equiv
1$ defined in Appendix D and the convergence of AM-QIS follows from the
results in Toth & Kelley (2015). We also use the Barzilai-Borwein (BB) method
(Barzilai & Borwein, 1988) for choosing the mixing parameter which turns out
to be effective and provides further accelerations.

The BFGS method and the limited memory variant L-BFGS are the most influential
among many quasi-Newton methods. They are the recommended choice for learning
graphical models in the classical machine learning literature Malouf (2002).
The BFGS method works with an unconstrained optimization problem
$\min_{x\in^{n}}f(x)$. The update in the BFGS algorithm has the form
$x_{k+1}=x_{k}-\eta_{k}H_{k}\nabla f(x_{k})$, where $\eta_{k}$ is the step
size which can usually be found by line search, and $H_{k}$ is a matrix that
is updated iteratively during the execution of the algorithm. We consider both
a fixed choice or the BB method for the initial approximation of the inverse
Hessian $H_{0}$ for a fair comparison with AM. The application of BFGS methods
to our problem is also straightforward as the dual problem is an unconstrained
optimization problem.

Figure 2: Comparison of QIS and GD algorithms. The loss in measured by the
error in the objective function of the maximum entropy problem.

AM and BFGS have different application scenarios. AM is an acceleration method
for solving fixed-point problems and the approximation $G_{t}$ of the inverse
Jacobian matrix is generally not symmetric. In contrast, BFGS is an
optimization method and constructs a symmetric approximation $H_{t}$ for the
inverse Hessian matrix. Since an optimization problem can usually be recast as
a fixed-point problem, AM also applies to solving optimization problems.
However, BFGS may be more efficient in some cases due to the maintained
symmetry structure compared with AM.

##  6 Experiments

We conducted numerical simulations to assess the comparative efficiency of
four approaches: the standard QIS, the standard GD algorithm, AM-QIS, and the
L-BFGS-GD algorithm applied to the dual problem. The Linux workstation we used
for the numerical experiment has a 16-core CPU (Intel(R) Xeon(R) Platinum
8369HB CPU @ 3.30GHz) and 64GB of memory.

In the experimental setup, we adopted a method involving generating random
Gibbs states for random local Hamiltonians, represented as
$H=\sum_{j}\lambda_{j}H_{j}$. Here, the local terms, denoted as $H_{j}$,
consist of tensor products of local Pauli operators, and the $\lambda_{j}$
parameters are the values to be learned. These Hamiltonians were then utilized
to create Gibbs states $\xi_{H}=\frac{1}{Z}e^{-\beta H}$. We feed the Gibbs
states and their corresponding local average values $\alpha_{j}=\left\langle
H_{j}\otimes\mathds{1},\xi_{H}\right\rangle$ to the algorithms. In this way,
we know the ground truth about the values of $\lambda_{j}$’s and the objective
value of the optimization programs in Figure 1, and we choose to evaluate the
algorithms’ performance by the error compared with the true objective value.

The results are summarized in Figures 2 and 3. In Figure 2, we compare the
performance of QIS and GD algorithms. We can see that QIS algorithm is more
efficient than GD algorithm regardless of whether we ensure the completeness
$\sum F_{j}=\mathds{1}$ or not. In Figure 3, we compare the performance of AM-
QIS and L-BFGS-GD, both with and without the Barzilai-Borwein method. We can
see that AM-QIS and L-BFGS-GD are comparable in general. The standard QIS
algorithm typically required approximately 1500 iterations to achieve an error
level of $10^{-6}\sim 10^{-8}$ (measured using the objective function of the
maximum entropy problem). In contrast, the AM-QIS and L-BFGS algorithms
achieved the same accuracy with only about 8 (or 20) iterations with (or
without) BB, showcasing a remarkable speedup of two orders of magnitude. The
efficiency of the AM-QIS algorithm is stable and does not change much when the
Hamiltonian is normalized and completed, while the efficiency of L-BFGS-GD
algorithm (in Figure 3(a)) is observed to be sensitive in this regard.

(a) With BB method. (b) Without BB method.

Figure 3: Comparisons of the AM-QIS and L-BFGS-GD algorithms. The Figure 3(a)
on the left uses the Barzilai-Borwein method to choose the step size and
Figure 3(b) on the right uses fixed step size. The dotted (red) and dashed
(green) lines represent the performance of the algorithms when the input
Hamiltonian terms are complete satisfying $\sum_{j}F_{j}=\mathds{1}$.

##  7 Discussions

In this study, we considered adaptive learning algorithms for the Hamiltonian
inference problem. We examined the convergence of the quantum iterative
scaling algorithm (QIS) and the gradient descent (GD) algorithm for the dual
problem. Furthermore, two quasi-Newton methods AM-QIS and L-BFGS-GD are
proposed.

The QIS algorithm iteratively updates the Hamiltonian parameters adaptively by
comparing $\left\langle H_{j}\otimes\mathds{1},\xi_{H(\lambda)}\right\rangle$
and the target value $\alpha_{j}$. Therefore it requires the ability to
prepare $\xi_{H(\lambda)}$ for the trial parameter $\lambda$ or estimate the
local information of the state. This is generally a computationally demanding
assumption, but if the physical system has exponentially decaying correlation
and satisfy certain Markov property, the preparation of the Gibbs state or its
local observations could be efficient (Brandão & Kastoryano, 2019; Kuwahara et
al., 2020). Furthermore, the issue may be solved or mitigated by combining
quantum belief propagation algorithms proposed in Hastings (2007); Leifer &
Poulin (2008); Poulin & Hastings (2011) which present possible ways of
computing the value $\left\langle
H_{j}\otimes\mathds{1},\xi_{H(\lambda)}\right\rangle$ approximately without
generating the full Gibbs state, thereby removing the use of the adaptive
Gibbs oracle. We leave the exploration of this possibility as future work.

In the proof of the upper bound of the Hessian of the log-partition function,
we developed a modified quantum belief propagation technique, which may be of
independent interest. It is an interesting problem to find more applications
of this new tool.

## References

  * Ackley et al. (1985) Ackley, D. H., Hinton, G. E., and Sejnowski, T. J.  A learning algorithm for Boltzmann machines.  _Cognitive Science_ , 9(1):147–169, 1985.
  * Amin et al. (2018) Amin, M. H., Andriyash, E., Rolfe, J., Kulchytskyy, B., and Melko, R.  Quantum Boltzmann Machine.  _Physical Review X_ , 8(2):021050, 2018.
  * Anderson (1965) Anderson, D. G.  Iterative Procedures for Nonlinear Integral Equations.  _Journal of the ACM_ , 12(4):547–560, 1965.
  * Anshu et al. (2021) Anshu, A., Arunachalam, S., Kuwahara, T., and Soleimanifar, M.  Sample-efficient learning of interacting quantum systems.  _Nature Physics_ , pp. 1–5, 2021.
  * Arora et al. (2012) Arora, S., Hazan, E., and Kale, S.  The Multiplicative Weights Update Method: a Meta-Algorithm and Applications.  _Theory of Computing_ , 8(6):121–164, 2012.
  * Bairey et al. (2019) Bairey, E., Arad, I., and Lindner, N. H.  Learning a Local Hamiltonian from Local Measurements.  _Physical Review Letters_ , 122(2):020504, 2019\.
  * Bairey et al. (2020) Bairey, E., Guo, C., Poletti, D., Lindner, N. H., and Arad, I.  Learning the dynamics of open quantum systems from their steady states.  _New Journal of Physics_ , 22(3):032001, 2020\.
  * Bakshi et al. (2023) Bakshi, A., Liu, A., Moitra, A., and Tang, E.  Learning quantum Hamiltonians at any temperature in polynomial time.  arXiv:2310.02243, 2023.
  * Barzilai & Borwein (1988) Barzilai, J. and Borwein, J. M.  Two-Point Step Size Gradient Methods.  _IMA Journal of Numerical Analysis_ , 8(1):141–148, 1988.
  * Berger et al. (1996) Berger, A. L., Della Pietra, S. A., and Della Pietra, V. J.  A Maximum Entropy Approach to Natural Language Processing.  _Computational Linguistics_ , 22(1):39–71, 1996\.
  * Brandão & Kastoryano (2019) Brandão, F. G. S. L. and Kastoryano, M. J.  Finite Correlation Length Implies Efficient Preparation of Quantum Thermal States.  _Communications in Mathematical Physics_ , 365(1):1–16, 2019.
  * Bresler (2015) Bresler, G.  Efficiently Learning Ising Models on Arbitrary Graphs.  In _Proceedings of the forty-seventh annual ACM symposium on Theory of Computing_ , STOC ’15, pp. 771–782, New York, NY, USA, 2015\. Association for Computing Machinery.
  * Brezinski et al. (2022) Brezinski, C., Cipolla, S., Redivo-Zaglia, M., and Saad, Y.  Shanks and Anderson-type acceleration techniques for systems of nonlinear equations.  _IMA Journal of Numerical Analysis_ , 42(4):3058–3093, 2022.
  * Cranmer et al. (2019) Cranmer, K., Golkar, S., and Pappadopulo, D.  Inferring the quantum density matrix with machine learning, 2019.
  * Csiszár & Shields (2004) Csiszár, I. and Shields, P. C.  Information Theory and Statistics: A Tutorial.  _Foundations and Trends® in Communications and Information Theory_ , 1(4):417–528, 2004.
  * Darroch & Ratcliff (1972) Darroch, J. N. and Ratcliff, D.  Generalized Iterative Scaling for Log-Linear Models.  _Annals of Mathematical Statistics_ , 43(5):1470–1480, 1972.
  * Davidon (1991) Davidon, W. C.  Variable Metric Method for Minimization.  _SIAM Journal on Optimization_ , 1(1):1–17, 1991\.
  * Della Pietra et al. (1997) Della Pietra, S., Della Pietra, V., and Lafferty, J.  Inducing features of random fields.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 19(4):380–393, 1997.
  * Evans et al. (2019) Evans, T. J., Harper, R., and Flammia, S. T.  Scalable Bayesian Hamiltonian learning, 2019.
  * Fang & Saad (2009) Fang, H.-r. and Saad, Y.  Two classes of multisecant methods for nonlinear acceleration.  _Numerical Linear Algebra with Applications_ , 16(3):197–221, 2009.
  * Garza & Scuseria (2012) Garza, A. J. and Scuseria, G. E.  Comparison of self-consistent field convergence acceleration techniques.  _The Journal of Chemical Physics_ , 137(5):054110, 2012.
  * Globerson & Jaakkola (2007) Globerson, A. and Jaakkola, T.  Approximate inference using conditional entropy decompositions.  In _Proceedings of the Eleventh International Conference on Artificial Intelligence and Statistics_ , pp. 131–138. PMLR, 2007.
  * Haah et al. (2022) Haah, J., Kothari, R., and Tang, E.  Optimal learning of quantum Hamiltonians from high-temperature Gibbs states.  In _2022 IEEE 63rd Annual Symposium on Foundations of Computer Science (FOCS)_ , pp. 135–146, 2022.
  * Hamilton et al. (2017) Hamilton, L., Koehler, F., and Moitra, A.  Information Theoretic Properties of Markov Random Fields, and their Algorithmic Applications.  In _Advances in Neural Information Processing Systems_ , volume 30. Curran Associates, Inc., 2017.
  * Hastings (2007) Hastings, M. B.  Quantum belief propagation: An algorithm for thermal quantum systems.  _Physical Review B_ , 76(20):201102, 2007.
  * Jaynes (1957) Jaynes, E. T.  Information Theory and Statistical Mechanics.  _Physical Review_ , 106(4):620–630, 1957.
  * Ji (2022) Ji, Z.  Classical and quantum iterative optimization algorithms based on matrix Legendre-Bregman projections.  arXiv:2209.14185, 2022.
  * Kindermann & Snell (1980) Kindermann, R. and Snell, J. L.  _Markov random fields and their applications_.  Number 1 in Contemporary mathematics. American mathematical society, Providence, 1980.
  * Klivans & Meka (2017) Klivans, A. and Meka, R.  Learning Graphical Models Using Multiplicative Weights.  In _2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS)_ , pp. 343–354, 2017.
  * Kuwahara et al. (2020) Kuwahara, T., Kato, K., and Brandão, F. G. S. L.  Clustering of Conditional Mutual Information for Quantum Gibbs States above a Threshold Temperature.  _Physical Review Letters_ , 124(22):220601, 2020\.
  * Lafferty et al. (2001) Lafferty, J. D., McCallum, A., and Pereira, F. C. N.  Conditional Random Fields: Probabilistic Models for Segmenting and Labeling Sequence Data.  In _Proceedings of the Eighteenth International Conference on Machine Learning_ , ICML ’01, pp. 282–289, San Francisco, CA, USA, 2001\. Morgan Kaufmann Publishers Inc.
  * Leifer & Poulin (2008) Leifer, M. S. and Poulin, D.  Quantum Graphical Models and Belief Propagation.  _Annals of Physics_ , 323(8):1899–1946, 2008\.
  * Li & Korb (2020) Li, Y. and Korb, K. B.  Learning Graphical Models.  In Phung, D., Webb, G. I., and Sammut, C. (eds.), _Encyclopedia of Machine Learning and Data Science_ , pp. 1–10. Springer US, New York, NY, 2020.
  * Liang et al. (2004) Liang, G., Yu, B., and Taft, N.  Maximum entropy models: convergence rates and applications in dynamic system monitoring.  In _International Symposium onInformation Theory, 2004. ISIT 2004. Proceedings._ , pp. 168–174, 2004.
  * Liu & Nocedal (1989) Liu, D. C. and Nocedal, J.  On the limited memory BFGS method for large scale optimization.  _Mathematical Programming_ , 45(1):503–528, 1989\.
  * Malouf (2002) Malouf, R.  A comparison of algorithms for maximum entropy parameter estimation.  In _proceeding of the 6th conference on Natural language learning - COLING-02_ , volume 20, pp. 1–7, 2002.
  * McCallum et al. (2000) McCallum, A., Freitag, D., and Pereira, F. C. N.  Maximum Entropy Markov Models for Information Extraction and Segmentation.  In _Proceedings of the Seventeenth International Conference on Machine Learning_ , ICML ’00, pp. 591–598, San Francisco, CA, USA, 2000. Morgan Kaufmann Publishers Inc.
  * Nocedal & Wright (2006) Nocedal, J. and Wright, S. J.  _Numerical optimization_.  Springer series in operations research. Springer, New York, 2nd ed edition, 2006.
  * Ostrowski (1966) Ostrowski, A. M.  _Solution of Equations and Systems of Equations_.  Academic Press, 2nd edition, 1966.
  * Poulin & Hastings (2011) Poulin, D. and Hastings, M. B.  Markov Entropy Decomposition: A Variational Dual for Quantum Belief Propagation.  _Physical Review Letters_ , 106(8):080403, 2011\.
  * Qiskit (2023) Qiskit, T.  Qiskit: An open-source framework for quantum computing, 2023.
  * Ravikumar et al. (2010) Ravikumar, P., Wainwright, M. J., and Lafferty, J. D.  High-dimensional Ising model selection using l1-regularized logistic regression.  _The Annals of Statistics_ , 38(3):1287–1319, 2010.
  * Sun et al. (2021) Sun, K., Wang, Y., Liu, Y., Zhao, Y., Pan, B., Jui, S., Jiang, B., and Kong, L.  Damped Anderson Mixing for Deep Reinforcement Learning: Acceleration, Convergence, and Stabilization.  In _Advances in Neural Information Processing Systems_ , volume 34, pp. 3732–3743. Curran Associates, Inc., 2021.
  * Toth & Kelley (2015) Toth, A. and Kelley, C. T.  Convergence Analysis for Anderson Acceleration.  _SIAM Journal on Numerical Analysis_ , 53(2):805–819, 2015.
  * Vuffray et al. (2016) Vuffray, M., Misra, S., Lokhov, A., and Chertkov, M.  Interaction Screening: Efficient and Sample-Optimal Learning of Ising Models.  In _Advances in Neural Information Processing Systems_ , volume 29. Curran Associates, Inc., 2016.
  * Wainwright & Jordan (2007) Wainwright, M. J. and Jordan, M. I.  Graphical Models, Exponential Families, and Variational Inference.  _Foundations and Trends® in Machine Learning_ , 1(1–2):1–305, 2007.
  * Wang et al. (2017) Wang, J., Paesani, S., Santagati, R., Knauer, S., Gentile, A. A., Wiebe, N., Petruzzella, M., O’Brien, J. L., Rarity, J. G., Laing, A., and Thompson, M. G.  Experimental quantum Hamiltonian learning.  _Nature Physics_ , 13(6):551–555, 2017.
  * Wiebe et al. (2014a) Wiebe, N., Granade, C., Ferrie, C., and Cory, D.  Quantum Hamiltonian learning using imperfect quantum resources.  _Physical Review A_ , 89(4):042314, 2014a.
  * Wiebe et al. (2014b) Wiebe, N., Granade, C., Ferrie, C., and Cory, D. G.  Hamiltonian Learning and Certification Using Quantum Resources.  _Physical Review Letters_ , 112(19):190501, 2014b.
  * Yuan (2015) Yuan, Y.-X.  Variable Metric Algorithms.  In Engquist, B. (ed.), _Encyclopedia of Applied and Computational Mathematics_ , pp. 1515–1519. Springer, Berlin, Heidelberg, 2015.

##  Appendix A Prime and Dual Divergence Minimization Problems

The primal and dual formulation of the maximum entropy program in Figure 1 is
a special case of the following duality result between two minimization
problems of the Kullback-Leibler divergence given in Figure 4. The duality
theorem of Ji (2022), or Jaynes’ principle, states that the following two
problems have the same minimizer which is the unique intersection point of the
linear family $\mathscr{L}(\rho_{0})$ and exponential family
$\mathrm{cl}\left(\mathscr{E}(\sigma_{0})\right)$ defined in Section 3. When
$\sigma_{0}$ is the maximally mixed state $\mathds{1}/d$, the linear family
minimization is the maximum entropy problem and the exponential family
minimization is a dual program in Figure 1.

Linear family minimization

minimize: $\displaystyle D(X,\sigma_{0})$ subject to: $\displaystyle
X\in\mathscr{L}(\rho_{0}).$

Exponential family minimization

minimize: $\displaystyle D(\rho_{0},Y)$ subject to: $\displaystyle
Y\in\mathscr{E}(\sigma_{0}).$

Figure 4: Two optimization problems of the Kullback-Leibler divergence that
are dual to each other.

##  Appendix B Proofs for Convergence Rate

This section proves the explicit formulas for the Jacobian matrix of the
iterations in QIS and GD algorithms.

###### Proof of Theorem 4.1.

We first prove two identities about the partial derivative of the function
$\ell$ and its natural logarithm.

$\displaystyle\pdv*{\ell(\lambda\cdot F)}{\lambda_{j^{\prime}}}$
$\displaystyle=\left\langle F_{j^{\prime}},\exp(\lambda\cdot F)\right\rangle,$
(2) $\displaystyle\pdv*{\ln\ell(\lambda\cdot F)}{\lambda_{j^{\prime}}}$
$\displaystyle=\left\langle F_{j^{\prime}},\xi\right\rangle,$ (3)

for $\xi$ defined as

$\xi=\dfrac{\exp(\lambda\cdot F)}{\operatorname{tr}\exp(\lambda\cdot F)}.$ (4)

In fact, we have

$\begin{split}\pdv*{\ell(\lambda\cdot
F)}{\lambda_{j^{\prime}}}&=\pdv*{\operatorname{tr}\,\exp(\lambda\cdot
F)}{\lambda_{j^{\prime}}}\\\
&=\operatorname{tr}\,\sum_{k=0}^{\infty}\pdv*{\frac{{(\lambda\cdot
F)}^{k}}{k!}}{\lambda_{j^{\prime}}}\\\
&=\operatorname{tr}\,\sum_{k=1}^{\infty}\sum_{j=0}^{k-1}\frac{{(\lambda\cdot
F)}^{k-j-1}F_{j^{\prime}}{(\lambda\cdot F)}^{j}}{k!}\\\
&=\Bigl{\langle}F_{j^{\prime}},\sum_{k=1}^{\infty}\frac{{(\lambda\cdot
F)}^{k-1}}{(k-1)!}\Bigr{\rangle}\\\ &=\left\langle
F_{j^{\prime}},\exp(\lambda\cdot F)\right\rangle,\end{split}$

where the second line follows from the Taylor expansion of the matrix
expansion and the fourth line is by the cyclic property of trace. This proves
Equation 2. Similarly, we have

$\begin{split}\pdv*{\ln\ell(\lambda\cdot
F)}{\lambda_{j^{\prime}}}&=\frac{1}{\ell(\lambda\cdot
F)}\pdv*{\ell(\lambda\cdot F)}{\lambda_{j^{\prime}}}\\\
&=\frac{1}{\operatorname{tr}\exp(\lambda\cdot F)}\left\langle
F_{j^{\prime}},\exp(\lambda\cdot F)\right\rangle\\\ &=\left\langle
F_{j^{\prime}},\xi\right\rangle,\end{split}$

where the second line follows from Equation 2. This completes the proof of
Equation 3.

Now, for $j=1,2,\ldots,k$, the update in the algorithm is

$\delta_{j}=\ln\left\langle F_{j},\rho_{0}\right\rangle-\ln\left\langle
F_{j},Y^{(t)}/\operatorname{tr}Y^{(t)}\right\rangle.$

Hence, for $j,j^{\prime}\in\\{1,2,\ldots,k\\}$ and $\xi$ defined in Equation
4,

$\begin{split}\pdv{\delta_{j}}{\lambda_{j^{\prime}}}=&-\frac{1}{\left\langle
F_{j},\xi\right\rangle}\pdv{\left\langle
F_{j},\xi\right\rangle}{\lambda_{j^{\prime}}}\\\ =&-\frac{1}{\left\langle
F_{j},\xi\right\rangle}\pdv*{\pdv*{\ln\ell(\lambda\cdot
F)}{\lambda_{j}}}{\lambda_{j^{\prime}}}\\\
=&-P_{j,j}^{-1}L_{j,j^{\prime}}.\end{split}$

For the second line, we used Equation 3. Equivalently, the Jacobian matrix

${\Bigl{(}\pdv{\delta_{j}}{\lambda_{j^{\prime}}}\Bigr{)}}_{j,j^{\prime}}=-P^{-1}L,$

for matrices $P$ and $L$ defined in the statement of the theorem. The Jacobian
$J_{\mathrm{QIS}}$ of the QIS iteration can be written as

$J_{\mathrm{QIS}}=\mathds{1}+{\Bigl{(}\pdv{\delta_{j}}{\lambda_{j^{\prime}}}\Bigr{)}}_{j,j^{\prime}}=\mathds{1}-P^{-1}L.$

This completes the proof the theorem. ∎

###### Proof of Theorem 4.2.

In an iteration of the algorithm, we have

$\delta_{j}=\eta\left\langle F_{j},\rho_{0}\right\rangle-\eta\left\langle
F_{j},Y^{(t)}/\operatorname{tr}Y^{(t)}\right\rangle.$

The Jacobian $J_{\mathrm{GD}}$ of each iteration has $(j,j^{\prime})$ entry

$\begin{split}&\mathds{1}-\eta\,\pdv*{\left\langle
F_{j},\xi\right\rangle}{\lambda_{j^{\prime}}}\\\
=\,&\mathds{1}-\eta\,\pdv*{\pdv*{\ln\operatorname{tr}\exp(\lambda\cdot
F)}{\lambda_{j}}}{\lambda_{j^{\prime}}}\\\ =\,&\mathds{1}-\eta L.\end{split}$

This completes the proof. ∎

We recall a theorem of Ostrowski which we will use to prove the convergence
rate by bounding the spectral radius of a Jacobian matrix.

######  Theorem B.1 (Ostrowski’s theorem (Ostrowski, 1966, Chapter 22)).

Assume function $f$ is differentiable at the neighborhood of a fixed point
$\zeta$. For an iterative algorithm $\zeta_{t+1}=f(\zeta_{t})$. A sufficient
condition for $\zeta$ to be a point of attraction is the spectral radius
$r(J_{f})<1$. Moreover, if $\zeta$ is an attraction point, the geometric
convergence rate of the iterative algorithm is given by

$\limsup_{t\to\infty}\frac{\left\lVert\zeta_{t+1}-\zeta\right\rVert}{\left\lVert\zeta_{t}-\zeta\right\rVert}=r(J_{f}).$

##  Appendix C Bounds on the Hessian Matrix

In this section, we prove the upper bound on the Hessian of the log-partition
function.

The proof uses a modified quantum belief propagation. The idea of quantum
belief propagation was studied in Hastings (2007) and we give a version of it
in the following lemma. It specifies how the matrix exponential function
changes with perturbations of the matrix.

######  Lemma C.1 (Quantum Belief Propagation (Hastings, 2007)).

Suppose $f_{\beta}(t)$ is the function whose Fourier transform is

$\tilde{f}_{\beta}(\omega)=\frac{\tanh(\beta\omega/2)}{\beta\omega/2}$

and $H(s)=H+sV$ for $s\in\interval{0}{1}$. Define the quantum belief
propagation operator

$\Phi_{H(s)}(V)=\int_{-\infty}^{\infty}\odif{t}\,f_{\beta}(t)\,e^{-iH(s)t}Ve^{iH(s)t}.$

Then

$\odv*{\exp(\beta H(s))}{s}=\frac{\beta}{2}\Bigl{\\{}\exp(\beta
H(s)),\Phi_{H(s)}(V)\Bigr{\\}}.$

Here, we introduce a modified version of it to prove that if the perturbation
is positive semidefinite, then so is the derivative of the matrix exponential
function. That is, the modified quantum belief propagation expresses the
derivative $\odv*{\exp(\beta H(s))}{s}$ so that its positivity is obvious for
positive $V$.

The proof uses the Bochner’s theorem and we give a simple version of it which
suffices for our purpose.

######  Lemma C.2 (Bochner’s Theorem).

A continuous function $f(x)$ on the real line with $f(0)=1$ is positive-
definite if and only if its Fourier transform is a probability measure on .

######  Lemma C.3.

Suppose $g_{\beta}$ is a function whose Fourier transform is

$\tilde{g}_{\beta}(\omega)=\frac{e^{\beta\omega/2}-e^{-\beta\omega/2}}{\beta\omega}$

and $H(s)=H+sV$ for $s\in\interval{0}{1}$. Define the modified quantum belief
propagation operator

$\Psi_{H(s)}(V)=\int_{-\infty}^{\infty}\odif{t}\,g_{\beta}(t)\,e^{-iH(s)t}Ve^{iH(s)t}.$

Then

$\odv*{\exp(\beta H(s))}{s}=\beta\,\exp\Bigl{(}\frac{\beta
H(s)}{2}\Bigr{)}\,\Psi_{H(s)}(V)\,\exp\Bigl{(}\frac{\beta H(s)}{2}\Bigr{)}.$
(5)

Furthermore, $g_{\beta}(t)$ is a probability density function over the real
line and $\Psi_{H(s)}$ is a completely positive trace-preserving map.

###### Proof of Lemma C.3.

Consider the spectrum decomposition of $H(s)$ as
$H(s)=\sum_{j}\lambda_{j}\lvert\mspace{0.5mu}\psi_{j}\mspace{0.5mu}\rangle\langle\mspace{0.5mu}\psi_{j}\mspace{0.5mu}\rvert$.

Using Duhamel’s formula, we have

$\begin{split}\odv*{\exp(\beta H(s))}{s}=\,&\int_{0}^{1}\odif{t}\,e^{t\beta
H(s)}\Bigl{(}\odv*{\beta H(s)}{s}\Bigr{)}e^{(1-t)\beta H(s)}\\\
=\,&\beta\int_{0}^{1}\odif{t}\,e^{t\beta H(s)}Ve^{(1-t)\beta
H(s)}.\end{split}$

Hence, the $(j,j^{\prime})$-th entry of $\odv*{\exp(\beta H(s))}{s}$ in the
basis of $\\{\lvert\mspace{0.5mu}\psi_{j}\mspace{0.5mu}\rangle\\}$ is

$\begin{split}\bigl{\langle}\mspace{0.5mu}\psi_{j}\mspace{0.5mu}\bigr{\rvert}\,\odv*{\exp(\beta
H(s))}{s}\,\bigl{\lvert}\mspace{0.5mu}\psi_{j^{\prime}}\mspace{0.5mu}\bigr{\rangle}=\,&\beta
V_{j,j^{\prime}}\int_{0}^{1}\odif{t}\,e^{t\beta\lambda_{j}+(1-t)\beta\lambda_{j^{\prime}}}\\\
=\,&\begin{cases}\;\dfrac{e^{\beta\lambda_{j}}-e^{\beta\lambda_{j^{\prime}}}}{\lambda_{j}-\lambda_{j^{\prime}}}V_{j,j^{\prime}}&\text{
if }\lambda_{j}\neq\lambda_{j^{\prime}}\\\ \;\beta
e^{\beta\lambda_{j}}V_{j,j^{\prime}}&\text{ o.w.}\end{cases}\end{split}$ (6)

where $V_{j,j^{\prime}}=\langle\mspace{0.5mu}\psi_{j}\mspace{0.5mu}\rvert
V\lvert\mspace{0.5mu}\psi_{j^{\prime}}\mspace{0.5mu}\rangle$.

Now we simplify the right-hand side of Equation 5. By the definition of
$\Psi_{H(s)}(V)$, the $(j,j^{\prime})$-th entry of $\Psi_{H(s)}(V)$ in the
basis $\\{\lvert\mspace{0.5mu}\psi_{j}\mspace{0.5mu}\rangle\\}$ is

$\int_{-\infty}^{\infty}\odif{t}\,g_{\beta}(t)\,e^{-i\lambda_{j}t}V_{j,j^{\prime}}e^{i\lambda_{j^{\prime}}t}=\tilde{g}_{\beta}(\lambda_{j}-\lambda_{j^{\prime}})\,V_{j,j^{\prime}}.$

Hence, the $(j,j^{\prime})$-th matrix entry of right-hand side in the basis
$\\{\lvert\mspace{0.5mu}\psi_{j}\mspace{0.5mu}\rangle\\}$ can be written as

$\beta\,e^{\beta(\lambda_{j}+\lambda_{j^{\prime}})/2}\,\tilde{g}_{\beta}(\lambda_{j}-\lambda_{j^{\prime}})\,V_{j,j^{\prime}}=\begin{cases}\;\dfrac{e^{\beta\lambda_{j}}-e^{\beta\lambda_{j^{\prime}}}}{\lambda_{j}-\lambda_{j^{\prime}}}V_{j,j^{\prime}}&\text{
if }\lambda\neq\lambda_{j^{\prime}}\\\ \;\beta
e^{\beta\lambda_{j}}V_{j,j^{\prime}}&\text{ o.w.}\end{cases}$

which is the same as the $(j,j^{\prime})$-th entry of the left-hand side by
Equation 6. This completes the proof of Equation 5. The fact that
$g_{\beta}(t)$ is a probability density function and the CPTP property of
$\Psi_{H(s)}$ follow from the Bochner’s theorem applied to $g_{\beta}$ and
$\tilde{g}_{\beta}$ and the fact that $\tilde{g}_{\beta}(0)=1$. ∎

To derive the bound in Theorem 4.3, we need two related results stated in
Theorems C.4 and C.5 which we now prove.

Define matrices

$\displaystyle\Delta$ $\displaystyle=\sum_{j}\pdv{\ell(\lambda\cdot
F)}{\lambda_{j}}\lvert\mspace{0.5mu}j\mspace{0.5mu}\rangle\langle\mspace{0.5mu}j\mspace{0.5mu}\rvert,$
$\displaystyle\Lambda$
$\displaystyle=\sum_{j,j^{\prime}}\pdv{\ell(\lambda\cdot
F)}{\lambda_{j},\lambda_{j^{\prime}}}\lvert\mspace{0.5mu}j\mspace{0.5mu}\rangle\langle\mspace{0.5mu}j^{\prime}\mspace{0.5mu}\rvert.$

$\Delta$ is a diagonal matrix and $\Lambda$ is the Hessian of the partition
function $Z=\ell(\lambda\cdot F)=\operatorname{tr}\exp(\lambda\cdot F)$.

######  Theorem C.4.

For $\Lambda$ and $\Delta$ defined above, we have $\Lambda\preceq\Delta$.

###### Proof of Theorem C.4.

Choose $H=\ln Y_{0}+\lambda\cdot F$, $\beta=1$, $s=0$, and $V=F_{j^{\prime}}$
in Lemma C.3, we have

$\begin{split}\Lambda_{j,j^{\prime}}&=\left\langle
F_{j},\pdv*{\exp(H)}{\lambda_{j^{\prime}}}\right\rangle\\\ &=\left\langle
F_{j},e^{H/2}\,\Psi_{H}(F_{j^{\prime}})\,e^{H/2}\right\rangle\geq
0.\end{split}$

That is, all entries of matrix $\Lambda$ (in the basis
$\bigl{(}\lvert\mspace{0.5mu}\psi_{j}\mspace{0.5mu}\rangle\bigr{)}$) are non-
negative.

Next, we prove that $\Delta-\Lambda$ is a diagonally dominant matrix. For all
$j^{\prime}$, the sum of the $j^{\prime}$-th column is

$\begin{split}\Delta_{j^{\prime},\,j^{\prime}}-\sum_{j}\Lambda_{j,\,j^{\prime}}=\,&\left\langle
F_{j^{\prime}},e^{H}\right\rangle-\left\langle\sum_{j}F_{j},e^{H/2}\Psi_{H}(F_{j^{\prime}})e^{H/2}\right\rangle\\\
\geq\,&\left\langle
F_{j^{\prime}},e^{H}\right\rangle-\left\langle\Psi_{H}(F_{j^{\prime}}),e^{H}\right\rangle\\\
=\,&\left\langle
F_{j^{\prime}},e^{H}\right\rangle-\left\langle\int_{-\infty}^{\infty}\odif{t}\,g_{1}(t)e^{-iHt}F_{j^{\prime}}e^{iHt},e^{H}\right\rangle\\\
=\,&\left\langle
F_{j^{\prime}},e^{H}\right\rangle-\int_{-\infty}^{\infty}\odif{t}g_{1}(t)\,\operatorname{Tr}\bigl{(}e^{-iHt}F_{j^{\prime}}e^{iHt}e^{H}\bigr{)}\\\
=\,&\left\langle
F_{j^{\prime}},e^{H}\right\rangle-\int_{-\infty}^{\infty}\odif{t}g_{1}(t)\,\operatorname{Tr}\bigl{(}F_{j^{\prime}}e^{H}\bigr{)}\\\
=\,&\left\langle F_{j^{\prime}},e^{H}\right\rangle-\left\langle
F_{j^{\prime}},e^{H}\right\rangle=0.\end{split}$

In the above, the inequality follows from $\sum_{j}F_{j}\preceq\mathds{1}$,
the positivity of $e^{H/2}\Psi_{H}(F_{j})e^{H/2}$ and the cyclic property of
the trace. The fifth line uses the commutativity of $e^{iHt}$ and $e^{H}$ and
the cyclic property of the trace. The last line follows from the fact that
$g_{1}(t)$ is the probability density function by the Bochner’s theorem.

The claim in the theorem now follows by the well-known matrix theory result
that diagonally dominant matrices are positive semidefinite. ∎

######  Lemma C.5.

For matrices $\Delta,\Lambda$ defined above and
$Q=\sum_{j,j^{\prime}}\left\langle F_{j}\right\rangle\left\langle
F_{j^{\prime}}\right\rangle\lvert\mspace{0.5mu}j\mspace{0.5mu}\rangle\langle\mspace{0.5mu}j^{\prime}\mspace{0.5mu}\rvert$,
we have the following identity

$\Lambda=Z(L+Q).$

###### Proof.

By definition, we have

$\begin{split}L_{j,j^{\prime}}&=\pdv*{\ln
Z}{\lambda_{j},\lambda_{j^{\prime}}}\\\
&=\frac{1}{Z^{2}}\Bigl{(}Z\,\pdv{Z}{\lambda_{j},\lambda_{j^{\prime}}}-\pdv{Z}{\lambda_{j}}\pdv{Z}{\lambda_{j^{\prime}}}\Bigr{)}\\\
&=\frac{1}{Z}\Lambda_{j,j^{\prime}}-\left\langle
F_{j}\right\rangle\left\langle F_{j^{\prime}}\right\rangle,\end{split}$

or equivalently $\Lambda=Z(L+Q)$ in the matrix form. ∎

We are now ready to prove the main result stated in Theorem 4.3.

###### Proof of Theorem 4.3.

By Lemma C.5, we have

$L=\frac{\Lambda}{Z}-Q.$

Hence, Theorem C.4 implies that

$L=\frac{\Lambda}{Z}-Q\preceq\frac{\Delta}{Z}-Q=P-Q.$

As $Q$ is outer product of vector $\sum_{j}\left\langle
F_{j}\right\rangle\lvert\mspace{0.5mu}j\mspace{0.5mu}\rangle$, it is a
rank-$1$ and positive semidefinite matrix. Therefore, we have

$L\preceq P-Q\preceq P$

which completes the proof. ∎

Finally, we will also need a lower bound on $L$, for which we recall a result
about the strong convexity of the log-partition function from Anshu et al.
(2021).

######  Theorem C.6 (Theorem 6 of Anshu et al. (2021)).

Let $H(\mu)=\sum_{j=1}^{m}\mu_{j}H_{j}$ be an $\ell$-local Hamiltonian over a
finite dimensional lattice. For a given inverse temperature $\beta$, there are
constants $c,c^{\prime}>3$ depending on the geometric property of the lattice
such that

$\nabla^{2}_{\mu}\ln\operatorname{tr}\bigl{(}e^{-\beta
H(\mu)}\bigr{)}\succeq\frac{e^{-O(\beta^{c})}\beta^{c^{\prime}}}{m}\mathds{1}.$

##  Appendix D Discussions on Quasi-Newton Methods

In this section, we give some details of the Anderson mixing method and L-BFGS
method.

The Anderson mixing method interpolates history information in order to speed
up a fixed-point iteration. More concretely, suppose $g:^{d}\rightarrow^{d}$
is a contraction and we are interested in finding the fix-point $x=g(x)$. The
standard fix-point iterative algorithm is to compute
$x_{t+1}=g\bigl{(}x_{t}\bigr{)}\text{ for }t=0,1,2,\ldots$, until a stopping
criteria is met. In Anderson mixing, a relatively small history size $m\geq 0$
is chosen and we define $m_{t}=\min\\{m,t\\}$. In our numerical
implementation, we use $m=10$. Define the residual
$r_{t}=g\bigl{(}x_{t}\bigr{)}-x_{t}$ and two matrices $X_{t},R_{t}\in^{d\times
m}$ storing the historical information

$\displaystyle X_{t}$ $\displaystyle=\bigl{(}\Delta x_{t-m_{t}},\Delta
x_{t-m_{t}+1},\ldots,\Delta x_{t-1}\bigr{)},$ $\displaystyle R_{t}$
$\displaystyle=\bigl{(}\Delta r_{t-m_{t}},\Delta r_{t-m_{t}+1},\ldots,\Delta
r_{t-1}\bigr{)},$

where $\Delta$ is the forward difference operator and $\Delta
x_{k}=x_{k+1}-x_{k}$. Then, the Anderson accelerated iteration can be written
succinctly as $x_{t+1}=x_{t}+G_{t}r_{t}$ where

$G_{t}=\beta_{t}I-(X_{t}+\beta_{t}R_{t}){\bigl{(}R_{t}^{T}R_{t}\bigr{)}}^{-1}R_{t}^{T}.$

Here, $\beta_{t}$ is the mixing parameter.

It is pointed out in Fang & Saad (2009) that $G_{t}$ approximates the inverse
of the Jacobian of $g$ and Anderson mixing method can be thought of as a
quasi-Newton method satisfying multi-secant equations. We note that there is a
matrix inverse in the above formula for $G_{t}$ which can be implemented using
Moore-Penrose pseudo-inverse. For stability and efficiency concerns, we found
that the AM algorithms have the best performance in our numerical simulations
when using a relative condition number of $1\mathrm{e}{-}7$ in the pseudo-
inverse, a cutoff threshold that sets small singular values of the matrix to
zero. This is easily implemented by setting the $\mathrm{rcond}$ parameter of
the $\mathrm{pinv}$ function in $\mathrm{numpy.linalg}$ package for Python
implementations.

The BFGS method is one of the most popular quasi-Newton methods that can be
applied to unconstrained optimization problems $\min_{x\in^{n}}f(x)$. It is
also known as the variable metric algorithm as first proposed by Davidon
(Davidon, 1991; Yuan, 2015). The algorithm maintains the approximate Hessian
$H_{k+1}$ of the optimization problem. The update rule of the algorithm is

$\displaystyle x_{t+1}=x_{t}-\eta_{t}H_{t}\nabla f(x_{t}),$

and the update rule of $H_{t}$ is

$\displaystyle
H_{t+1}=\Bigl{(}\mathds{1}-\frac{s_{t}y_{t}^{T}}{y_{t}^{T}s_{t}}\Bigr{)}H_{t}\Bigl{(}\mathds{1}-\frac{y_{t}s_{t}^{T}}{y_{t}^{T}s_{t}}\Bigr{)}+\frac{s_{t}s_{t}^{T}}{y_{t}^{T}s_{t}},$

where $s_{t}=x_{t+1}-x_{t}$, $y_{t}=\nabla f(x_{t+1})-\nabla f(x_{t})$ and
$H_{0}$ is a predefined initial approximation of the inverse Hessian matrix.
We refer readers to Nocedal & Wright (2006) for a discussion on how the BFGS
update is derived. To optimize the memory usage, the limited memory version of
BFGS called L-BFGS (Liu & Nocedal, 1989) is used. Since using line search for
choosing $\eta_{t}$ can incur additional function evaluations in each
iteration, we use $\eta_{t}\equiv 1$ and only tune $H_{0}$ in the numerical
simulation.

Both the AM and the BFGS methods employed our numerical simulation can be
further strengthened by using a heuristics invented by Barzilai and Borwein
(Barzilai & Borwein, 1988) to choose the $\beta_{t}$ in AM and $H_{0}$ in
BFGS. In AM, we set

$\displaystyle\beta_{t}=-\frac{\Delta r_{t-1}^{T}\Delta x_{t-1}}{\Delta
r_{t-1}^{T}\Delta r_{t-1}},$

which solves $\min_{\beta}\left\lVert\Delta x_{t-1}+\beta\Delta
r_{t-1}\right\rVert_{2}$. In BFGS, we set

$\displaystyle H_{0}=\frac{y_{t-1}^{T}s_{t-1}}{y_{t-1}^{T}y_{t-1}}\mathds{1}$

to be the initialization of the approximate inverse Hessian in the $t$-th
iteration (Nocedal & Wright, 2006, Pages 143 and 178).

##  Appendix E Experimental Details

In this section, we provide details of the numerical experiments.

###  E.1 Hamiltonian Generation

Three types of Hamiltonians are used in our numerical experiments. These
include the Ising family of Hamiltonians, the transversal $2$-local
Hamiltonians in 1D (Transversal1D), and the $2$-local Hamiltonians in 1D
(Local1D). The figures in the main text demonstrate the results of the
experiments on the Local1D Hamiltonian in $6$-qubit system.

#### Ising Hamiltonian

The Ising Hamiltonian is a sum of $2$-local terms, which can be written as

$H=\sum_{i=1}^{n}a_{i}\sigma_{x}^{(i)}+\sum_{i=1}^{n-1}b_{i}\sigma_{z}^{(i)}\sigma_{z}^{(i+1)},$

where $\sigma_{z}^{(i)}$ is a Pauli $Z$ operator acting on the $i$-th qubit.
The coefficients $a_{i}$ and $b_{i}$ are randomly generated from the standard
normal distribution.

#### Transversal1D Hamiltonian

The transversal 1D Hamiltonian is a sum of $2$-local terms, and it can be
written as

$H=\sum_{p\in\\{x,y,z\\}}a_{p}\sum_{i=1}^{n}\sigma_{p}^{(i)}+\sum_{p,q\in\\{x,y,z\\}}b_{p,q}\sum_{i=1}^{n}\sigma_{p}^{(i)}\sigma_{q}^{(i+1)}.$

The coefficients $a_{j}$ and $b_{p,q}$ are randomly generated from the
standard normal distribution.

Notice that, in this scenario, there are $12$ terms in the Hamiltonian, and
each term in the Hamiltonian is a sum of Pauli operators acting on different
qubits. Thus, the operator norm of each term is bounded by the number of
qubits instead of $1$. As a result, in the implementation of our algorithm, we
need to normalize the Hamiltonian by dividing the number of qubits.

#### Local1D Hamiltonian

The local 1D Hamiltonian is a sum of $2$-local terms, including the Ising
Hamiltonian as a special case, and it can be written as

$H=\sum_{i=1}^{n}\sum_{p\in\\{x,y,z\\}}a_{p,i}\sigma_{p}^{(i)}+\sum_{i=1}^{n}\sum_{p,q\in\\{x,y,z\\}}b_{p,q,i}\sigma_{p}^{(i)}\sigma_{q}^{(i+1)},$

where $\sigma^{(i)}$ is a Pauli operator acting on the $i$-th qubit for $i\leq
n$, and $\sigma^{(n+1)}$ is defined as $\sigma^{(n+1)}=\sigma^{(1)}$. The
coefficients $a_{p,i}$ and $b_{p,q,i}$ are randomly generated from the
standard normal distribution.

###  E.2 Initialization, Learning Rate, and Stopping Criteria

The initial guess for the parameters in the QIS algorithm and the gradient
descent algorithm is set to $0$. The learning rate in the gradient descent
algorithm is set to the number of terms in the Hamiltonian. In the
experiments, we first generate the Hamiltonian and then normalize it by
dividing it by the number of qubits. Then, the inverse temperature $\beta$ is
set to $1$. Therefore, the optimal solutions and optimal objective values can
be directly computed and compared with the results of the algorithms. With
these settings, the stopping criteria can be set to the error of the objective
value. The stopping criteria for the QIS, gradient descent algorithm, and
quasi-Newton accelerations is set to have error at most $10^{-12}$. The
maximum number of iterations for the quasi-Newton accelerations is set to
$40$.

###  E.3 Results on Different Hamiltonians

In this section, we provide the results of the numerical experiments on
different Hamiltonians. We have conducted experiments on the Ising,
Transversal1D, and Local1D Hamiltonians, with different numbers of qubits. The
results are shown in the following tables and figures. We first present the
results in Table 1 and then provide the figures for each family of
Hamiltonians.

| QIS | GD | AM-QIS | L-BFGS-GD
---|---|---|---|---
6-qubit-Ising | 319 | 479 | 7 | 6
6-qubit-Transversal | 1743 | 2562 | 5 | 7
6-qubit-Local | 1913 | 2377 | 5 | 8
7-qubit-Ising | 370 | 551 | 6 | 6
7-qubit-Transversal | 2006 | 2940 | 5 | 7
7-qubit-Local | 2236 | 2721 | 5 | 7
8-qubit-Ising | 421 | 623 | 6 | 5
8-qubit-Transversal | 2267 | 3311 | 4 | 7
8-qubit-Local | 2525 | 3020 | 5 | 7
Table 1: Number of steps to achieve $10^{-7}$ precision for different
Hamiltonians.

The results of the Ising Hamiltonian are shown in Figure 5 for a $7$-qubit
system with $13$ Hamiltonian terms. In such cases, the QIS algorithm converges
to the optimal solution with error $10^{-9}$ in roughly $400$ iterations with
less loss than the GD algorithm, while the quasi-Newton methods converge to
the optimal solution with similar error in less than $20$ iterations.

(a) QIS and GD results. (b) Quasi-Newton methods without BB method. (c) Quasi-
Newton methods with BB method.

Figure 5: Results of the Ising Hamiltonian for $7$ qubits.

The results of the Transversal1D Hamiltonian are shown in Figure 6 for a
$7$-qubit system with $12$ terms. In such cases, the QIS algorithm converges
to the optimal solution with error $10^{-7}$ in roughly $1000$ iterations with
less loss than the GD algorithm, while the quasi-Newton methods converge to
the optimal solution with similar error in less than $20$ iterations.

(a) QIS and GD results. (b) Quasi-Newton methods without BB method. (c) Quasi-
Newton methods with BB method.

Figure 6: Results of the Transversal1D Hamiltonian for $7$ qubits.

The results of the Local1D Hamiltonian are shown in Figure 7 for a $7$-qubit
system with $96$ terms. In such cases, the QIS algorithm converges to the
optimal solution with error $10^{-9}$ in roughly $2000$ iterations with less
loss than the GD algorithm, while the quasi-Newton methods converge to the
optimal solution with similar error in less than $40$ iterations.

(a) QIS and GD results. (b) Quasi-Newton methods without BB method. (c) Quasi-
Newton methods with BB method.

Figure 7: Results of the Local1D Hamiltonian for $7$ qubits.

###  E.4 Reproducibility Details

The numerical experiments were conducted using Python 3.11.5 in Python
scripts. Running the script requires the following packages: numpy and scipy
for linear algebra operations, matplotlib for data visualization, sys for
outputting results in txt format, time for measuring the runtime of the
algorithm, and Qiskit (Qiskit, 2023) for generating Hamiltonians. For
reproducibility, we set the random seed to $100$ for QIS and GD, and random
seed array from $1$ to the number of experiments for quasi-Newton
accelerations. The average runtime of the script is less than 6 minutes for
systems of five to eight qubits (of dimensions $32\times 32$ to $256\times
256$).