Quantum Natural Gradient with
Efficient Backtracking Line Search

Touheed Anwar Atif,  Uchenna Chukwu

Jesse Berwald and Raouf Dridi

Quantum Computing Inc.

We consider the Quantum Natural Gradient Descent (QNGD) scheme which was recently proposed to train variational quantum algorithms. QNGD is Steepest Gradient Descent (SGD) operating on the complex projective space equipped with the Fubini-Study metric.  Here we present an adaptive implementation of  QNGD based on Armijo's rule, which is an efficient backtracking line search that enjoys a proven convergence.
The proposed algorithm is tested using noisy simulators on three different models with various initializations. Our results show that Adaptive QNGD dynamically adapts the step size and consistently outperforms the original QNGD, which requires knowledge of optimal step size to perform competitively. In addition, we show that the additional complexity involved in performing the line search in Adaptive QNGD is minimal, ensuring the gains provided by the proposed adaptive strategy dominates any increase in complexity. Additionally, our benchmarking demonstrates that a simple SGD algorithm (implemented in the Euclidean space) equipped with the adaptive scheme above, can yield performances similar to the QNGD scheme with optimal step size.

Our results are yet another confirmation of the importance of differential geometry in variational quantum computations.  As a matter of fact, we foresee advanced mathematics to play a prominent role in the NISQ era in guiding the design of faster and more efficient algorithms.

§ INTRODUCTION

Quantum processors currently suffer from low qubit counts and short coherence times, which restricts the number (and the depth) of quantum algorithms that can be executed. An approach to deal with this limitation is the  variational approach [1, 2, 3], where the short-lived qubits are recycled multiple times. A classical optimizer  (e.g., SGD) is entrusted with administering part of the dynamics. The quantum component runs a parametric circuit, an ansatz, given the parameter values provided by the classical controller.  Conversely, the  classical optimizer, updates the circuit parameters by minimizing the expectation value of the parametric circuit. The term hybrid quantum-classical neural network is often used and puts variational algorithms in the wider context of quantum machine learning [4, 5].

A recent and attractive proposal to train variational algorithms was outlined in [6], which employs a SGD, preconditioned with the inverse of the  Fubini-Study tensor. For a geometer,  the parameters update takes place on the Kähler manifold of the pure quantum states (i.e., the complex projective space).  The idea was used before in [7] in classical machine learning. Subsequent results ([6, 8, 9]) showed that this preconditioning does indeed lead to a better performance and faster training.

Here we take this geometrical approach further and provide
an adaptive implementation of the QNGD scheme, thus mitigating the dependence of its performance  on the step size. For this, we have employed the so-called Armijo rule, which is an efficient backtracking line search scheme with a proved convergence [10].

We have tested our scheme on three different models with various initializations. Our results show that the Adaptive QNGD scheme dynamically adapts the step size and consistently outperforms QNGD which requires the knowledge of the optimal step size to perform competitively.  Remarkably, our results also show that
a simple steepest gradient descent  equipped with the adaptive scheme above, can yield performances similar to the QNGD with optimal step size.

This work is organized as follows. In Section <ref>, we briefly review the concepts behind the QNGD technique, and  highlight
its implementation on quantum processors. In  Section <ref>, we introduce our scheme  and  compare its complexity with that of the QNGD scheme. We show that the additional circuit evaluations needed  are insignificant.

Section <ref> contains the results of our benchmarking conducted on three different problems. We conclude with a brief summary.

§ ACRONYMS

AdaptQNGD       Adaptive Quantum Natural Gradient Descent

qngd       Quantum Natural Gradient Descent

sgd       Steepest
Gradient Descent

realAmplitude     ansatz alternating $R_Y(\theta)$ and cnot layers

FullFubiniMetric     full Fubini-Study metric calculation procedure

§ QUANTUM NATURAL GRADIENT DESCENT

The  complex projective space
$\mathbb C\mathbb P^N$ is identified with the quotient space $S^{2N+1}/U(1)$ where $U(1)$ is the circle group.

Given two (normalized) pure quantum states $|\psi\rangle$ and $|\varphi\rangle$ in $\mathbb C\mathbb P^N$,  the Fubini-Study distance between  the two is given by

\begin{equation}
d_{FS} (|\psi\rangle, |\varphi\rangle)
= \mathrm{arccos} \sqrt{ {\langle \psi| \varphi\rangle \langle\varphi |\psi\rangle}{}}.
% {\langle \psi|\psi\rangle\langle \varphi|\varphi\rangle}},
\end{equation}

Let $N=2^n$. For convenience,  we place   ourselves in $\mathbb C\mathbb P^{N-1}$. Let us assume    $|\varphi\rangle = U(\boldTheta) |0\rangle,$  where the ansatz   $ \{U(\boldTheta), \, \boldTheta\in  \mathbb R^p\}$ is a smooth submanifold of the Lie group $U(2^n)$.
Passing to  the infinitesimals gives

the quantum geometric tensor  [6]

\begin{equation}\label{FS}
G_{ij} = \left\langle\partial_{\boldTheta_i}\varphi | \partial_{\boldTheta_j}\varphi \right\rangle
- \left\langle \partial_{\boldTheta_i}\varphi | \varphi \right\rangle
\left\langle \varphi| \partial_{\boldTheta_j}\varphi \right\rangle.
\end{equation}

The real part $F=\mathcal{R}e(G)$, called the Fubini-Study metric, is positive semi-definite, Hermitian, and is the unique metric compatible with the quotient structure above. Formally, $\mathbb C\mathbb P^N$ equipped with the Fubini-Study metric is a Kähler manifold.  This rich geometrical structure implies that the gradient descent update has the form

\begin{equation}
\boldTheta  \leftarrow \boldTheta - \lambda \nabla f(\boldTheta),
\end{equation}

\begin{equation}
\nabla f(\boldTheta) =  F(\boldTheta)^{-1}\nabla_{\mathrm{Eucl}} f(\boldTheta),
\end{equation}

for any smooth real-valued function $f$ on $\mathbb C\mathbb P^{N-1}$.

All ansatzes considered in this paper are
structured with layers. An example of this is given by the following parametric circuit

$\ket{0}$     R_Y(2a)     1     R_Y(2c)

$\ket{0}$     R_Y(2b)          R_Y(2d)

which has two parametric layers: a first layer of parametric rotations, $R_Y(\theta) = e^{-i\frac{\theta}{2}Y}$, followed by a non-parametric layer consisting of a single cnot gate,  followed by a final parametric layer of $Y$ rotations. This type of ansatze (alternating is $R_Y(\theta)$ cnot gates) is referred to as a realAmplitude ansatz in the literature.

Given this layer-based structure, the Fubini-Study metric is structured into block-diagonal (a block for each layer)  and off-diagonal parts. The former are readily computed as in  [6]. For the off-diagonal parts, which are slightly less straightforward to implement, we have used
a method similar to the one proposed in [12] which adds an extra ancilla qubit to the original circuit. For instance, continuing with the ansatz above,  the  circuit for the first component of the off-diagonal term   $G_{ac} = \left\langle\partial_a\varphi | \partial_c\varphi \right\rangle$ is given by

$\ket{0}$         Y        R_Y(2a)     1    Y

$\ket{0}$             R_Y(2b)

$\ket{0}$    H    -2           -2   H

Here, the second component in equation <ref> vanishes, since the prepared quantum state has only real amplitudes.

We shall refer  to the full Fubini-Study metric calculation procedure as FullFubiniMetric.

§ ADAPTIVE QUANTUM NATURAL GRADIENT DESCENT

§.§ Illustrative Example

Before we give the details of the proposed scheme, let us illustrate how the new adaptive scheme compares with traditional  QNGD scheme.

We will do so using the  problem of
finding the electronic ground state energy  (after Born-Oppenheimer approximation) of  the
Hydrogen molecule [13].  Its fermionic Hamiltonian is mapped (using Jordan–Wigner transformation) into the two-qubit Hamiltonian

\begin{align} \label{eq:h2hamiltonian}
H = \alpha_0 I_1 I_2 + \alpha_1 Z_1I_2 + \alpha_2 I_1Z_2 + \alpha_3 Z_1Z_2 + \alpha_4 Y_1Y_2 + \alpha_5 X_1X_2.
\end{align}

$I_i$,  $X_i, Y_i$ and $Z_i$ are resp., the identity operator and the Pauli operators acting on the $i$th qubits. The coefficients $\alpha_j$ are real-valued functions of the inter-atomic distance between the two Hydrogen nuclei.
For this Hamiltonian, we consider the simplified ansatz

$\ket{0}$     R_Y(θ_1)     1

which prepares the state $|{\psi(\theta_1,\theta_2)}\rangle$. The Fubini-Study metric is

\begin{align*}
F(\theta_1,\theta_2) = \frac{1}{4}\left[\begin{array}{cc}
1 & 0
\\
0 & 1
\end{array}\right],
\end{align*}

which is, in this simple case, independent of both $\theta_1$ and $\theta_2$.

Figure <ref> shows the training trajectories taken by QNGD with different learning rates $\eta$ and by the proposed adaptive scheme, all starting from the same initial point   $\theta_1 =-0.1$ and $\theta_2 = -0.2$.   The results show that,
not only does AdaptQNGD  converge much faster than the three other schemes, but it also uses a large step-size in the early stages of the algorithm (when the surface allows for this behaviour) and decreases the steps size closer to the minima, effectively slowing down which avoids the arbitrary bouncing as observed with step size $1$.

Upper Left: The contour plot of the Hydrogen Hamiltonian with demo ansatz.  Upper Right: Comparison of training paths taken by different schemes to approach the minima, with maximum number of epochs set to 20. The QNGD scheme with step size $\eta =1 $ does not converge. The QNGD with step sizes $=0.5$ and $0.25$ and the AdaptQNGD schemes all converge, however, AdaptQNGD  converges much faster. Lower left: The figure shows the convergence of average value of the Hamiltonian with iteration counts for each of the schemes discussed in the center plot. Lower right: The figure depicts the step sizes chosen by the AdaptQNGD over the course of the algorithm. These figures demonstrate that the algorithm is able to exploit larger step sizes when the surface allows, and tune it lower upon reaching close to the minima, resulting in faster overall convergence.

§.§ Main Algorithm

AdaptQNGD  algorithm

Input parameters: $(\alpha,\beta,k_m,tol,\epsilon)$ – See table below for description.

Step 1 Initialize: $\boldTheta_0, i=0.$

Step 2 Compute the Euclidean gradient:

$\nabla_{\mathrm{Eucl}} f(\boldTheta_i) = $ ParameterShift$(\boldTheta_i)$.

Step 3  Compute the Fubini-Study gradient:

$F(\boldTheta_i)$ =  FullFubiniMetric$(\boldTheta_i)$,

$F(\boldTheta_i)^{-1}$ = pseudoInvert  $\left(F(\boldTheta_i), \epsilon\right)$,

$ \nabla f(\boldTheta_i) :=  F(\boldTheta_i)^{-1}\nabla_{\mathrm{Eucl}} f(\boldTheta_i)$.

Step 4 Terminate if the stopping criteria (defined below) is met, else

(1)  Compute $k_i$:

\begin{align}
k_i &\deq \min\left\{k\in[0,k_{m}]\big|
f(\boldTheta_i) - f(\boldTheta_i - \frac{\beta}{2^k}\nabla f(\boldTheta_i)) \geq \alpha\frac{\beta}{2^{k}}\|\nabla f(\boldTheta_i)\|_2^2) \right\}, \label{eq:armijos}
\end{align}

where $\|\|_2$ is the Euclidean norm.

(2) Update step size to $\lambda_i = \beta/2^{k_i}$,

(3) Update parameters as $\boldTheta_{i+1} = \boldTheta_i - \lambda_i \nabla f(\boldTheta_i),$ and proceed to Step 2.

The five input parameters are defined as follows:

* $\alpha \in (0,1) $ is a constant to tune the sensitivity of the Armijo's line search principle. Lower values of $\alpha$ skews the rule toward preferring higher step sizes and higher skews it toward lower step sizes.

* $\beta > 0$ is the maximum rate used by the algorithm.

* $k_{m} > 0$ is the maximum number of steps searched by the line search algorithm in finding the best step size.

* $tol > 0$  governs the stopping criteria defined below.

* $\epsilon > 0$ a tolerance parameter used in the pseudoInvert function described below. For all our simulations, $\epsilon$ is set to $10^{-3}$.

The ParameterShift procedure in Step 2 is well known and defined in different places, including [14].  The pseudoInvert function is a variant of the Moore–Penrose pseudo-inverse  where apart from the zero eigenspace, we also ignore the eigenspace having eigenvalues less than $\epsilon$.
It doesn't  satisfy all Moore–Penrose properties[15] but works very well in the simulations considered here, avoiding the unreasonably large jumps caused at the singularities of the Fubini metric $F$. This can also be seen as a variant of Tikhonov regularization [16], which has shown to circumvent the effect of barren plateaus [9].

The stopping criteria is defined by the difference between the computed  and the exact energies i.e., $|f(\boldTheta)-f_{\mathrm{Exact}}(\boldTheta)|\leq tol$. In cases when computing the exact energies is infeasible, one can resort to other techniques developed in the gradient descent literature, for instance

$||\nabla f(\boldTheta_i)|| \leq tol$, or even $\|f(\boldTheta_{i-1}) - f(\boldTheta_i))\| \leq tol$.

The key point to note here is that step 4 shows that the line search is essentially a binary search. This means that our  scheme AdaptQNGD compares  reasonably with the “vanilla" QNGD which requires
$O(pM+p^2M) = O(p^2M)$ QPU calls–where $M$ is the total number of iterations to achieve the desired tolerance, and $p$ is  the total number of parameters.  For AdaptQNGD the total number of QPU calls is $O(p^2M'+M'\log_2(k_m))$ or  $O(p^2M')$ since
$\log_2(k_m) < p^2$, where $M'$ (which is less or equal to $M$)

is the number of iterations the proposed scheme takes to converge. This favourable asymptotic behaviour is key for the future applicability of our adaptive scheme.

§ EXPERIMENTATION

§.§ Setup

We benchamarked our algorithm on three models, the Hydrogen and Lithium Hydride molecules and the transverse field Ising model. We have  used the RealAmplitude ansatze introduced in Section <ref>. All procedures used were implemented with TensorFlow Quantum library.

§.§.§ The  Hydrogen Model
We took a closer look into the Hydrogen molecule discussed in the illustrative example in Subsection <ref>.  In particular, we wanted to compare against the results presented in [8]. Therefore,  we have used the same Hamiltonian (with $\alpha_0 = \alpha_3 = \alpha_4 = 0, \alpha_1  = \alpha_2 = 0.4, $ and $\alpha_5 = 0.2$)

\begin{equation}
H = 0.4 (ZI + IZ) + 0.2 XX,
\end{equation}

which has the spectrum $\{-\sqrt{17}/5, -1/5, 1/5, \sqrt{17}/5 \}$, with  $-\sqrt{17}/5 \approx -0.82462$ being the lowest energy (our target).  We have also considered the same ansatz

$\ket{0}$     R_Y(2a)     1     R_Y(2c)

$\ket{0}$     R_Y(2b)          R_Y(2d)

We randomly initialize its parameters by

choosing 100 different initial points, uniformly distributed over $[0,2\pi]$. For each of these initial points, we first simulate the regular QNGD scheme for different rates, and then simulate the AdaptQNGD scheme with $\alpha = 0.01$ and ${\beta = 0.5}$.

§.§.§ Lithium Hydride (LiH) Model

The Hamiltonian for LiH operates on four qubits and has 99 4-local Pauli terms. The complete list of these local Pauli  terms can be found in [17].
We employ again the RealAmplitude ansatz, now   with six parameterized layers, with a full entangling layer between every two of the former. The number of layers, compared to the H$_2$ model, is greater here to make the ansatz sufficiently expressive [17].   To initialize the parameters of the ansatz, a similar random strategy as above is used, with simulations for 100 different initial points. For each of these initial points, the QNGD scheme is simulated with different step sizes, and finally the AdaptQNGD   is simulated using $\alpha = 0.01$ and $\beta = 1$.

§.§.§ Transverse Field Ising (TFI) Model

The Hamiltonian is used in [9] and  is given by

\begin{equation}
H_{TFIM} = - \left(\sum_{1\leq i\leq N} Z_i Z_{i+1} + t \sum_{1\leq i\leq N} X_i
\right)
\end{equation}

with $t$, a real positive parameter, and the periodic condition, $N+1=1$, on the two sums (the $N$ here is a different notation than the one used in the Section <ref>).

We fix $t$ to $0.1$ and  vary $N$ in the interval $[2,10]$. We have used two setups:

* In the first setup, we choose the ansatz to have a minimum number of parameters while being expressive enough. Our simulations show that we need at least two parameterized layers with a full entanglement layer in-between.

* For the second setup, we choose the ansatz to study the behaviour of the proposed scheme in handling the over-parameterization problem. For this, we choose the number of parameterized layers equivalent to the size of the problem $(N)$, while enclosing an entanlging layer between each parameterized layers.  For instance, for $N=5$, we use five parameterized layers and four entangling layers.

§.§ Results

In this section, we report our results and findings obtained using AdaptQNGD   for the different models introduced above. Below we have set $\alpha = 0.01, \beta = 0.5$ and $k_{m} = 6$, independent of the model, the instance size and the ansatz depth.  The common thread below is that slight deviations in the optimal step size substantially degrades the performance of the various QNGD algorithms. At the same time, our AdaptQNGD scheme was able to automatically adjust its step size and match the optimal QNGD. Moreover, the AdaptQNGD was able to capture the correct dynamical change of the step size throughout the training and thus outperforming QNGD even with the optimal initial step size.

§.§.§ The  Hydrogen Model

Left: Variation of the median number of epochs needed to converge for 10 random initializations of the H$_2$ ansatz parameters with learning rate. The box plots indicates the 25th and the 75th percentile of the epochs required to converge. Right: Comparison of CDFs of the number of epochs needed for the ADAPT QNGD  and the QNGD, with varying step-sizes, to reach the lowest energy of the Hydrogen Hamiltonian.

sAs discussed earlier, we chose $100$ different initial points. For each of these initial points, we simulate the QNGDS for step sizes $0.125, 0.25, 0.5$ and $1$.

We set the termination tolerance $tol$ to  $0.01$. We refer to  the number of epochs, for a given scheme to terminate, as Epochs to Terminate (EoT).

Figure <ref> (Left) shows the variation of the median number of epochs needed to converge to the minimum eigenvalue of the Hydrogen Hamiltonian with the learning rate (step size) of the QNGD scheme. The median is taken across 10 random initial runs. It is clear from the plot that only for a range of learning rates, the QNGD scheme actually converges, and the rate of convergence highly depends on the chosen learning rate.

Figure <ref> (Right) shows the Cumulative Distribution Functions (CDF) of EoT for the Hydrogen model using the QNGDs with different step sizes and the AdaptQNGD. In addition, the dotted black line at ordinate $0.5$ is drawn to compare the performance of different schemes at their medians. Figure <ref>  clearly shows the dependence of the performance of the QNGD scheme on the choice of the step-size, with the step size of $0.5$ performing the best among others simulated. Doubling the step size to $1$ degrades the performance extensively causing almost all the runs to be requiring more than $100$ EoTs. On the other hand, decreasing the step sizes to $0.25$ and $0.125$  leads to sub-optimal performance. In general,
it  is  difficult to predict beforehand the optimal step size, which requires extensive numerical simulations.

Figure <ref> also shows  that the AdaptQNGD, with the parameters mentioned above, performs similar to the QNGDs with the optimal step size $0.5$. Note that the AdaptQNGD   only requires the step sizes it has to choose from, eliminating the need for the knowledge of the optimal step size.

§.§.§ Lithium Hydride Model

Comparison of CDFs of the number of epochs needed for the AdaptQNGD  and the QNGDS to reach the lowest energy of the LiH Hamiltonian, with (a) $tol = 0.01$ and (b) $tol = 0.001$.

The simulations were ran for $100$ different random initial points, as discussed earlier, and for each initial point, the QNGDS was simulated with step sizes $0.125, 0.25, 0.5$ and $1$.

Figures <ref> and <ref> shows the distribution of EoT for each of the schemes, using $tol = 0.01$ and $tol = 0.001$, respectively.

Both figures   confirm the dependence of the performance of  the  QNGD scheme on the choice of the step-size. Notice that in this model, the optimal step size is $0.5$.

In the same time,
the AdaptQNGD  was able to perform close to the best performance obtained by the QNGD scheme.
In fact, the performance in Figure <ref>, of the former is slightly better than the latter. This is due to the fact that, although the step size of $0.5$ seems to be a better choice among the others simulated here, it may be the case that it is not the best for the entire duration of convergence. The landscape of the cost function may at times allow for a faster descent (or convergence), but the QNGD scheme is always limited to using only a fixed step size for the entire duration of the algorithm. Here again, the AdaptQNGD learns this dynamically and proceeds with choosing larger step sizes than $0.5$ whenever possible, and as a result can converge faster.

§.§.§ Transverse Field Ising (TFI) Model

Comparison of median number of epochs needed of the AdaptQNGD  with the QNGD scheme   using Boxplots. The box around the median points indicates the containment of 25th to 75th percentiles .

In this last model, the simulations were ran for $25$ different random initial points, and for each initial point the QNGD scheme was simulated with step sizes $0.03125, 0.0625, 0.125,$ and $0.25$ and, as in the previous models, the AdaptQNGD scheme was simulated with $\alpha= 0.01, \beta = 1$ and $k_{m} = 6.$

For each size $N$, we first plot the median of the distribution of epochs to terminate, and then around the median, plot a box containing the realizations between the $25$th and the $75$th percentiles to obtain a comprehensive view of the performance.

Figure <ref> shows our first result for the TFI model obtained using ${tol = 0.01}$.

This set of simulations again show that QNGD scheme depends on the choice of step size, in which the optimal step size in this case equals
$0.125$. The performance degrades  considerably
with step size $0.25$ or higher. Expectedly, the AdaptQNGD   overcomes this and adaptively chooses the best step size.  Note that as the problem size increases, the number of epochs needed to terminate also increases, a phenomena also seen in [9].

Comparison of median number of epochs needed of the ADAPT SGD scheme with the QNGD scheme  using Boxplots. The box around the median points indicates the containment of 25th to 75th percentiles .

As a next result, we provide a comparison between the QNGD scheme and the Steepest Gradient Descent (SGD) equipped with the adaptive procedure used in AdaptQNGD (i.e., setting the tensor $F$ to identity in AdaptQNGD). This is an interesting comparison since the vanilla SGD  is blind to the underlying geometry (the complex projective space with the Fubini-Study metric), which results in
much slower convergence than the QNGD scheme [8]. Remarkably, by using the adaptive step size technique, which defaults to choosing a large step size when permissible by the landscape, we are able to obtain similar performance as the optimal QNGD scheme. Figure <ref> demonstrates this behaviour. This can prove to be very useful, as avoiding the computation of Fubini matrix leads to a substantial decrease in the number of circuit evaluations and henceforth the cost of simulating the algorithm can be significantly smaller.

§.§.§ Over-parameterization

Comparison of median number of epochs needed of the ADAPT QNGD scheme with the QNGDS for varying rates, using Boxplots. The box around the median points indicates the containment of 25th to 75th percentiles.

We conclude this section with another interesting result.We have simulated the TFI model in the context of  over-parametrization. Figure <ref> shows the results obtained. We observe here a concerning aspect of the QNGD scheme. Note that, the rate which performs superior to all other rates for lower $N$, eventually degrades in performance as $N$ is increased, when employing the QNGD scheme. For instance, the rate of $0.125$ performs the best, on average, relative to other rates, for $N$ less than $5$, but significantly deteriorates for any $N$ larger than $5$. As a result, a new rate dominates amongst the set of rates, in terms of performance, as $N$ changes. This makes it exhausting to compute the best rate for every problem size (i.e., $N$) individually. AdaptQNGD  was able to choose the most favourable rate for every problem size and performs equivalent to the best QNGD performance without being aware of the optimal rate.

§ CONCLUSION

In this paper, we have proposed an adaptive and efficient training scheme for variational quantum algorithms (such as QAOA and VQE) based on the Fubini-Study metric, in conjunction with Armijo line search rule.  Our results show that the new adaptive scheme outperforms QNGD and mitigates its dependence  on the optimal step size. Our scheme captures the optimal dynamical change in the step size as the training proceeds. In addition, our results also show that
a simple gradient descent scheme (with respect to the Euclidean metric) equipped with the adaptive scheme, can yield performances similar to the QNGD with optimal step size.

More importantly, the results presented provide further confirmation of the importance of differential geometry in variational quantum computations. This is not surprising given the fact that this type of computation is  optimization of smooth functions on Riemannian manifolds. More generally, we foresee that  advanced mathematics will feature centrally in  the NISQ era in guiding the design of faster and more efficient practical algorithms.

§ ACKNOWLEDGEMENT

We gratefully acknowledge the very helpful discussions with Nick Chancellor, and the rest of the Advanced Technology Team at QCI.

[1]
Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.
A quantum approximate optimization algorithm, 2014.

[2]
Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou,
Peter J. Love, Alán Aspuru-Guzik, and Jeremy L. O’Brien.
A variational eigenvalue solver on a photonic quantum processor.
Nature Communications, 5(1), Jul 2014.

[3]
Jarrod R McClean, Jonathan Romero, Ryan Babbush, and Alán Aspuru-Guzik.
The theory of variational hybrid quantum-classical algorithms.
New Journal of Physics, 18(2):023023, 2016.

[4]
Jarrod R McClean, Jonathan Romero, Ryan Babbush, and Alán Aspuru-Guzik.
The theory of variational hybrid quantum-classical algorithms.
New Journal of Physics, 18(2):023023, Feb 2016.

[5]
Michael Broughton, Guillaume Verdon, Trevor McCourt, Antonio J. Martinez,
Jae Hyeon Yoo, Sergei V. Isakov, Philip Massey, Murphy Yuezhen Niu, Ramin
Halavati, Evan Peters, Martin Leib, Andrea Skolik, Michael Streif, David Von
Dollen, Jarrod R. McClean, Sergio Boixo, Dave Bacon, Alan K. Ho, Hartmut
Neven, and Masoud Mohseni.
Tensorflow quantum: A software framework for quantum machine
learning, 2020.

[6]
James Stokes, Josh Izaac, Nathan Killoran, and Giuseppe Carleo.
Quantum natural gradient.
Quantum, 4:269, May 2020.

[7]
Shun-ichi Amari.
Neural learning in structured parameter spaces - natural riemannian
In M. C. Mozer, M. Jordan, and T. Petsche, editors, Advances in
Neural Information Processing Systems, volume 9. MIT Press, 1997.

[8]
Naoki Yamamoto.
On the natural gradient for variational quantum eigensolver, 2019.

[9]
David Wierichs, Christian Gogolin, and Michael Kastoryano.
Avoiding local minima in variational quantum eigensolvers with the
natural gradient optimizer.
Physical Review Research, 2(4), Nov 2020.

[10]
P. A. Absil, R. Mahony, and B. Andrews.
Convergence of the iterates of descent methods for analytic cost
SIAM Journal on Optimization, 16(2):531–547, 2005.

[11]
Germund G. Dahlquist.
A special stability problem for linear multistep methods.
BIT Numerical Mathematics, 3(1):27–43, 1963.

[12]
Jonathan Romero, Ryan Babbush, Jarrod R. McClean, Cornelius Hempel, Peter Love,
and Alán Aspuru-Guzik.
Strategies for quantum computing molecular energies using the unitary
coupled cluster ansatz, 2018.

[13]
Attila Szabo and Neil S. Ostlund.
Modern Quantum Chemistry: Introduction to Advanced Electronic
Structure Theory.
Dover Publications, Inc., Mineola, first edition, 1996.

[14]
Maria Schuld, Ville Bergholm, Christian Gogolin, Josh Izaac, and Nathan
Evaluating analytic gradients on quantum hardware.
Physical Review A, 99(3), Mar 2019.

[15]
Stephen L Campbell and Carl D Meyer.
Generalized inverses of linear transformations.
SIAM, 2009.

[16]
Gene H Golub, Per Christian Hansen, and Dianne P O'Leary.
Tikhonov regularization and total least squares.
SIAM journal on matrix analysis and applications,
21(1):185–194, 1999.

[17]
Abhinav Kandala, Antonio Mezzacapo, Kristan Temme, Maika Takita, Markus Brink,
Jerry M. Chow, and Jay M. Gambetta.
Hardware-efficient variational quantum eigensolver for small
molecules and quantum magnets.
Nature, 549(7671):242–246, Sep 2017.

§ APPENDIX: CONNECTION WITH $A-$STABILITY THEORY [11]

Let $f(\boldTheta) = \langle \varphi(\boldTheta)|H|\varphi(\boldTheta) \rangle$ be the expectation value where $|\varphi(\boldTheta)\rangle = U(\boldTheta) |0\rangle$ and $H$ is a Hamiltonian of interest (for instance, $U(\theta)$ is the QAOA ansatz and $H$ is a classical Hamiltonian i.e., a QUBO).  Let us also assume that $\boldTheta^*$ is a local minimum for the expectation value function $f(\boldTheta)$. One can  think of the dynamical update of the ansatz parameters  as the flow of the following dynamical system

d()d    =    - ∇f (()),

()    =    _0.

In the vicinity of $\boldTheta^*$, where $\boldTheta_0$ is assumed to belong,  we linearize the gradient flow above into the following linear ODE system (Dahlquist equations)

d()d    =    - H_f(^*) ·(-^*),

()    =    _0,

with $H_f(\boldTheta^*)$ the Hessian of $f$, which is positive semi-definite because $\boldTheta^*$ is a local minimum.  If we apply the forward Euler method, we get

\begin{equation}
\boldTheta_{k+1} = (Id -h H_f(\boldTheta_0) )\boldTheta_k =  (Id -h H_f(\boldTheta_0) )^k\boldTheta_0.
\end{equation}

The term $R = (Id -h H_f(\boldTheta_0) )$ is the stability function of the forward Euler method and is required to satisfy $||R||<1$  for the method to converge.  Staying in   the vicinity of a local minimum, the optimal rate for the forward Euler method is
$h = {\sum \lambda_i}/{\sum \lambda_i^2 },$  where   $\lambda_i$ are the eigenvalues of the Hessian matrix $H_f(\boldTheta_0).$ This comes from the fact that
the eigenvalues of the stability function $R$ are exactly the real numbers $1-h \lambda_i$, and subsequently, the optimal rate can be obtained with

\begin{equation}
argmin_{h>0} ||R|| = argmin_{h>0}  \sum_i (1-h \lambda_i)^2.
\end{equation}

Since Fubini-Study tensor is a positive semi-definite, after  the change of variables $\tau = F^{-1}(\boldTheta_0) \boldTheta$, the formula  above can also be applied to the new Hessian $H_f(\tau) = {F^{-1}}^t H_f F^{-1}$.
Notice this change of coordinates changes the Riemannian manifold from the Euclidean space to the Fubini-Study manifold, and the dynamical update (the gradient flow (<ref>)) changes accordingly.

Now consider depth one QAOA solving the max-cut problem for a triangle

\begin{eqnarray}\nonumber
U(\boldTheta)=U(\beta, \gamma) &=& e^{-i \Frac{\beta}{2} (X_1 + X_2 + X_3)}\\\nonumber
&& \times e^{-i \Frac{\gamma}{2} (1-Z_1Z_2 - Z_1Z_3 - Z_2Z_3)},   \nonumber
\end{eqnarray}

with $(\beta, \gamma)\in [0, \pi/2]\times [0, \pi]$.
Starting from the uniform superposition, we obtain the expectation value depicted in Figure <ref>, and the Fubini-Study tensor:

\begin{equation}
F(\boldTheta) = \left[\begin{array}{cc}
2-2 \cos \! \left(2 \boldTheta_{2}\right) & 0
\\
0 & \frac{1}{4}
\end{array}\right]
\end{equation}

A direct application of the optimal rate formula above gives $h = 0.093$ with $R=0.602$ when Fubini-Study is used,  and a slower rate $h = 0.022$ and $R=0.762$. In other words, Fubini-Study metric captures the correct geometry of the parameter space and makes the  gradient descent converge faster.

Expectation value landscape for the simple example of QAOA depth 1 solving the max cut problem for a triangle. The level sets around the two symmetric global minima are not exactly circular, which makes the gradient descent to slightly slowdown.


# The double hypergeometric series for the partition function of the 2D
anisotropic Ising model

G. M. Viswanathan  Department of Physics and National Institute of Science and
Technology of Complex Systems, Universidade Federal do Rio Grande do Norte,
59078-970 Natal–RN, Brazil

###### Abstract

In 1944 Lars Onsager published the exact partition function of the
ferromagnetic Ising model on the infinite square lattice in terms of a
definite integral. Only in the literature of the last decade, however, has it
been recast in terms of special functions. Until now all known formulas for
the partition function in terms of special functions have been restricted to
the important special case of the isotropic Ising model with symmetric
couplings. Indeed, the anisotropic model is more challenging because there are
two couplings and hence two reduced temperatures, one for each of the two axes
of the square lattice. Hence, standard special functions of one variable are
inadequate to the task. Here, we reformulate the partition function of the
anisotropic Ising model in terms of the Kampé de Fériet function, which is a
double hypergeometric function in two variables that is more general than the
Appell hypergeometric functions. Finally, we present hypergeometric formulas
for the generating function of multipolygons of given length on the infinite
square lattice, for isotropic as well as anisotropic edge weights. For the
isotropic case, the results allow easy calculation, to arbitrary order, of the
celebrated series first found by Cyril Domb.

##  1 Introduction

Lars Onsager solved the ferromagnetic Ising model with nearest neighbor
interactions on the infinite square lattice in 1942 and published the solution
in 1944 in a breakthrough paper [1]. Let $J_{1}$ and $J_{1}$ denote the
couplings, $K_{i}=J_{i}/kT$ the reduced temperatures, and $Z(K_{1},K_{2})$ the
partition function per site of the ferromagnetic Ising model on the square
lattice. Onsager’s solution is then given by

$\log(Z/2)=\frac{1}{2\pi^{2}}\underset{~{}0}{\overset{~{}~{}\mbox{\small$\pi$}}{\iint}}\log[\cosh
2K_{1}\cosh 2K_{2}-\sinh 2K_{1}\cos\theta_{1}-\sinh
2K_{2}\cos\theta_{2}]~{}d\theta_{1}d\theta_{2}~{}.~{}$ (1)

Furthermore, let

$\displaystyle 2\kappa_{1}$ $\displaystyle=$ $\displaystyle\tanh 2K_{1}~{}\rm
sech\,2K_{2}$ $\displaystyle 2\kappa_{2}$ $\displaystyle=$ $\displaystyle\tanh
2K_{2}~{}\rm sech\,2K_{1}~{}.$ (2)

Then, Onsager showed that

$\displaystyle\log Z-\tfrac{1}{2}\log(4\cosh 2K_{1}\cosh 2K_{2})$
$\displaystyle=\frac{1}{2\pi^{2}}\underset{0}{\overset{~{}~{}\mbox{\small$\pi$}}{\iint}}\log[1-2\kappa_{1}\cos\theta_{1}-2\kappa_{2}\cos\theta_{2}]~{}d\theta_{1}d\theta_{2}$
(3)
$\displaystyle=-\frac{1}{2}\sum_{j+k>0}^{\infty}{(2j+2k-1)!\over(j!)^{2}(k!)^{2}}\kappa_{1}^{2j}\kappa_{1}^{2k}~{}.$
(4)

This last Eq. (4) is Eq. (109b) in the 1944 paper [1]. It is easily obtained
from Eq. (3) as follows: (i) first Taylor expand the logarithm in the Mercator
series, (ii) then apply the binomial theorem to expand powers of
$(2\kappa_{1}\cos\theta_{1}+2\kappa_{2}\cos\theta_{2})$, and (iii) finally
perform the double integral of the now-separable integrands as an iterated
integral, term by term. For the case of the isotropic Ising model with
symmetric couplings $J_{1}=J_{2}=J$, we can write $K_{1}=K_{2}=K$ and
$\kappa_{1}=\kappa_{2}=\kappa$ and the above series simplifies considerably
(see Eq. (109c) of Ref. [1]).

In the literature of the last decade, the partition function for the isotropic
Ising model has been reformulated in terms of special functions (see below).
However, the anisotropic model is much more challenging because there are two
reduced temperatures — one for each coupling. Here we address this problem by
more carefully studying the series (4). We thereby obtain a double
hypergeometric reformulation of the partition function for the 2D anisotropic
Ising model.

Hypergeometric functions belong to the family of special functions [2, 3].
Much insight and understanding can be gained when a definite integral or an
infinite series is re-expressed in terms of special functions, because the
latter have been widely studied and are often better understood than the
former. The analogy with elementary functions can be helpful to make this
point clear:

$\displaystyle\sum_{n=1}^{\infty}{(1-x)^{n}\over n}$
$\displaystyle=-\log(x)~{},\quad\quad$ $\displaystyle|x|<1$
$\displaystyle\int_{a}^{b}{dx\over x}$ $\displaystyle=\log\left(b\over
a\right)~{},\quad\quad$ $\displaystyle a,b>0~{}.$

The left-hand sides of the above are precisely equal to the right-hand sides.
However, the appearance of the logarithm on the right gives us immediate
insight and intuition that may not be so obvious in the expressions to the
left. The same kind of insight can be gained from special functions. For
further details on the motivation for expressing definite integrals and
infinite series in terms of special functions, see Refs. [2, 6, 5, 3, 7, 4] .

In Section 2 we review the definition of the prerequisite hypergeometric
functions. In Section 3 we state our main results, together with the proofs.
Section 4 ends with discussion and conclusions.

##  2 The Kampé de Fériet double hypergeometric function

The generalized hypergeometric function ${}_{p}F_{q}$ of a variable $x$ has
series expansion $\sum c_{n}x^{n}$ such that the ratio $c_{n+1}/c_{n}$ of
successive coefficients is a ratio of polynomials in $n$. In other words, the
ratio of coefficients is a a rational function of $n$. The degrees of the
polynomials of the numerator and denominator are $p$ and $q+1$ respectively.
Let the Pochhammer symbol $(x)_{n}$ denote the rising factorial for
$n=1,2,3\ldots$ as follows:

$\displaystyle(x)_{0}$ $\displaystyle=1~{},$ $\displaystyle(x)_{n}$
$\displaystyle=(x)_{n-1}(x+n-1)$ (5)
$\displaystyle={\Gamma(x+n)\over\Gamma(x)}~{}.$

Then ${}_{p}F_{q}$ is defined as follows [2]:

$_{p}F_{q}\left[\begin{array}[]{c}{a_{1},a_{2},\dots,a_{p}}\\\
{b_{1},b_{2},\dots,b_{q}}\end{array};x\right]=\sum_{n=0}^{\infty}{(a_{1})_{n}(a_{2})_{n}\dots(a_{p})_{n}\over(b_{1})_{n}(b_{2})_{n}\dots(b_{q})_{n}}~{}{x^{n}\over{n!}}~{}.$
(6)

For $p<q+1$ the function is entire. In contrast, for $p>q+1$ it has zero
radius of convergence. For $p=q+1$ the radius of convergence is 1.

The isotropic Ising model can be expressed in terms of ${}_{p}F_{q}$ functions
(see below). In contrast, the expression (4) for the anisotropic Ising model
contains a power series in two different variables, hence we will need a more
general and more powerful hypergeometric function.

The Kampé de Fériet function $F^{p,q}_{r,s}$ is a double hypergeometric
function of higher order in 2 variables [8, 4] that generalizes the
${}_{p}F_{q}$ function. It is defined as

$\displaystyle F^{p,q}_{r,s}\left[~{}\begin{array}[]{c|c|}a_{1}\ldots
a_{p}&~{}b_{1},b^{\prime}_{1}~{};~{}\ldots~{};~{}b_{q},b^{\prime}_{q}\\\
c_{1}\ldots
c_{r}&~{}d_{1},d^{\prime}_{1}~{};~{}\ldots~{};~{}d_{s},d^{\prime}_{s}\end{array}~{}~{}x,y\right]$
(9)
$\displaystyle~{}~{}=~{}\sum_{m,n=0}^{\infty}{~{}\\!\prod_{k=1}^{p}(a_{k})_{m+n}~{}~{}\prod_{k=1}^{q}\big{[}(b_{k})_{m}(b^{\prime}_{k})_{n}\big{]}\over~{}\prod_{k=1}^{r}(c_{k})_{m+n}~{}~{}\prod_{k=1}^{s}\big{[}(d_{k})_{m}(d^{\prime}_{k})_{n}\big{]}}~{}{x^{m}y^{n}\over
m!n!}~{}~{}.$ (10)

We note in passing that for specific low values of the parameters
$p,~{}q,~{}r,~{}s$, the Kampé de Fériet function can be reduced to the much
simpler Appell hypergeometric functions, of which there are four. For the
particular values of the parameters $p$ and $q$ for the anisotropic Ising
model (see below), this standard reduction to the Appell functions is not
possible, so we do not further discuss this point here.

##  3 Results and proofs

Our first result is the natural continuation of previous results. In the
1970s, Glasser and Onsager obtained the following expression for the partition
function of the 2D isotropic Ising model [6, 9]:

$\ln Z(K)=\ln(2\cosh 2K)-{1\over
2}+{1\over\pi}{\mathsf{E}(4\kappa)}+\kappa^{2}~{}{{}_{4}}F_{3}\left[\begin{array}[]{c}{{1\over
2},1,1,{3\over 2}}\\\ {2,2,2}\end{array};16\kappa^{2}\right]~{}.$ (11)

where ${\mathsf{E}(\cdot)}$ denotes the complete elliptic integral of the
second kind that can also be written as a ${}_{2}F_{1}$ function [6, 2]. In
2015, a simpler hypergeometric formula was discovered [6]:

$\ln Z=\ln(2\cosh
2K)-\kappa^{2}~{}{{}_{4}}F_{3}\left[\begin{array}[]{c}{1,1,{3\over 2},{3\over
2}}\\\ {2,2,2}\end{array};16\kappa^{2}\right]~{}.$ (12)

The above two expressions are related via a hypergeometric identity (see Eq.
(33) in Ref. [6]). In 2016, the exact low temperature series was expressed as
an infinite series in terms of complete and partial Bell polynomials [7]. All
the above results are restricted to the isotropic case and do not apply to the
anisotropic Ising model. Here we address this limitation and report an advance
for the anisotropic case:

###### Theorem 1

Let $\kappa_{1}$ and $\kappa_{2}$ be given by (2). Then the partition function
$Z(K_{1},K_{2})$ of the anisotropic Ising model as defined in (1) is given by

$\displaystyle\log Z-\tfrac{1}{2}\log(4\cosh 2K_{1}\cosh 2K_{2})$
$\displaystyle=-3~{}\kappa_{1}^{2}\kappa_{2}^{2}~{}{F^{2,1}_{0,2}\left[~{}\begin{array}[]{c|c|}2,5/2&~{}1,1{}\\\
\mbox{--}&~{}2,2;2,2\end{array}~{}~{}4\kappa^{2}_{1},4\kappa^{2}_{2}~{}\right]}$
(15)
$\displaystyle~{}~{}~{}-\frac{1}{2}~{}\kappa_{1}^{2}~{}_{3}F_{2}\left[\begin{array}[]{c}{1,1,3/2}\\\
2,2\end{array};4\kappa_{1}^{2}~{}\right]~{}-\frac{1}{2}\kappa_{2}^{2}~{}_{3}F_{2}\left[\begin{array}[]{c}{1,1,3/2}\\\
2,2\end{array};4\kappa_{2}^{2}~{}\right]~{}.$ (20)

###### Corollary 1

$\displaystyle\log Z-\tfrac{1}{2}\log(4\cosh 2K_{1}\cosh 2K_{2})$
$\displaystyle=-3~{}\kappa_{1}^{2}\kappa_{2}^{2}~{}{F^{2,1}_{0,2}\left[~{}\begin{array}[]{c|c|}2,5/2&~{}1,1{}\\\
\mbox{--}&~{}2,2;2,2\end{array}~{}~{}4\kappa^{2}_{1},4\kappa^{2}_{2}~{}\right]}$
(23)
$\displaystyle~{}~{}~{}+\frac{1}{2}\left[\,\log\left(\frac{1}{2}\left(1+\sqrt{1-4\kappa_{1}^{2}}\right)\right)+\log\left(\frac{1}{2}\left(1+\sqrt{1-4\kappa_{2}^{2}}\right)\right)\,\right]~{}.$
(24)

To prove the above, we will use the following supporting lemmas:

######  Lemma 1 (Basic identities for rising factorials)

Let $(a)_{n}$ be defined according to (5). Then,

$\displaystyle n!$ $\displaystyle=~{}$ $\displaystyle(1)_{n}~{},$ (25)
$\displaystyle(n+1)!$ $\displaystyle=~{}$ $\displaystyle(2)_{n}~{},$ (26)
$\displaystyle a+n~{}$ $\displaystyle=~{}$ $\displaystyle
a{(a+1)_{n}\over(a)_{n}}~{},$ (27) $\displaystyle(2n+1)!~{}$
$\displaystyle=~{}$ $\displaystyle 2^{2n}(3/2)_{n}(1)_{n}~{},$ (28)
$\displaystyle(2j+2k+3)!~{}$ $\displaystyle=~{}$ $\displaystyle{6\cdot
2^{2(j+k)}(2)_{j+k}(5/2)_{j+k}}~{}.$ (29)

’

  * Proof

The first three above follow immediately from the definition of the rising
factorial. The claim (28) can be proven by splitting the terms in the
factorial into odd and even factors and dividing each factor by 2:

$\displaystyle(2n+1)!$ $\displaystyle=$
$\displaystyle[1]\times(3)(5)(7)\ldots(2n+1)\times(2)(4)(6)\ldots(2n)$ (30)
$\displaystyle=$ $\displaystyle 2^{2n}~{}(3/2)_{n}(1)_{n}~{}.$

The claim (29) is slightly more complicated but follows from the previous 4
claims:

$\displaystyle(2j+2k+3)!$ $\displaystyle=(2j+2k+1)!~{}(2j+2k+2)(2j+2k+3)$ (33)
$\displaystyle=2^{2j+2k}(3/2)_{j+k}(1)_{j+k}~{}(2j+2k+2)(2j+2k+3)$
$\displaystyle=2^{2j+2k+2}(3/2)_{j+k}(1)_{j+k}(j+k+1)(j+k+3/2)$
$\displaystyle={2^{2j+2k+2}(3/2)_{j+k}(1)_{j+k}(2)_{j+k}(3/2)(5/2)_{j+k}\over(1)_{j+k}(3/2)_{j+k}}$
$\displaystyle={6\cdot 2^{2(j+k)}(2)_{j+k}(5/2)_{j+k}}~{},$

which completes the proof. $\square$

Our next lemma deals with the infinite double series in (4). For convenience
let us write the partition function in (4) as

$\log Z-\tfrac{1}{2}\log(4\cosh 2K_{1}\cosh
2K_{2})=-\frac{1}{2}\Phi(\kappa_{1},\kappa_{2})~{},$ (34)

with $\Phi$ thus defined as

$\Phi(x_{1},x_{2})=\sum_{j+k>0}^{\infty}{(2j+2k-1)!\over(j!)^{2}(k!)^{2}}x_{1}^{2j}x_{1}^{2k}~{}.$
(35)

###### Lemma 2

$\displaystyle\Phi(x_{1},x_{2})$ $\displaystyle=$ $\displaystyle
6~{}x_{1}^{2}x_{2}^{2}~{}{F^{2,1}_{0,2}\left[~{}\begin{array}[]{c|c|}2,5/2&~{}1,1{}\\\
\mbox{--}&~{}2,2;2,2\end{array}~{}~{}4x^{2}_{1},4x^{2}_{2}~{}\right]}$ (43)
$\displaystyle+x_{1}^{2}~{}_{3}F_{2}\left[\begin{array}[]{c}{1,1,3/2}\\\
2,2\end{array};4x_{1}^{2}~{}\right]+x_{2}^{2}~{}_{3}F_{2}\left[\begin{array}[]{c}{1,1,3/2}\\\
2,2\end{array};4x_{2}^{2}~{}\right]~{}.$

  * Proof

If the sum in (35) started from $j=k=0$, then the task would be easier, but
the sum excludes the point $j=k=0$, which is problematic. So we will rewrite
$\Phi$ as a sum of a double series for $j,k>1$ and separate single series
corresponding to $j=0$ and $k=0$:

$\displaystyle\Phi(x_{1},x_{2})$ $\displaystyle=$
$\displaystyle\sum_{j,k=1}^{\infty}{(2j+2k-1)!\over(j!)^{2}(k!)^{2}}x_{1}^{2j}x_{1}^{2k}+\sum_{j=1}^{\infty}{(2j-1)!\over(j!)^{2}}(x_{1}^{2j}+x_{2}^{2j})$
(45)

We next change the summation indices to start at zero and invoke the
identities in lemma 1:

$\displaystyle\Phi(x_{1},x_{2})$ $\displaystyle=$
$\displaystyle\sum_{j,k=0}^{\infty}{(2(j+k))+3)!\over((j+1)!)^{2}((k+1)!)^{2}}x_{1}^{2(j+1)}x_{1}^{2(k+1)}$
(46)
$\displaystyle+\sum_{j=0}^{\infty}{(2j+1)!\over((j+1)!)^{2}}(x_{1}^{2(j+1)}+x_{2}^{2(j+1)})$
$\displaystyle=$ $\displaystyle x_{1}^{2}x_{2}^{2}\sum_{j,k=0}^{\infty}{6\cdot
2^{2(j+k)}~{}(2)_{j+k}~{}(5/2)_{j+k}\over((2)_{j})^{2}~{}((2)_{k})^{2}}x_{1}^{2j}x_{1}^{2k}$
$\displaystyle+x_{1}^{2}\sum_{j=0}^{\infty}{2^{2j}(3/2)_{j}(1)_{j}\over((2)_{j})^{2}}x_{1}^{2j}$
$\displaystyle+x_{2}^{2}\sum_{j=0}^{\infty}{2^{2j}(3/2)_{j}(1)_{j}\over((2)_{j})^{2}}x_{2}^{2j}~{}.$

Finally, from the definitions (6) of the generalized hypergeometric function
and (10) of the Kampé de Fériet function, (46) becomes

$\displaystyle{\Phi(x_{1},x_{2})}$ $\displaystyle=$ $\displaystyle
6~{}x_{1}^{2}x_{2}^{2}~{}{F^{2,1}_{0,2}\left[~{}\begin{array}[]{c|c|}2,5/2&~{}1,1{}\\\
\mbox{--}&~{}2,2;2,2\end{array}~{}~{}4x^{2}_{1},4x^{2}_{2}~{}\right]}$ (54)
$\displaystyle+~{}x_{1}^{2}~{}_{3}F_{2}\left[\begin{array}[]{c}{1,1,3/2}\\\
2,2\end{array};4x_{1}^{2}~{}\right]+~{}x_{2}^{2}~{}_{3}F_{2}\left[\begin{array}[]{c}{1,1,3/2}\\\
2,2\end{array};4x_{x}^{2}~{}\right]~{},~{}$

which completes the proof. $\square$

  * Proof of theorem 1

The claim follows from using lemma 2 in (34). $\square$

  * Proof of corollary 1

To prove the claim (24) starting from theorem 1, it suffices to show that

$x~{}_{3}F_{2}\left[\begin{array}[]{c}{1,1,3/2}\\\
2,2\end{array};x~{}\right]=-{4\log\left(\frac{1}{2}\left(\sqrt{1-x}+1\right)\right)}~{}.$
(56)

We begin with the Mercator series and Newton’s generalized binomial theorem:

$\displaystyle\log(1+\sqrt{1-x})$ $\displaystyle=$
$\displaystyle-\sum_{n=1}^{\infty}{(-1)^{n}(1-x)^{n/2}\over n}$ (57)
$\displaystyle=$ $\displaystyle-\sum_{n=1}^{\infty}{(-1)^{n}\over
n}\sum_{j=0}^{\infty}{n/2\choose j}(-x)^{j}$ (58) $\displaystyle=$
$\displaystyle-\sum_{j=0}^{\infty}(-x)^{j}\sum_{n=1}^{\infty}{(-1)^{n}\over
n}{n/2\choose j}~{}.$ (59)

Next, we rewrite the binomial coefficient more conveniently as

${n/2\choose j}={(n/2+1-j)_{j}\over j!}~{}.$ (60)

We then get

$\displaystyle\log(1+\sqrt{1-x})$ (61) $\displaystyle=$
$\displaystyle-\sum_{j=0}^{\infty}(-x)^{j}\sum_{n=1}^{\infty}{(-1)^{n}\over
n}{(n/2+1-j)_{j}\over j!}$ (62) $\displaystyle=$
$\displaystyle-\sum_{j=0}^{\infty}{(-x)^{j}\over
j!}\sum_{n=1}^{\infty}{(-1)^{n}\over n}{(n/2+1-j)_{j}}$ (63) $\displaystyle=$
$\displaystyle-\sum_{j=0}^{\infty}{(-x)^{j}\over
j!}\sum_{n=1}^{\infty}\left[{(n+1-j)_{j}\over 2n}-{(n+1/2-j)_{j}\over
2n-1}\right]~{}.~{}~{}~{}~{}~{}~{}~{}~{}~{}$ (64)

The sum over $n$ can be done after transforming the rising factorials to gamma
functions:

$\displaystyle\sum_{n=1}^{\infty}\left[{(n+1-j)_{j}\over
2n}-{(n+1/2-j)_{j}\over 2n-1}\right]$ (65)
$\displaystyle=\left[\frac{\,\Gamma(-j)}{2\,\Gamma(1-j)^{2}}-\frac{\sqrt{\pi}\,\Gamma(-j)}{2\,\Gamma\left(\frac{1}{2}-j\right)\,\Gamma(1-j)}\right]$
(66) $\displaystyle=\left[\frac{\sqrt{\pi}\cos(\pi
j)\,\Gamma\left(j+\frac{1}{2}\right)-\sin(\pi j)\,\Gamma(j)}{2\pi
j}\right]~{}.$ (67)

where in the last step we have used Euler’s reflection formula. The terms for
$j>0$ in (64) can be evaluated as usual, but the term for $j=0$ is
indeterminate due to the $j$ in the denominator in (67), so it needs to be
evaluated as a limit. We take the limit $z\to j=0$ as follows:

$\lim_{z\to 0}\frac{\sqrt{\pi}\cos(\pi
z)\,\Gamma\left(z+\frac{1}{2}\right)-\sin(\pi z)\,\Gamma(z)}{2\pi z}=-\log
2~{}.$ (68)

For $j\neq 0$ the term with $\sin(\pi j)$ vanishes. We can then simplify the
expression using Legendre’s duplication formula for the gamma function and
applying lemma 1:

$\displaystyle\log(1+\sqrt{1-x})$ $\displaystyle=$ $\displaystyle\log
2-\sum_{j=1}^{\infty}{1\over j!}\frac{\sqrt{\pi}\cos(\pi
j)\,\Gamma\left(j+\frac{1}{2}\right)}{2\pi j}(-x)^{j}$ (69)
$\displaystyle\log(\tfrac{1}{2}(1+\sqrt{1-x}))$ $\displaystyle=$
$\displaystyle-\sum_{j=1}^{\infty}{1\over
j!}\frac{\sqrt{\pi}\,\Gamma\left(j+\frac{1}{2}\right)}{2\pi j}x^{j}$ (70)
$\displaystyle=$
$\displaystyle-\sum_{j=1}^{\infty}\frac{2^{1-2j}\,\Gamma(2j)}{2j\,\Gamma(j)j!}x^{j}$
(71) $\displaystyle=$
$\displaystyle-x\sum_{j=0}^{\infty}\frac{2^{-2j-2}\,\Gamma(2j+2)}{(j+1)!^{2}}x^{j}$
(72) $\displaystyle=$
$\displaystyle-\frac{x}{4}\sum_{j=0}^{\infty}\frac{2^{-2j}(2j+1)!}{(j+1)!^{2}}x^{j}$
(73) $\displaystyle=$
$\displaystyle-\frac{x}{4}\sum_{j=0}^{\infty}\frac{(3/2)_{n}(1)_{n}}{(2)_{n}^{2}}x^{j}$
(77)
$\displaystyle=-\frac{x}{4}~{}{{}_{3}F_{2}}\left[\begin{array}[]{c}{1,1,3/2}\\\
2,2\end{array};x~{}\right]~{},$

thereby establishing (56), which completes the proof. $\square$

Theorem 1 allows the generating function of multipolygons on the infinite
square lattice to be reformulated in terms of hypergeometric functions. By
multipolygon on the infinite square lattice is meant a connected or
disconnected graph all of whose nodes have even degrees. Let $t_{1}$ be the
weight of edges in the first lattice axis direction and $t_{2}$ be the edge
weight for the other direction. For convenience, we will say that $t_{1}$ is
the weight of edges of the first type (say, horizontal edges) and $t_{2}$ the
weight of edges of the second type (say, vertical edges).

Let $\Lambda_{N}(t_{1},t_{2})$ be the generating function for the the number
of multipolygons on a finite square lattice with $N=M^{2}$ sites, such that
the coefficient of the term of degree $t_{1}^{m}t_{2}^{n}$ gives the number of
multipolygons with $m$ edges of the first type and $n$ edges of the second
type. The case $t_{1}=t_{2}$, which we will call the isotropic case, has been
widely studied [10]. The anisotropic case $t_{1}\neq t_{2}$ has received
somewhat less attention [11].

We define the generating function $\Lambda(t_{1},t_{2})$ for anisotropic
multipolygons on the infinite square lattice in terms of a limit when the
number of sites goes to infinity:

$\Lambda(t_{1},t_{2})=\lim_{N\to\infty}[\Lambda_{N}(t_{1},t_{2})]^{1/N}~{}.$
(78)

The connection between the function
$\Lambda(t_{1},t_{2})=1+t_{1}^{2}t_{2}^{2}+t_{1}^{2}t_{2}^{4}+t_{1}^{4}t_{2}^{2}\ldots$
and the low-temperature series for the partition function of the 2D Ising
model is well known. Moreover, the self-dual property of the Ising model on
the square lattice leads to a similar relation for the high-temperature
series. Let

$\displaystyle u_{i}$ $\displaystyle=$ $\displaystyle\tanh(K_{i})~{},$ (79)
$\displaystyle v_{i}$ $\displaystyle=$ $\displaystyle\exp(-2K_{i})~{},$ (80)

for $i=1,2$ be high- and low-temperature variables respectively. Then the
following are well known [10, 11]:

$\displaystyle Z$ $\displaystyle=$ $\displaystyle 2\cosh K_{1}\cosh
K_{2}~{}\Lambda(u_{1},u_{2})$ (81) $\displaystyle Z$ $\displaystyle=$
$\displaystyle\exp(K_{1})\exp(K_{2})~{}\Lambda(v_{1},v_{2})~{}.$ (82)

From either of these two expressions, we can obtain the double hypergeometric
formulation of $\Lambda(t_{1},t_{2})$:

###### Theorem 2

Let $\Lambda(t_{1},t_{2})$ be the generating function for multipolygons with
anisotropic weights $t_{1}$ and $t_{2}$ for the two directions on the infinite
square lattice, defined according to (78). Let $\Phi$ be given by (54). Then,

$\displaystyle\log\Lambda(t_{1},t_{2})$ $\displaystyle=$
$\displaystyle~{}\frac{1}{2}\log\left[{\left({t_{1}}^{2}+1\right)\left({t_{2}}^{2}+1\right)}\right]$
(83)
$\displaystyle-\frac{1}{2}\Phi\left[\frac{{t_{1}}\left(1-{t_{2}}^{2}\right)}{\left({t_{1}}^{2}+1\right)\left({t_{2}}^{2}+1\right)},\frac{{t_{2}}\left(1-{t_{1}}^{2}\right)}{\left({t_{1}}^{2}+1\right)\left({t_{2}}^{2}+1\right)}\right]~{}.$

  * Proof

We can arrive at the claim through either the low- or the high-temperature
formulations, due to the self-dual property. For completeness, we show both
approaches. First note that

${\sqrt{4\cosh 2K_{1}\cosh 2K_{2}}\over 2\cosh K1\cosh
K1}=\sqrt{\left({u_{1}}^{2}+1\right)\left({u_{2}}^{2}+1\right)}~{}.$ (84)

Similarly,

${\sqrt{4\cosh 2K_{1}\cosh
2K_{2}}\over\exp(2K_{1})\exp(2K_{1})}=\sqrt{\left({v_{1}}^{2}+1\right)\left({v_{2}}^{2}+1\right)}~{}.$
(85)

Next observe that $\kappa_{1}$ and $\kappa_{2}$ can be expressed in terms of
the high- and low-temperature variables as

$\displaystyle\kappa_{1}=\frac{{u_{1}}\left(1-{u_{2}}^{2}\right)}{\left({u_{1}}^{2}+1\right)\left({u_{2}}^{2}+1\right)}=\frac{{v_{2}}\left({v_{1}}^{2}-1\right)}{\left({v_{1}}^{2}+1\right)\left({v_{2}}^{2}+1\right)}~{},$
(86)
$\displaystyle\kappa_{2}=\frac{{u_{2}}\left(1-{u_{1}}^{2}\right)}{\left({u_{2}}^{2}+1\right)\left({u_{1}}^{2}+1\right)}=\frac{{v_{1}}\left({v_{2}}^{2}-1\right)}{\left({v_{2}}^{2}+1\right)\left({v_{1}}^{2}+1\right)}~{}.$
(87)

The claim (83) follows from substituting these expressions for $\kappa_{i}$
into (34) and then using either (81) and (84) for the high-temperature
variable or else (82) with (85) for the low-temperature variable. $\square$

The above result simplifies considerably for the isotropic case
$t_{1}=t_{2}=t$. The next result is known to experts in the field, but for
some reason never seems to have been published as such:

###### Corollary 2

Let $\Lambda(t)$ denote the generating function $\Lambda(t,t)$ of
multipolygons on the infinite square lattice with isotropic weights $t$. Then,

$\log\Lambda(t)=\log\left({t}^{2}+1\right)-\frac{t^{2}\left(t^{2}-1\right)^{2}}{\left(t^{2}+1\right)^{4}}~{}{{}_{4}}F_{3}\left[\begin{array}[]{c}{1,1,{3\over
2},{3\over 2}}\\\
{2,2,2}\end{array};\frac{16t^{2}\left(t^{2}-1\right)^{2}}{\left(t^{2}+1\right)^{4}}\right]~{}.$
(88)

  * Proof

The claim follows from (12), (34) and theorem 2 with $t_{1}=t_{2}=t$.
$\square$

##  4 Discussion and conclusion

We first comment on a possible avenue for further research. There is a deep
connection between the isotropic Ising model on the one hand and spanning
trees on the infinite square lattice on the other hand. Guttmann and Rogers
have defined a generating function that generalizes the spanning tree constant
[12]. Subsequently, it was shown that this spanning tree generating function
is related to the partition function of the isotropic Ising model in a precise
way [13]. The connection is due to the fact that the same Mahler measure
appears in both, which in turn is due to the random walk structure function
for the square lattice. A natural question now arises. Let us assume that one
allows the spanning tree generating function to have anisotropic weights in
some suitably defined manner. Then, does the known connection between spanning
trees and the isotropic Ising model generalize to the anisotropic model? We
tentatively believe that the answer is yes, and hope that others will take
interest in finding a conclusive answer.

We also comment on how the results were originally intuited prior to the
discovery of the actual proofs. The results above were not first proven or
even intuited using computer algebra systems, theorem-proving software,
artificial intelligence, machine learning, and so forth. Wolfram Mathematica
software, for example, cannot evaluate the integral in (3) and is unable to
recognize the series for the definition $\Phi$ in (35). In fact, the Kampé de
Fériet function (10) is not part of the repertoire of special functions
included in Mathematica (at the time of this writing). Instead, theorem 1 was
the result of traditional mathematical detective work, starting from (4).
After re-expressing the infinite double series in the form (35) of $\Phi$, we
calculated the ratio of successive terms of the series, in the two variables
separately and also together. We thus saw that these ratios are rational
functions of the degrees of the two variables — precisely the mathematical
signature of hypergeometric series. It was then just a matter of finding the
suitable hypergeometric function. We first looked at the Appell functions but
they are not of high enough order. The Kampé de Fériet function was the next
natural candidate. What made the task slightly more difficult is that both the
Appell and Kampé de Fériet functions are double hypergeometric and therefore
not found in the textbooks. For example, the Special Functions by Andrews,
Askey and Roy [2] and Special Functions: A Graduate Text by Beals and Wong [3]
are both considered to be authoritative texts, but neither book deals with the
Appell or Kampé de Fériet functions. Once the relevant definitions were found,
however, the proofs came quite naturally.

Finally, we note that theorem 2 and corollary 2 allow very easy explicit
evaluation of the multipolygon generating function to arbitrary order. For the
anisotropic case we get from (83)

$\displaystyle\Lambda(t_{1},t_{2})$ $\displaystyle=$ $\displaystyle
1~{}+~{}t_{1}^{2}t_{2}^{2}~{}+~{}t_{1}^{4}t_{2}^{2}~{}+~{}t_{1}^{2}t_{2}^{4}$
(89)
$\displaystyle~{}+~{}t_{1}^{2}t_{2}^{6}~{}+~{}3t_{1}^{4}t_{2}^{4}~{}+~{}t_{1}^{6}t_{2}^{2}$
$\displaystyle~{}+~{}t_{1}^{2}t_{2}^{8}~{}+~{}6t_{1}^{4}t_{2}^{6}~{}+~{}6t_{1}^{6}t_{2}^{4}~{}+~{}t_{1}^{8}t_{2}^{2}$
$\displaystyle~{}+~{}t_{1}^{2}t_{2}^{10}~{}+~{}10t_{1}^{4}t_{2}^{8}~{}+~{}22t_{1}^{6}t_{2}^{6}~{}+~{}10t_{1}^{8}t_{2}^{4}~{}+~{}t_{1}^{10}t_{2}^{2}+\ldots$

For the isotropic case, the formula (88) generates, to arbitrary order, the
celebrated series for the generating function of multipolygons first found by
Cyril Domb [14] in 1949 and whose coefficients have been incorporated into the
On-Line Encyclopedia of Integer Sequences (OEIS) [15]:

$\displaystyle\Lambda(t)$ $\displaystyle=$ $\displaystyle
1+t^{4}+2t^{6}+5t^{8}+14t^{10}+44t^{12}+152t^{14}+566t^{16}$ (90)
$\displaystyle+2234t^{18}+9228t^{20}+39520t^{22}+174271t^{24}+787246t^{26}$
$\displaystyle+3628992t^{28}+17019374t^{30}+O(t^{32})~{}.$

This work was supported by CNPq (grant no. 302051/2018-0). We thank Wolfram
Research and Daniel Pryjma for a courtesy license for Mathematica. We thank H.
D. Jennings for helpful comments.

## References

## References

  * [1] Onsager L 1944 Crystal Statistics. I. A Two-Dimensional Model with an Order-Disorder Transition, Phys. Rev.  65 117
  * [2] Andrews G E, Askey R, and Roy R 1999 Special functions (Encyclopedia of Mathematics and Its Applications vol 71) ed G C Rota (Cambridge: Cambridge University Press)
  * [3] Beals R and Wong R 2010 Special Functions: A Graduate Text (Cambridge: Cambridge University Press)
  * [4] Exton, H 1978 Handbook of hypergeometric integrals, Mathematics and its Applications (Chichester: Ellis Horwood Ltd.)
  * [5] Cipra B A 1998 A new testament for special functions? SIAM News 31(2)
  * [6] Viswanathan G M, 2015 The hypergeometric series for the partition function of the 2D Ising model J. Stat. Mech. 2015 P04007
  * [7] Siudem G, Fronczak A and Fronczak P 2016 Exact low-temperature series expansion for the partition function of the zero-field Ising model on the infinite square lattice Sci. Rep. 6 33523
  * [8] Kampé de Fériet, J 1937 La fonction hypergéométrique Mémorial des sciences mathématiques 85 (Paris: Gauthier-Villars)
  * [9] Glasser M L, personal e-mail communication, 12 November 2014
  * [10] Feynman R P, 1972 Statistical Mechanics. A set of lectures (Reading: Benjamin and Cummings Publishing)
  * [11] Newell G F and Montroll E W 1953 On the Theory of the Ising Model of Ferromagnetism, Rev. Mod. Phys. 25 353
  * [12] Guttmann A J and Rogers M D 2012 Spanning tree generating functions and Mahler measures J. Phys. A: Math. Theor. 45 494001
  * [13] Viswanathan G M 2017 Correspondence between spanning trees and the Ising model on a square lattice Phys. Rev. E 95 062138
  * [14] Domb C 1949 Order-disorder statistics. ii. a two-dimensional model Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 199 199
  * [15] OEIS Foundation Inc. 2019 The On-Line Encyclopedia of Integer Sequences, https://oeis.org/A002890