Determining system Hamiltonian from eigenstate measurements without correlation functions

TL;DR

This paper introduces a novel algorithm that reconstructs local Hamiltonians solely from local eigenstate measurements, bypassing the need for nonlocal correlation data, and demonstrates its effectiveness through numerical tests.

Contribution

It presents a new method to determine local Hamiltonians using only local measurements, reformulating the problem as an optimization task and employing machine learning techniques.

Findings

Successful numerical reconstruction of random local Hamiltonians.

Local measurements suffice to recover the full Hamiltonian in generic cases.

The method avoids reliance on nonlocal correlation functions.

Abstract

Local Hamiltonians arise naturally in physical systems. Despite its seemingly `simple' local structure, exotic features such as nonlocal correlations and topological orders exhibit in eigenstates of these systems. Previous studies for recovering local Hamiltonians from measurements on an eigenstate $|\psi\rangle$ require information of nonlocal correlation functions. In this work, we develop an algorithm to determine local Hamiltonians from only local measurements on $|\psi\rangle$, by reformulating the task as an unconstrained optimization problem of certain target function of Hamiltonian parameters, with only polynomial number of parameters in terms of system size. We also develop a machine learning-based-method to solve the first-order gradient used in the algorithm. Our method is tested numerically for randomly generated local Hamiltonians and returns promising reconstruction in the desired accuracy. Our result shed light on the fundamental question on how a single eigenstate can encode the full system Hamiltonian, indicating a somewhat surprising answer that only local measurements are enough without additional assumptions, for generic cases.

Figures

Tables

Table 1. Table 1: Time consumed of each trial and the number of trials. For general operator cases and the 4-qubit local operators case, we used 200 examples to achieve the averages. For 7-qubit local operators, The result is obtained from 40 examples.

	General operators				Local operators
Number of qubits	4	5	6	7	4	7
Time for single trial $t_0$(s)	0.0487	0.0518	0.0539	0.0579	0.2246	7.938
Number of trials $\overline{n}$	32.36	42.34	45.89	161.6	54.84	539.4
Average time cost $\left(\overline{t}\right)$ (s)	1.576	2.193	2.377	2.4735	12.32	4282

Table 2. Table 2: The types of the variables and values of the intermediate functions.

Node	Function	Derivative	Type of variable(s)	Type of function value
$v_2$	$v_2=\sum c_i\left(A_i-a_{i}I\right)$	$A_i-a_{i}I$	number	matrix
$v_3$	$v_3=v_2\cdot v_2$	$\frac{\partial v_2}{v_{1,k}}\tilde{H}+\tilde{H}\frac{\partial v_2}{v_{1,k}}$	matrix	matrix
$v_4$	$v_4=e^{-v_3}$	${\int }_0^1{e}^{\alpha v_3}\frac{\partial v_3}{v_{1,k}}{e}^{\left(1-\alpha \right)v_3}𝑑\alpha$	matrix	matrix
$v_5$	$v_5=\text{tr}{v}_4$	$\text{tr}\left(\frac{\partial v_4}{v_{1,k}}\right)$	matrix	number
$v_6$	$v_6=\frac{v_4}{v_5}$	$\frac{\partial v_4}{v_{1,k}}\frac{1}{v_5}-\frac{v_4}{v_5^2}\frac{\partial v_5}{\partial v_{1,k}}$	matrix & number	matrix
$v_7$	$v_7,j=\text{tr}\left(A_j{v}_6\right)-a_j$	$\text{tr}\left(A_j\frac{\partial v_6}{\partial v_{1,k}}\right)$	matrix	number
$v_8$	$v_8=\sum v_{7,j}^2$	$2v_{7,j}\frac{\partial v_{7,j}}{\partial v_{1,i}}$	number	number
$v_9$	$v_9=\text{tr}\left(v_3{v}_6\right)$	$\text{tr}\left(\frac{\partial v3}{\partial v_{1,i}}{v}_6+v_3\frac{\partial v6}{\partial v_{1,i}}\right)$	matrix	number
$v_{10}$	$v_{1}0=v_8+v_9$	$\frac{\partial v8}{\partial v_{1,i}}+\frac{\partial v9}{\partial v_{1,i}}$	number	number

Equations

$$|\psi\rangle\mbox{$\langle\psi|$}=\frac{e^{-\beta H}}{\textrm{tr}(e^{-\beta H})}$$

$$\tilde{H}=H-\lambda I=\sum _i{c}_i(A_i-a_{i}I),$$

$$\mbox{$\langle\psi|$}H|\psi\rangle=\lambda=\sum_{i}c_{i}\mbox{$\langle\psi|$}A_{i}|\psi\rangle=\sum_{i}c_{i}a_{i}.$$

$$\rho(\vec{c})=|\psi\rangle\mbox{$\langle\psi|$}=\frac{e^{-\beta\tilde{H}^{2}}}{\textrm{tr}(e^{-\beta\tilde{H}^{2}})}$$

$$\text{tr}(A_i\rho (\vec{c}))=a_i,\text{tr}({\tilde{H}}^2\rho (\vec{c}))=0.$$

$$f(\vec{x})=\sum _{i=1}^m[\text{tr}(A_i\rho (\vec{x}))-a_i{]}^2+\text{tr}({\tilde{H}}^2\rho (\vec{x})),$$

$$\rho (\vec{c})=\frac{e^{-{\tilde{H}}_{\beta }^2}}{\text{tr}(e^{-{\tilde{H}}_{\beta }^2})}.$$

$$f(\vec{x})=\sum _{i=1}^m[\text{tr}(A_i\cdot \rho (\vec{x})]-a_i{)}^2+\text{tr}({\tilde{H}}_{\beta }^2\cdot \rho (\vec{x})).$$

$$f(H_{al},H_{rd})=\frac{\text{tr}{H}_{al}{H}_{rd}}{\sqrt{\mathrm{Tr}{H}_{al}^2}\sqrt{\mathrm{Tr}{H}_{rd}^2}}.$$

$$\text{tr}(A_i{A}_j)=d{\delta }_{ij},$$

$$f(H_{al},H_{rd})=\frac{\mathrm{Tr}{H}_{al}{H}_{rd}}{\sqrt{\mathrm{Tr}{H}_{al}^2}\sqrt{\mathrm{Tr}{H}_{rd}^2}}=\frac{\vec{c}\cdot {\vec{c}}^{\prime }}{\parallel \vec{c}\parallel \parallel {\vec{c}}^{\prime }\parallel },$$

$$H=c_1{A}_1+c_2{A}_2+c_3{A}_3.$$

$$H=\sum _{1\le j,i\le 4}\sum _{i<j}{c}_{ij}{A}_{ij},$$

$$H=\sum _{1\le i\le 6}{c}_{i,i+1}{A}_{i,i+1},$$

$$\bar{t}=\bar{n}{t}_0.$$

$$H=\sum _i{c}_i{L}_i,$$

$$M_{ij}=\frac{1}{2}⟨\psi |\{L_i,L_j\}|\psi ⟩-⟨\psi |L_i|\psi ⟩⟨\psi |L_j|\psi ⟩,$$

$$H_{cor}=\sum {\omega }_i{L}_i$$

$$H=\sum _{i=1}^4{a}_i{\sigma }_z^i+\sum _{i=1}^3{b}_i{\sigma }_x^i{\sigma }_x^{i+1},$$

$$⟨\psi |\{L_i,L_j\}|\psi ⟩,$$

$$P(E=E_g)=\frac{e^{-E_g}}{\mathrm{Tr}{e}^{-{\tilde{H}}_{\beta }^2}}$$

$$P(E=E_g)=\frac{1}{\mathrm{Tr}{e}^{-{\tilde{H}}_{\beta }^2}}.$$

$$\frac{\partial e^{X(c_k)}}{\partial c_k}={\int }_0^1{e}^{\alpha X(c_k)}\frac{\partial X(c_k)}{\partial c_k}{e}^{(1-\alpha )X(c_k)}d\alpha ,$$

$$A(\beta )=\frac{\partial e^{X(c_k)}}{\partial c_k}={\int }_0^{\beta }{e}^{\alpha X(c_k)}\frac{\partial X(c_k)}{\partial c_k}{e}^{-\alpha X(c_k)}d\alpha ,$$

$$B(\beta )=e^{-\beta X(c_k)}A(\beta ),$$

$$C(\beta )=e^{-\beta X(c_k)}.$$

$$\frac{d}{d\beta }\left(\begin{array}{c}B(\beta )\\ C(\beta )\end{array}\right)=\left(\begin{array}{cc}-X(c_k)& \frac{\partial X(c_k)}{\partial c_k}\\ 0& -X(c_k)\end{array}\right)\left(\begin{array}{c}B(\beta )\\ C(\beta )\end{array}\right).$$

$$G=\left(\begin{array}{cc}-X(c_k)& \frac{\partial X(c_k)}{\partial c_k}\\ 0& -X(c_k)\end{array}\right),$$

$$\left(\begin{array}{c}B(1)\\ C(1)\end{array}\right)=e^G\left(\begin{array}{c}B(0)\\ C(0)\end{array}\right),$$

Full text

I Introduction

The principle of locality, arising naturally in physical systems, states that objects are only affected by their nearby surroundings. Locality is naturally embedded in numerous physical systems characterized by local Hamiltonians, which play a critical role in various quantum physics topics, such as quantum lattice models Tarasov et al. (2016); Tarasov (1985); Zanardi (2002), quantum simulation Feynman (1982); Cirac and Zoller (2012); Lloyd (1996); Buluta and Nori (2009), topological quantum computation Freedman et al. (2003), adiabatic quantum computation Nagaj (2008); Aharonov et al. (2008); Jordan et al. (2006), and quantum Hamiltonian complexity  Kempe et al. (2006); Kempe and Regev (2003); Bravyi et al. (2006). In the past few years, rapidly developing machine learning techniques allow us to study these topics in a new manner Xin et al. (2018); Bairey et al. (2019); Deng et al. (2017); Lu et al. (2019); Pudenz and Lidar (2013); Liu et al. (2017). Empowered by traditional optimization methods and contemporary machine learning techniques, we map the task of revealing the information encoded in a single eigenstate of a local Hamiltonian to an optimization problem.

For a local Hamiltonian $H=\sum_{i}c_{i}A_{i}$ with $A_{i}$ being some local operators, it is known that a single (non-degenerate) eigenstate $|\psi\rangle$ can encode the full Hamiltonian $H$ in certain cases Garrison and Grover (2018); Qi and Ranard (2019); Chen et al. (2012a), such as when the expectation value of $A_{i}$s on $|\psi\rangle$, a_{i}=\mbox{\langle\psi|}A_{i}|\psi\rangle, are given and further assumptions are satisfied. A simple case is that when $|\psi\rangle$ is the unique ground state of $H$; thus the corresponding density matrix of $|\psi\rangle$ can be represented in the thermal form as

[TABLE]

for sufficiently large $\beta$. This implies that the Hamiltonian $H$ can be directly related to $|\psi\rangle$, and hence can be determined by the measurement results $a_{i}$, using algorithms developed in the literature  Zhou (2008); Niekamp et al. (2013) for practical cases. Because $A_{i}$s are local operators, the number of parameters of $H$ (i.e., the number of $c_{i}$s) is only polynomial in terms of system size. We remark that the problem of finding $H$ is also closely related to the problem of determining quantum states from local measurements Xin et al. (2018); Bairey et al. (2019); Kalev et al. (2015); Linden et al. (2002); Linden and Wootters (2002); Diósi (2004); Chen et al. (2012b, c, 2013); Zeng et al. (2015), and also has a natural connection to the study of quantum marginal problem Klyachko (2006); Liu (2006); Liu et al. (2007); Wei et al. (2010), as well as its bosonic/fermionic version that are called the $N$-representability problem Coleman (1963); Erdahl (1972); Klyachko (2004, 2006); Altunbulak and Klyachko (2008); Schilling et al. (2013); Walter et al. (2013); Sawicki et al. (2014).

For a wavefunction $|\psi\rangle$ that is an eigenstate (i.e. not necessarily a ground state), one interesting situation is related to the eigenstate thermalization hypothesis (ETH) Deutsch (1991); Srednicki (1994, 1996, 1999); D’Alessio et al. (2016). When the ETH is satisfied, the reduced density matrix of a pure and finite energy density eigenstate for a small subsystem becomes equal to a thermal reduced density matrix asymptotically Garrison and Grover (2018). In other words, Eq. (1) will hold for some eigenstate $|\psi\rangle$ of the system in this case, and one can use a similar algorithm Zhou (2008); Niekamp et al. (2013) to find $H$ from $a_{i}$s, as in the case of ground states. Another situation previously discussed is that if the two-point correlation functions \mbox{\langle\psi|}A_{i}A_{j}|\psi\rangle are known, one can reproduce $H$ without satisfying ETH Chen et al. (2012a); Qi and Ranard (2019). Once the two-point correlation functions are known, one can again use an algorithm to recover $H$ from the correlation functions, for the case of ground states. However, in practice, the nonlocal correlation functions \mbox{\langle\psi|}A_{i}A_{j}|\psi\rangle are not easy to obtain Li et al. (2017).

In this paper, we answer a simple but significant question: can we determine a local Hamiltonian ($c_{i}$s) from only the local information ($a_{i}$s) of any one of the eigenstates ($|\psi\rangle$) without further assumptions. Obviously, there are cases that an eigenstate can not determine the system Hamiltonian, such as the product state and the eigenstates of a frustration-free Hamiltonians. However, for a randomly chosen physical system of which the Hamiltonian has no special structures, we show that only the knowledge of $a_{i}$s is sufficient to determine $H$.

Based on the knowledge in hand - a set of possible Hermitian operators $\{A_{i}\}$ for the system Hamiltonian and the measurement outcomes $\{a_{i}\}$, we formulate a positive-semidefinite function $f(\vec{x})$ (Eq. 6), with the $f(\vec{c})=0$ for the desired $\vec{c}$. The problem of finding the exact Hamiltonian converts to a unconstrained optimization problem of finding $\vec{x}$ that minimizing $f$. A small real number $\epsilon$ is chosen to control the precision of the result. We update the sampled $\vec{x}$ until $f(\vec{x})$ is less than $\epsilon$. Fig. 1 depicts the procedure of our algorithm.

We test the algorithm in two scenarios: one with four randomly generated operator $A_{i}$s acting on the whole system, and one with two random local operators $A_{i}$s. Our algorithm almost perfectly reproduces $c_{i}$s - the average fidelities for both cases are close to $1$. Since the algorithm recovers Hamiltonians with almost perfect fidelities Based on the reconstructed Hamiltonian $H$ (i.e. $c_{i}$s), one can also recover the eigenvalue of $|\psi\rangle$ from $c_{i}$s and $a_{i}$s, and the wave function itself from the eigenvectors of Hamiltonian $H$. In the case when $A_{i}$s are local operators, our method can find $c_{i}$s from only local measurement results, hence shed light on the correlation structures of eigenstates of local Hamiltonians.

II Algorithm

We start to discuss our method in a general situation, where the Hamiltonian $H$ can be expressed in terms of a set of known Hermitian operators $\{A_{1},A_{2},\ldots,A_{m}\}$: $H=\sum c_{i}A_{i}$ with $(c_{1},c_{2},\ldots,c_{m})=\vec{c}$. For an eigenstate $|\psi\rangle$ of $H$ with unknown eigenvalue $\lambda$, which satisfies $H|\psi\rangle=\lambda|\psi\rangle$, we can denote the measurement results as $a_{i}=\langle\psi|A_{i}|\psi\rangle$. With only knowing the measurement results $a_{i}$s, our goal is to find the coefficients $\vec{c}$ to determine $H$.

We observe that, even if $|\psi\rangle$ is not the ground state of $H$, it can be the ground state of another Hamiltonian $\tilde{H}^{2}$ with $\tilde{H}$ given by

[TABLE]

since

[TABLE]

Then the density matrix of $|\psi\rangle$ (which is in fact of rank 1) can be written in the form of a thermal state:

[TABLE]

for sufficiently large $\beta$, satisfying

[TABLE]

With these conditions in mind, we are ready to reformulate our task as an optimization problem with the following objective function:

[TABLE]

where $\vec{x}$ is the estimation of $\vec{c}$, and $\vec{x}=\vec{c}$ when $f(\vec{x})$ is minimized.

Notice that the first term of $f(\vec{x})$ is minimized by $\textrm{tr}(A_{i}\rho(\vec{c}))=a_{i}$, which guarantees that the state $\rho(\vec{c})$ obtained is the state that produces the desired measurement outcomes on $A_{i}$s. However, there may be many (thermal) states which also yield such outcomes. By simply minimizing this first term, optimization algorithms tend to return thermal states $\rho(\vec{c})$ with nonzero entropy, which is not the eigenstate of $H$ (with entropy zero) that we are willing to find. In order to fix this issue and ensure that the optimization returning a (near) rank $1$ state $\rho(\vec{c})$, we add the second term, which is only zero unless $\rho(\vec{c})$ is the ground state of $\tilde{H}^{2}$, hence an eigenstate of $H$. Combining these two terms together, we make sure that when the minimum value of Hamiltonian $f(\vec{x})$ is reached, we will obtain a $\rho(\vec{c})$ corresponding to measurement outcomes $a_{i}$, and at the same time an eigenstate of $H$. In practice, we set up a parameter $\epsilon$, such that with $f(\vec{x})<\epsilon$, we find a result with high fidelity to the desired value of $\{c_{i}\}$.

For the convenience of numerical implementation, we let $\tilde{H}_{\beta}=\sqrt{\beta}\tilde{H}$, then the thermal state $\rho(\vec{c})$ becomes

[TABLE]

Consequently, the objective function Eq. 6 can be rewritten in an equivalent form

[TABLE]

We aim to solve for $f(\vec{x})=0$ by minimizing $f(\vec{x})$. In practice, we terminate our iterations when $f(\vec{c})$ is smaller than a fixed small value $\epsilon$. The corresponding optimization result is denoted as $\vec{c}_{\text{opt}}$. As we reformed the objective function by using $\tilde{H}_{\beta}$, the result is actually $\vec{c}_{\text{opt}}=\sqrt{\beta}\cdot\vec{c}$.

Theoretically, we need that $\beta$ turns to infinity to pick up the ground state of $\tilde{H}_{\beta}^{2}$. In practice, however, we only require that $\beta$ is some “large number". What is more, to pick up the ground state of $\tilde{H}_{\beta}^{2}$ for $\rho(\vec{c})$, what really matters is in fact the gap between the first excited state and the ground state of $\tilde{H}_{\beta}^{2}$. Therefore, we just simply set $\beta=1$ and let the optimizer automatically amplify the energy gap during iteration, when $\rho(\vec{c})$ is approaching the desired state. A more detailed discussion regarding the choice of $\beta$ can be found in  Appendix A.

Since all the constraints are written in Eq. 8, minimizing $f(\vec{c})$ is an unconstrained minimization problem. There are plenty of standard algorithms for this task. In our setting, computing the second-order derivative information of $f(\vec{c})$ is quite complicated and expensive. Therefore, instead of using Newton method which requires the Hessian matrix of the second derivatives of the objective function, we choose the quasi-Newton method with BFGS-formula to approximate the Hessians BROYDEN (1970); Fletcher (1970); Goldfarb (1970); Shanno (1970). The MATLAB function fminunc, which uses the quasi-Newton algorithm, is capable of realizing this algorithm starting from an initial random guess of $c_{i}$s. When the quasi-Newton algorithm fails (converges into a local minimum), we start with a different set of random initial value $c_{i}$s and optimize again until we obtain a good enough solution.

The BFGS algorithm is a typical optimization algorithm which requires gradients of the objective function on $c_{i}$s. The form of the objective function $f(\vec{c})$ is so complicated such that computing the gradient is a difficult task. To solve this issue, we borrow the methods of computational graph and auto differentiation from machine learning. The computational graph is shown in Fig. 10, in which we show the intermediate functions and variables. Mathematically, the final gradients could be calculated via chain rules. However, since some of the intermediate variables are matrices and complex-valued, the automatic differentiation toolboxes, which deal with real variables, can not be applied directly. To obtain the gradients, we have to careful handle the derivation of the intermediate functions, especially those with matrices as their variables. More details about the our gradient method can be found in Appendix B.

We use Matlab to implement our algorithm, and the detailed implementation can be found in Appendix C.

III Results

In this section, we test our algorithm in three steps as follows. First, we randomly generate several $A_{i}$s and $c_{i}$s, hence the Hamiltonian $H_{rd}$ and its eigenstates. Second, for each Hamiltonian, we randomly choose one eigenstate $|\psi\rangle$, therefore we have $a_{i}=\langle\psi|A_{i}|\psi\rangle$ and $\rho=|\psi\rangle\langle\psi|$. Hereby we can run our algorithm to find the Hamiltonian $H_{al}$. Comparing $H_{al}$ and $H_{rd}$, we then know that how well the approach works. The algorithm has been tested for two cases: $A_{i}$ being the generic operator and local operator.

To compare $H_{al}$ and $H_{rd}$, we need a measure to characterize the similarity, or distance between these two Hamiltonians. The metric we used here is the following fidelity as discussed in Fortunato et al. (2002):

[TABLE]

To see the meaning of this metric, notice that $\textrm{tr}H_{al}H_{rd}$ is the inner product of the two Hamiltonians, while $\sqrt{\operatorname{Tr}H_{al}^{2}}$ and $\sqrt{\operatorname{Tr}H_{rd}^{2}}$ are the two normalization constant. Therefore, $f(H_{al},H_{rd})\in[0,1]$. If $H_{al}$ and $H_{rd}$ describe the same system up to a constant $b\in\mathbb{R}$, then $f(H_{al},H_{rd})=1$. Smaller $f(H_{al},H_{rd})$ then indicates that the two Hamiltonians are more far apart. Moreover, notice that in our settings the Hamiltonians are represented by vectors $\vec{c}$ and $\vec{c}^{\prime}$ in $m$-dimensional real space. If the chosen $A_{i}$s are normalized and orthogonal, which means

[TABLE]

where $d$ is the normalization constant given by the system dimension (e.g. for Pauli matrices, $d=2$), this fidelity definition is exactly the cosine loss function of $\vec{c}$ and $\vec{c}^{\prime}$, where $\vec{c}^{\prime}$ is generated from our algorithm, that is, for normalized $A_{i}$s,

[TABLE]

where $\parallel\vec{c}\parallel$ means the 2-norm of the vector $\vec{c}$.

III.1 Results with general operators

When implementing our method on generic operators, we randomly generate three Hermitian operators $A_{i}$s and fix them. We also create several real constants $c_{i}$s randomly, then assemble $A_{i}$s and $c_{i}$s into Hamiltonians

[TABLE]

After diagonalizing $H$, we choose one eigenvector $|\psi\rangle$ and calculate the expectation values $a_{i}=\langle\psi|A_{i}|\psi\rangle$.

For $n$-qubit systems, the dimension of the system is $d=2^{n}$. We tested the cases for $n=4,5,6$ and $7$. For each $n$, we generate $200$ data points. Each test set is $\{c_{i},A_{i},a_{i}|i=1,2,3\}$, where $c_{i}\in(0,1)$ for $i=1,2,3$ and $A_{i}$s are randomly generated Hermitian matrices of dimension $d=2^{n}$ 111To generate a random Hermitian matrix, we first generate $d^{2}$ complex number $\{e_{ij}|0\leq i,j\leq d\}$ as the entries of a matrix $E$, then we construct $A=E+E^{\dagger}$. It can be easily seen $A$ is an Hermitian matrix.. With these data, $H_{rd}=\sum_{i}c_{i}A_{i}$ is obtained. Diagonalizing $H_{rd}$ and randomly choosing one eigenvector $|\psi\rangle$, we obtain $\{a_{1},a_{2},a_{3}\}$. We show the results for applying our algorithm to all cases in Fig. 2.

We find that the final fidelities are larger than 99.8% for all tested cases. Although the fidelity slightly decreases with the increasing number of qubits, the lowest fidelity (7-qubit case) is still higher than 99.8%. Because the eigenvector $|\psi\rangle$ is randomly chosen, our method is not dependent on the energy level of eigenstates.

III.2 Results with local operators

In this section, we report our results on the systems with a $2$-local interaction structure. We tested two different structures, shown in Fig. 3(a) and Fig. 3(d). Each circle represents a qubit on a lattice, and each line represents an interaction between the connected two qubits.

The fully-connected $4$-qubit system is shown in Fig. 3(a). The Hamiltonian can be written as

[TABLE]

where $c_{ij}$s are real parameters and $A_{ij}$s are random generated $2$-local operators. One eigenstate out of 16 is randomly chosen as the state $|\psi\rangle$. Our algorithm has been tested on 800 such Hamiltonians.

We then analyze the $7$-qubit chain model shown in Fig. 3(d). Similarly, we can write the Hamiltonian as

[TABLE]

where $c_{i,i+1}$s and $A_{i,i+1}$s are the parameters and 2-local interactions. We randomly generate $40$ such Hamiltonians and applied our algorithm.

The results of these two $2$-local Hamiltonians are shown in Fig. 3. Our algorithm recovered Hamiltonians with high fidelities for both cases. The average fidelity for our 4-qubit (7-qubit) system is 99.99% (99.73%). As the dimension of the system increases, the fidelity between $H_{rd}$ and $H_{al}$ is slightly decreased. The histogram of the fidelities shows that, for most data points, the fidelities are very close to $1$. Our algorithm almost perfectly recovered these $2$-local Hamiltonian from measurement outcomes of a randomly picked eigenvector.

We examine the effectiveness of our algorithm according to different eigenvectors of the same Hamiltonians. Fig. 4 demonstrates average fidelities between $H_{rd}$ and $H_{al}$ of $4$-qubit case by chosen different eigenvectors. These average fidelities are higher than 99.9% for all $16$ eigenvectors. Hence, the effectiveness of our algorithm is independent of energy levels.

IV Further analysis of the algorithm

In this section, we analyze the error tolerance and the performance of the algorithm, based on the results of our numerical experiments.

IV.1 Error tolerance analysis

The numerical tests in previous sections deal with noiseless theoretical data. In practical scenarios, however, data is always noisy. Here we provide analysis of error tolerance for our algorithm.

As an example, we consider a 4-qubit system with Hamiltonian $H=\sum_{i=1}^{3}c_{i}A_{i}$ where $A_{i}$’s are random generated $16$ by $16$ Hermitian operators. Choose one eigenstate $|\psi\rangle$ of $H$, the noiseless measurements of the eigenstate are denoted as \{a_{i}|a_{i}=\mbox{\langle\psi|}A_{i}|\psi\rangle\}. Noises used here are randomly drawn from normal distributions $\gamma\sim\mathcal{N}(0,\sigma^{2})$. Adding the generated noise $\gamma$ to measurements $a_{i}$, the noisy data $a_{i,m}$ follows the normal distribution $\mathcal{N}(a_{i},\sigma^{2})$.

We can control how noisy $a_{i,m}$ is by changing the standard deviation $\sigma$. Clearly the noisiness of data is relative to the magnitude of true values, thus $\sigma/a_{i}$ can be used as an indicator. We test ten different $\sigma$’s of which $\sigma/a_{i}$’s ranging from $0.01$ (1%) to $0.1$ (10%). For each $\sigma$, $1000$ data points have been generated.

The Hamiltonians attained by our algorithm with these noisy data $a_{i,m}$ are denoted as $H_{\text{err}}$. The fidelity between each pair of the true Hamiltonian $H$ and $H_{\text{err}}$ is shown in Fig. 5. Though increasing the noise cause the fidelities of a few points significantly decreased, all median fidelities of different $\sigma$ are above 99.8%. Even for $\sigma/a_{i}=0.1$ (added 10% of noise), only 2.1% of data points have fidelities lower than 99.0%. In other words, singular points which have significant low fidelities are rare. This is demonstrates that the performance of our algorithm is stable.

IV.2 Performance analysis

The convergence and time cost of our method are analyzed in this subsection. We study the relation between the halting condition $\epsilon$ and the convergence of our algorithm. We also observed that the time cost tends to differ for each Hamiltonian configuration.

The parameter $\epsilon$, which is the halting condition of the algorithm, effects the convergence of our algorithm. It is the consequence of the object function $f(\vec{x})$ has many local minima. In the schematic diagram Fig. 7, if the initial point is chosen as $A$ or $A^{\prime}$, the gradient method will return the value $A_{0}$, which is a local minimum (so as $B$ or $B^{\prime}$ to $B_{0}$, and $C$ or $C^{\prime}$ to $C_{0}$).

When $\epsilon$ is greater then $\|f(A_{0})-f(B_{0})\|$ or $\|f(C_{0})-f(B_{0})\|$, the algorithm may recognize $A_{0}$ or $C_{0}$ as the optimal solution instead of finding the true global minimum $B_{0}$. An appropriate choice of $\epsilon$ is necessary to eliminate certain local minima.

The object function Eq. 6 used in our method has high dimension and complicated landscape. Its properties also subject to the particular class of $H$ that we work with. More analysis could be done in terms of finding an appropriate $\epsilon$. Empirically, we choose $\epsilon$ as $10^{-6}$ for all our numerical experiments. It is numerically proved to be applicable as shown in Fig. 7. The figure depicts the relation between $\epsilon$ and the output fidelity of the 3-qubit general operators case. Chosen $\epsilon\leq 10^{-3}$ grantees the average fidelity almost equals to 1.

By setting the halting parameter $\epsilon$ to $10^{-6}$, we can now discuss the time cost for each numerical experiment. The total time cost $t$ for optimization depends on the number of trials $n$ for each numerical example and the time $t_{0}$ of each trial. The initial values are chosen randomly, $n$ is different from case to case. We consider the average number of trials $\bar{n}$ for each Hamiltonian class. Let $t_{0}$ be the duration of a single trial which is mainly depending on the Hamiltonian configuration and number of qubits. The total average time $\bar{t}$ could be estimated as

[TABLE]

We test our algorithm on a workstation with Intel i7-8700K and 32 GB RAM, the results are listed in Table 1. The table demonstrates that the time cost does not grow rapidly for the general operator cases, while it changes dramatically with the system size for the local operator cases. The $\bar{t}$ of 7-qubit general operator instance is almost 2 times more then $\bar{t}$ of the 4-qubit general operator case. On the contrast, the average time cost $\bar{t}$ of 7-qubit lattice is almost 350 times more then the $\bar{t}$ of 4-qubit lattice.

V Comparison with the correlation matrix method

In Ref. Qi and Ranard (2019), a method is proposed to recover Hamiltonians from correlation functions. Their method works as follows: with an eigenstate $|\psi\rangle$, the Hamiltonian of the system is defined as

[TABLE]

where $L_{i}$s are normalized and orthogonal Hermitian matrices (e.g. the Pauli matrices). They defined a correlation matrix, of which the matrix elements are

[TABLE]

where $\{L_{i},L_{j}\}=L_{i}L_{j}+L_{j}L_{i}$ is the anti-commutator. Then diagonalize matrix $M$ and find the eigenvector $\vec{\omega}=(\omega_{1},\omega_{2},\cdots,\omega_{i})^{T}$ corresponding to the eigenvalue [math]. The Hamiltonian

[TABLE]

is the desired Hamiltonian of the system $H$. We refer to this algorithm as the correlation matrix method (CMM).

We conduct numerical experiments to compare the performance of CMM and our method. We calculate such 4-qubit Hamiltonians

[TABLE]

where $\sigma_{k}^{i}$ is the $k-$th Pauli matrix acting on the $i-$th qubit. It turns out that CMM and our algorithm both render good estimations. The results of CMM possess greater accuracy–the error rate, defined as $r(H_{\text{rd}},H_{\text{al}})=1-f(H_{\text{rd}},H_{\text{al}})$, is less than $10^{-10}$. Error rates of our results range from $10^{-4}$ to $10^{-2}$.

CMM is deterministic, more accurate as well as faster than our method. It only involves diagonalization of matrix $M$, of which the dimension is polynomial of the number of qubits for local systems. CMM also provides a criteria to determine whether the eigenstate is uniquely determine the Hamiltonian. That is, if $M$ only has one eigenvalue equals to [math], then the Hamiltonian is uniquely determined. The advantage of CMM arises from more restrictions and more information required for applying it: (1) All $L_{i}$s in CMM are orthogonal, while in our tests the corresponding $A_{i}$s are randomly generated; (2) CMM requires non-local correlation functions. The matrix element includes the term

[TABLE]

where $\{L_{i},L_{j}\}=L_{i}L_{j}+L_{j}L_{i}$. Although $L_{i}$s may be local, $L_{i}L_{j}$ for all $i$s and $j$s could be global. Therefore, as indicated in Qi and Ranard (2019), if only partial knowledge of the system is available, the question that whether the Hamiltonian can be reconstructed remains unclear.

Another problem of applying CMM is similar to the problem of the halting parameter $\epsilon$ in our method. The recovery depends on the existence of one eigenvalue equals to [math]. However, the equivalence of tow numbers in a computer is different from in theory. Even when working with noiseless theoretical data, the finite length data storage (e.g. the precision of double float-point format is $10^{-16}$) and data processing can introduce certain errors, not to mention the noisy data from real quantum devices. Namely, determine whether a number is equal to [math] is up to a certain precision. Therefore, one needs to set a threshold $\delta$ to determine equivalence before applying CMM. How to appropriately chose a $\delta$ is empirical and case dependent. From this perspective, both methods are not similarly complicated.

We are also willing to check the consistency of CMM and our methods. Lacking a systemic way of choosing the threshold $\delta$, we need to seek another way to bridge these two method. From the numerical experiments, we find in CMM that the lowest absolute eigenvalue (denotes as $\lambda_{0}$) of $M$ is about $10^{-16}\leavevmode\nobreak\ 10^{-18}$ and the second lowest absolute eigenvalue (denotes as $\lambda_{1}$) ranges from $10^{-2}$ to $10^{-12}$. If $\lambda_{1}$ is too small, it would be hard to tell whether there exists one or more 0 eigenvalues, that is, the farther $\lambda_{0}$ and $\lambda_{1}$ are, the more accurate the CMM’s outcome is. It will be an evidence of consistency if we can witness the same tendency from our method. We, therefore, define the condition number as the ratio of the second lowest absolute eigenvalue to the lowest absolute eigenvalue ($\lambda_{1}/\lambda_{0}$). We use the estimations provided by our method to construct the correlation matrix $M$, obtain the $\lambda_{0}$ and $\lambda_{1}$ from $M$, then calculate the condition number $\lambda_{1}/\lambda_{0}$. The results is shown in Fig. 8. We find, generally, the larger the condition number, the better our algorithm performs (higher fidelity). It is consistent with the behaviour of CMM.

In summary, for all the tests performed, CMM and our method has similar behaviour in terms of the condition number, which shows that both method return the desired results (for reconstruction the desired Hamiltonian with high fidelity). Regarding to algorithm performance, when all sorts of correlations are available, the CMM is numerically more efficient and accurate. In contrast, however, our algorithm use less information, therefore can be used in a much wider situations, especially when certain global correlation information is hard to obtain in practice.

VI Discussion

In this work, we discuss the problem of reconstructing system Hamiltonian $H=\sum_{i}c_{i}A_{i}$ using only measurement data a_{i}=\mbox{\langle\psi|}A_{i}|\psi\rangle of one eigenstate $|\psi\rangle$ of $H$ by reformulating the task as an unconstrained optimization problem of some target function of $c_{i}$s. We numerically tested our method for both randomly generated $A_{i}$s and also the case that $A_{i}$s are random local operators. In both cases, we obtain good fidelities of the reconstructed $H_{al}$. Our results are somewhat surprising: only local measurements on one eigenstate are enough to determine the system Hamiltonian for generic cases, no further assumption is needed. Though theoretically it is beautiful, our algorithm is not scalable since the calculations of exponentiation and gradient of matrices become expensive when the system size gets large. Improvements can be done by more efficiently approximate them.

We also remark that, in the sense that our method almost perfectly recovered the Hamiltonian $H$, the information encoded in the Hamiltonian $H$, such as the eigenstate $|\psi\rangle$ itself (though described exponentially many parameters in system size) can also in principle be revealed. This builds a bridge from our study to quantum tomography and other related topics of local Hamiltonians. Empowered by traditional optimization methods and machine learning techniques, our algorithm could be applied in various quantum physics problems, such as quantum simulation, quantum lattice models, adiabatic quantum computation, etc.

Our discussion also raises many new questions. For instance, a straightforward one is whether other methods can be used for the reconstruction problem, and their efficiency and stability compared to the method we have used in this work. One may also wonder how the information of “being an eigenstate" helps to determine a quantum state locally, which is generically only the case for the ground state, and how this information could be related to help quantum state tomography in a more general setting.

Aknowlegement

We thank Feihao Zhang, Lei Wang, Shilin Huang, and Dawei Lu for helpful discussions. S.-Y.H is supported by National Natural Science Foundation of China under Grant No. 11847154. N.C. and B.Z. are supported by NSERC and CIFAR.

Appendix A The value of $\beta$

It is easy to observe that $\tilde{H}^{2}$ and $\tilde{H}_{\beta}^{2}$ have the same eigenstates. The constant $\beta$ only contributes a constant factor to the magnitude of eigenvalues, i.e., the eigenenergies of the given system. Theoretically, a thermal state tends to the ground state of the given Hamiltonian if and only if the temperature is zero (or, for a numerical algorithm, close to zero), which means $\beta=\frac{1}{kT}$ goes to infinite. Numerically, the $\beta$ only needs to be a sufficiently large positive number.

Let us denote the $i$-th eigenvalue of $\tilde{H}_{\beta}^{2}$ as $E_{i}$ and the energy gaps as $\Delta_{i}=E_{i+1}-E_{i}$. From the definition of $\tilde{H}_{\beta}^{2}$, the ground energy $E_{g}=E_{1}$ is always 0. As we observed, during the optimization process, the eigenenergies, as well as the energy gaps of the Hamiltonian, gets larger. From Fig. 9(a), we can see the energy gaps grow as the optimization goes. The corresponding $\beta$ is a finite sufficient large number.

On the other hand, we can also examine the probability of the ground state of $\tilde{H}^{2}_{\beta}$ as the iteration goes. This probability can be expressed as

[TABLE]

Since the ground state energy of $\tilde{H}_{\beta}^{2}$ is always have zero, the probability is

[TABLE]

The change of the probability $P(E=E_{g})$ with iterations is shown in Fig. 9(b), from which we can see that finally, the probability goes to $1$, which means that the thermal state form is (almost) ground state when $\beta$ is a relatively large positive constant.

Appendix B Calculation of the gradient of $f(\vec{c})$

B.1 General method

The computational graph of this technique is shown in Fig. 10. Each node of the graph represents an intermediate function. In principle, the chain rule can be applied. It seems that we can compute the gradient by using the existing deep learning frameworks. However, those auto differentiation tools embedded in the frameworks such as PyTorch Paszke et al. (2017) and TensorFlow Abadi et al. (2015) cannot harness our problem due to the following two reasons: first, most of the tools only deal with real numbers while in quantum physics we often deal with complex numbers; second and the most important, some of the intermediate functions in the computational graph use matrices (or in physics, operators) as variables, while in neural networks, the variables are real numbers. One should be careful while deriving matrices since a matrix usually does not commute with its derivative. For instance, even for some simple conditions, the node $v_{3}=\tilde{H}^{2}=v_{2}^{2}$, the derivative $\frac{\partial v_{3}}{\partial v_{1,k}}$ can not be simply considered as $2v_{2}\frac{\partial v_{2}}{\partial v_{1,k}}$ because $[\tilde{H},\frac{\partial\tilde{H}}{c_{k}}]\neq 0$. Instead, $\frac{\partial v_{3}}{\partial v_{1,k}}=v_{2}\frac{\partial v_{2}}{\partial v_{1,k}}+\frac{\partial v_{2}}{\partial v_{1,k}}v_{2}$.

The derivatives in the computational graph shown in Fig. 10 are listed in Table 2. Here, due to the reason mentionedabove, we use forward propagation to calculate the gradients. In Table 2, the derivative could be separated into the following categories: 1. The variable(s) of the function is a number. This corresponds to the simplest case, and the derivatives could be obtained simply using chain rules; 2. The variables are matrices, but the function is the trace function tr. The derivatives could be obtained by applying $\frac{d\textrm{tr}(A)}{dx}=\textrm{tr}\frac{dA}{dx}$; 3. Some elementary functions take matrices as variables, i.e., $v_{3}=v_{2}^{2}$. In this case, one should be careful about the commutators when dealing with it; 4. $v_{4}=e^{\tilde{H}^{2}}$, of which the gradients are difficult to calculate, which will be discussed in detail.

B.2 Derivative of matrix exponentials

The derivative of function $f(X)=\exp(X(c_{k}))$ with respect to $c_{k}$ can be written as

[TABLE]

where, in our cases, $X(c_{k})=-\tilde{H}(c_{k})^{2}=-v_{3}$. Generally, the commutator $[X,\frac{\partial X}{\partial c_{k}}]\neq 0$. Therefore, the integration Eq. 21 can not be simply calculated. In Ref. Niekamp et al. (2013), the authors approximate the integration using the value of the upper and lower limits. This approach does not work for our case. Thus we introduce a new way to calculate it. Let

[TABLE]

[TABLE]

and

[TABLE]

Note that $A(1)=B(1)\exp[X(c_{k})]$. It can be derived that

[TABLE]

Let

[TABLE]

The derivative equation (25) can be solved as

[TABLE]

where $B(0)$ is a matrix with all entries 0 and $C(0)$ is the identity matrix. With $B(1)$, we can obtain $A(1)$ as well as the integration Eq. 21. This completes the gradient calculation for our algorithm.

Appendix C Software implementations

According to the computational graph Fig. 10 and the Table 2, we can define a function which accepts $\vec{c}$, the operators $A$, and returns the value of $f(\vec{c})$ and the gradient $\nabla f(\vec{c})$:

⬇

function [f,g]=fun(c, A)

% f is the value of the objective function

% g is the gradient

The function fminunc in MATLAB/Octave can accept the function and the gradient and doing the optimization:

⬇

f_handle=@(c)fun(c,A)

The default algorithm is BFGS. To make the algorithm using gradient we provide:

⬇

options=optimoptions(’fminunc’,’SpecifyObjectiveGradient’,’true’);

Then, we can start the optimization

⬇ [c1, fval]=fminunc(f_handle,c0)

Here, c1 is the final points, fval is the final value of the objective function and c0 is the initial value.