State preparation and evolution in quantum computing: a perspective from Hamiltonian moments

One straightforward way is to apply Pauli reduction and commutativity to reduce the number of terms, and then do term-by-term measurements. The Pauli reduction simply utilizes the commutation and anti-commutation relations between Pauli matrices

where $\sigma_{i},\sigma_{j},\sigma_{k}$ are Pauli matrices, $\epsilon_{ijk}$ is the structure constant following the Levi-Civita symbol Tyldesley (1975), and $\delta_{ij}$ is the Kronecker delta. From (12) it is straightforward to see the product of Pauli strings is another Pauli string with an appropriate phase factor. Therefore, in comparison to the evaluation of $\langle\hat{\mathcal{H}}\rangle$, the evaluation of $\langle\hat{\mathcal{H}}^{n}\rangle$ does not increase the circuit depth. Furthermore, we notice that the measurement of a single Pauli string could be used to determine classes of contributions to Hamiltonian moments of arbitrary order. For example, since every Pauli string is unitary, the evaluation of $\langle P_{i}\rangle$ can simultaneously contribute all the $\langle H^{2k+1}\rangle$’s ($k\geq 0,k\in\mathbb{N}$) Kowalski and Peng (2020). Suchsland et al. systematically analyzed the actual number of Pauli strings corresponding to $\hat{\mathcal{H}}^{n}$ ($n=1,2,3$) for a collection of 11 molecules, including H_n ($n=2,4,6,8,10$), LiH, NaH, H₂O, NH₃, and N₂, that require up to $N_{q}=18$ qubits Suchsland et al. (2021). It was found that the actual numbers of Pauli strings in these systems after applying (12) exhibit at most $\mathcal{O}(N_{q}^{3})$, $\mathcal{O}(N_{q}^{5})$, and $\mathcal{O}(N_{q}^{7})$ scalings in comparison to the predicted $\mathcal{O}(N_{q}^{4})$, $\mathcal{O}(N_{q}^{8})$, and $\mathcal{O}(N_{q}^{12})$ scalings for $\hat{\mathcal{H}}^{n}$ ($n=1,2,3$) respectively.

On the top of the above simplification, a more significant term reduction can be achieved by utilizing the commutativity of the Pauli strings. To see this method works mathematically, consider operators $\hat{\mathcal{A}}$ and $\hat{\mathcal{B}}$ that commute, i.e. $[\hat{\mathcal{A}},\hat{\mathcal{B}}]=0$, then there exists a complete orthonormal eigenbasis ${|u_{i}\rangle}$ ($i=0,1,\cdots$) that simultaneously diagonalizes both operators, and we can write

with $\{\lambda_{i},i=0,1,\cdots\}$ and $\{\kappa_{i},i=0,1,\cdots\}$ being the eigenvalues of $\hat{\mathcal{A}}$ and $\hat{\mathcal{B}}$, respectively. Now the expectation values of $\hat{\mathcal{A}}$ and $\hat{\mathcal{B}}$ with respect to a trial state $|\phi\rangle$ can be expressed as

Therefore, if we measure each $|\langle u_{i}|\phi\rangle|^{2}$’s, we can deduce $\langle\hat{\mathcal{A}}\rangle$ and $\langle\hat{\mathcal{B}}\rangle$ based on (14).

The simplest commutativity among Pauli strings is qubitwise commutativity (QWC) Kandala et al. (2017); McClean et al. (2016), where two Pauli strings qubitwise commute if two single-qubit Pauli matrices at each index commute, i.e.

For example, $\{\sigma_{x}\otimes\sigma_{y},\sigma_{x}\otimes I,I\otimes\sigma_{y},I\otimes I\}$ is a QWC class, since any pair of Pauli strings in this class commute (a more pictorial QWC grouping procedure is shown in Fig. 2). For the observables grouped in one QWC class, one can easily find a shared eigenbasis that simultaneously diagonalizes all the observables. Take the previous example, if we want to measure the QWC class in the computational basis, utilizing the Clifford representation of single Pauli matrices $\sigma_{x}$ and $\sigma_{y}$

we will only need perform the following measurement

Now, given a large amount of Pauli strings that constitute $\hat{\mathcal{H}}^{n}$, how should we find an optimal grouping of these Pauli strings such that the total number of measurements is minimized? Unfortunately, it turns out that finding the optimal grouping is NP-hard Karp (2010); Gokhale et al. (2020), and can be converted to minimum clique cover or minimum graph coloring problem, which in practice can be approximately solved through heuristic approaches that scale quadratically with the number of the Pauli strings Yen, Verteletskyi, and Izmaylov (2020). In terms of term reduction, for small molecules the QWC grouping can bring roughly 75% saving for evaluating $\langle\hat{\mathcal{H}}\rangle$, and even bigger savings for for $\langle\hat{\mathcal{H}}^{n}\rangle$ ($n>1$). Early studies of applying QWC bases to evaluating $\langle\hat{\mathcal{H}}^{4}\rangle$ for the Heisenberg models represented by 10 to 40 qubits exhibits a sub-linear scaling of the number of measurements in the number of the qubits Vallury et al. (2020). Furthermore, for some small molecular systems, it has been shown that the number of QWC bases grouping the Pauli strings that constitute $\hat{\mathcal{H}}^{n}$ eventually reach a plateau, regardless of the power $n$ Claudino et al. (2021).

Essentially, the QWC grouping has its roots in mutually unbiased bases (MUB) Schwinger (1960); Klappenecker and Rötteler (2004) from quantum information theory associated with maximizing the information learned from a single measurement. In general, MUBs group the potentially $4^{N_{q}}-1$ $N_{q}$-qubit Pauli strings (excluding identity string) into commuting groups of maximal size, of which using QWCs would allow reducing this requirement to $3^{N_{q}}$ Altepeter, James, and Kwiat (2004). Further, it is known that there in principle exists an MUB grouping for $N_{q}$-qubit that can further reduce this requirement to even $2^{N_{q}}+1$ with maximum $2^{N_{q}}-1$ distinct Pauli strings in each commuting group, however, the entanglement in the MUBs Lawrence, Brukner, and Zeilinger (2002) makes MUB quantum state tomography challenging. In practice, beside QWC, other grouping techniques have alternatively been proposed, including general commutativity Gokhale et al. (2020); Yen, Verteletskyi, and Izmaylov (2020), unitary partitioning Izmaylov et al. (2020), and Fermionic basis rotation grouping Huggins et al. (2021). Generally speaking, the applications of these grouping rules for evaluating $\langle\hat{\mathcal{H}}\rangle$ have shown that, at a cost of introducing additional one-/multi-qubit unitary transformation before the measurement, the total number of terms can be significantly reduced from $\mathcal{O}(N_{q}^{4})$ to $\mathcal{O}(N_{q}^{2\sim 3})$. In particular, for simpler cases where the entire Hamiltonian or its power could even be transformed by a single unitary, the evaluation of the corresponding expectation values can be done in a single set of measurements.

It is worth mentioning that the above discussion is limited to the number of terms that can be efficiently reduced through the groupings that feature the commutativity of Pauli strings, while the total number of measurements required from the number of groups also critically depends on the covariance between Pauli strings, $\text{cov}\big{(}P_{i},P_{j}\big{)}$, and the desired precision $\epsilon$. Given that we can always write the Hamiltonian powers as linear combination of Pauli strings similar to Eq. (2), the total number of measurements for evaluating the moments can then be expressed as Gonthier et al. (2020); Rubin, Babbush, and J. (2018); Wecker, Hastings, and Troyer (2015)

with $G$ indexing the groups. Therefore, it is likely that the number of groups decreases at the cost of introducing larger covariances that could essentially increase the total number of measurements required to achieve a desired precision.

There are also many other advanced measurement schemes proposed recently. One example is to simultaneously obtain expectation values of multiple observables by randomly measuring and projecting the quantum state into classical shadows Huang, Kueng, and Preskill (2020, 2021); Aaronson (2018); Struchalin et al. (2021); Chen et al. (2021); Zhao, Rubin, and Miyake (2021); Acharya, Saha, and Sengupta (2021); Hadfield (2021); Hillmich et al. (2021); Zhang et al. (2021). The algorithm in principle enables the measurements of $M$ low-weight observables using only $\mathcal{O}(\log_{2}M)$ samples. The practical performance of the algorithm for model and molecular Hamiltonians on NISQ device, in terms of accuracy and efficiency, is still under intense study.

Another way to exactly implement $\hat{\mathcal{H}}^{n}$ is by Chebyshev polynomial expansion. We can exactly express a monomial $x^{n}$ on $\left[-1,1\right]$ as a finite sum of Chebyshev polynomials of the first kind, i.e.

where

and $T_{k}(x)$ is the $k^{th}$ degree Chebyshev polynomial of the first kind Subramanian, Brierley, and Jozsa (2019). This expansion is useful because a “quantum walk” can exactly produce the effect of Chebyshev polynomials in $\hat{\mathcal{H}}/d$, where $d$ is the sparsity of the matrix. A quantum walk is a process in which the Hilbert space is enlarged, and then a walk operator that works in the enlarged Hilbert space is repeatedly implemented. For details on quantum walks and their implementations, see e.g. Refs. 68; 69. With amplitude amplification to obtain a constant success probability of preparing a quantum state and a flag qubit to indicate if the preparation was successful, the query complexity of this method is $\mathcal{O}\left(nd^{n-1}\lambda_{0}^{-n}\right)$ Subramanian, Brierley, and Jozsa (2019). Using this method without amplitude amplification, the query complexity of this implementation is $\mathcal{O}\left(n\right)$. Here, $\lambda_{0}$ denotes the smallest (magnitude) eigenvalue of $H$, with the assumption that $\hat{\mathcal{H}}$ does not have 0 as an eigenvalue. However, this method has its own drawbacks. The implementation of the quantum walk can be expensive, requiring a number of ancillary qubits that grows linearly with the system size Subramanian, Brierley, and Jozsa (2019). The Chebyshev polynomial expansion’s potential application for computing the imaginary time evolution of an unnormalized quantum state, and in general $e^{-\beta\hat{\mathcal{H}}}$ for evaluating partition function has been discussed in Ref. 70, where they use the recursion relation for Chebyshev polynomials.

Recently, Seki and Yunoki proposed a quantum power method to evaluate $\hat{\mathcal{H}}^{n}|\phi\rangle$ Seki and Yunoki (2021), where the Hamiltonian power $\hat{\mathcal{H}}^{n}$ is first expressed as the $n$-th time derivative of $\hat{\mathcal{U}}(t)=e^{i\hat{\mathcal{H}}t}$ at $t=0$, then the time derivative is approximated through the central finite difference (CFD) which can be alternatively expressed as a linear combination of $\hat{\mathcal{U}}(\Delta_{t})$’s at different time variables $\Delta_{t}$’s in close proximity to $t=0$, i.e.

To implement the approach to quantum circuit, each $\hat{\mathcal{U}}(\frac{\Delta_{t}}{2})$ is further decomposed using symmetric Suzuki-Trotter decomposition (SSTD) Trotter (1959); Suzuki (1990). There are two sources of error in this approach, the CFD error and SSTD error. The systematic CFD error scales quadratically with $\Delta_{t}$, i.e. $\mathcal{O}(\Delta_{t}^{2})$, while the $2m$-th order SSTD brings the systematic error of $\mathcal{O}(\Delta_{t}^{2m})$. However, there are several advantages in the proposed approach. First, both aforementioned errors can be systematically suppressed through Richardson extrapolation Richardson and Gaunt (1927); Temme, Bravyi, and Gambetta (2017) at the cost of involving more terms in the linear combination in Eq. (20), which on the other hand offers the space for using low order SSTD. Second, the number of gates required for approximate $\hat{\mathcal{H}}^{n}$ for an $N_{q}$-qubit Hamiltonian $\hat{\mathcal{H}}$ scales as $\mathcal{O}(nN_{q})$. In principle, the quantum power method can be combined with many classical approaches, which will be elaborated in Section III for targeting the ground/excited states and properties of a quantum many-body system. As shown in Ref. 71, numerical noiseless simulations employing the quantum power method in Krylov-subspace diagonalization for targeting ground state of model systems demonstrate systematically improved accuracy over the conventional VQE.

Inspired by the classical inverse power iteration approach for finding the dominant eigenstate of a given hermitian matrix with a more favorable logarithmic complexity in the iteration depth Sachdeva and Vishnoi (2014), Kyriienko recently proposed a quantum inverse iteration algorithm Kyriienko (2020) for approximately computing the inverse of $n$-th power of the Hamiltonian, $\hat{\mathcal{H}}^{-n}$. The approach extends the discretized Fourier approximation of the Hamiltonian inverse Childs, Kothari, and Somma (2017) to the $n$-th power,

where $N_{n}$ is a normalization factor, $y_{j}\equiv j\Delta_{y}$ and $z_{k}\equiv k\Delta_{z}$ with the integers $j\in[0,J]$ and $k\in[-K,K]$ and intervals $\Delta_{y,z}\in\mathbb{R}$ defining a discretization grid space, and unitary $\hat{\mathcal{U}}_{j,k}$ is given by $e^{-iy_{j}z_{k}\hat{\mathcal{H}}}$. Then $\langle\phi|\hat{\mathcal{H}}^{-n}|\phi\rangle$ is evaluated through the SWAP test Ekert et al. (2002); Higgott, Wang, and Brierley (2019), overlap measurement Mitarai and Fujii (2019), or Bell-type-like measurement Kyriienko (2020). However, this method requires a number of calls to a time-evolution oracle that is dependent on the condition number of the system, as well as a requirement that the Hamiltonian has been shifted such that all of the eigenvalues are positive. This latter requirement requires the user to have some knowledge of the ground state energy of the system a priori. A suitable guess can be found, for example, using the minimum label finding algorithm of Ref. 81.

Alternatively, if considering a fault-tolerant implementation with favorable resource scaling $\hat{\mathcal{H}}^{n}|\phi\rangle$ or $\hat{\mathcal{H}}^{-n}|\phi\rangle$ can be done through amplitude amplification approach Brassard et al. (2002), Hamiltonian simulation Berry et al. (2015), qubitization Low and Chuang (2017), or the direct block-encoding Gilyén et al. (2019) methods at the cost of introducing deeper circuit and implementing controlled-$\hat{\mathcal{U}}$ operations which, however, are still challenging to be implemented in the NISQ devices. Some of these methods have deep connection to linear combination of unitaries (LCU) technique Childs and Wiebe (2012). The basic idea of LCU can be briefly demonstrated from the following toy example. Suppose we want to implement $\hat{\mathcal{U}}_{0}+\hat{\mathcal{U}}_{1}$ on $|\phi\rangle$. We can start by introducing an ancilla qubit and preparing $|+\rangle|\phi\rangle$ with $|+\rangle\equiv\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$. Then by performing the controlled unitary $|0\rangle\langle 0|\otimes\hat{\mathcal{U}}_{0}+|1\rangle\langle 1|\hat{\mathcal{U}}_{1}$ on $|+\rangle|\phi\rangle$ we are able to obtain the state $\frac{1}{\sqrt{2}}(|0\rangle\hat{\mathcal{U}}_{0}|\phi\rangle+|1\rangle\hat{\mathcal{U}}_{1}|\phi\rangle)$. Now, if we measure the ancilla qubit in the $\{|+\rangle,|-\rangle\}$ basis and obtain the outcome associated with $|+\rangle$, then we obtain a state proportional to $(\hat{\mathcal{U}}_{0}+\hat{\mathcal{U}}_{1})|\phi\rangle$.

One of the most common approaches that can be exploited is the Lanczos approach Lanczos (1950). In the Lanczos approach, an orthonormal set of quantum states $\{|\phi_{i}\rangle,i=1,2,\cdots\}$ is generated through a three-term recursion

where $\alpha_{i}=\langle\phi_{i}|\hat{\mathcal{H}}|\phi_{i}\rangle$, $\beta_{i}=\langle\phi_{i}|\hat{\mathcal{H}}|\phi_{i+1}\rangle$. The three-term recursion essentially transforms the Hamiltonian $\hat{\mathcal{H}}$ into a tridiagonal form

where the matrix elements $\alpha_{i}$ and $\beta_{i}$ can essentially be represented recursively in terms of Hamiltonian moments. Explicitly, it can be shown

Now the strict upper bound of the ground state energy of the Hamiltonian can be obtained by directly diagonalizing $\hat{\mathcal{T}}_{i}$. It has been shown that the classical Lanczos scheme can directly be applied even employing the noisy $\langle\hat{\mathcal{H}}^{n}\rangle$ Suchsland et al. (2021), which effectively provides a way to improve the energy estimate from real quantum hardware without being subject to specific hardware and an explicit description of the underlying noise. Preliminary tests on H₂, H₃, and four-site tetrahedral Heisenberg model demonstrate that combining the Lanczos scheme with VQE algorithm helps correct the noise and enhance the quality of the ansatz.

Alternatively, the matrix elements $\alpha_{i}$ and $\beta_{i}$ can be expressed in terms of the so-called connected moments (or cumulant) Horn and Weinstein (1984). An $n$-th order connected moments is defined as

from which (28) can be generalized as $z$-expansions of $\alpha$ and $\beta$ Hollenberg (1993); Witte and Hollenberg (1994); Hollenberg and Witte (1996)

where $z\equiv i/V$ with $i$ the recursion index and $V$ the volume of the system. From the bound analysis of orthogonal polynomials van Doorn (1987); Ismail and Li (1992), it has been shown that the true ground state energy of the Hamiltonian can be approximated through the greatest lower bound (i.e. infimum) of the expansion (32) Hollenberg and Witte (1996),

For example, the first order in $z$ gives Hollenberg and Witte (1994)

Approach (34) has been employed to demonstrate the energy estimate of 2D quantum magnetism model on lattices up to 25 qubits on IBM-Q quantum hardware Vallury et al. (2020), where the results show a consistent improvement of the estimated energies in comparison to the VQE results for the same ansatz. Another application of the Lanczos method is the continued fractions expansion of the Green’s function. The Green’s function can be calculated in the Lehmann representation using the $\alpha$ and $\beta$ parameters from the Lanczos method. For full details, see e.g. Chapter 8 of Ref. 95.

It is worth mentioning that standard Lanczos is a variational approach while the infimum approach (33) does not offer a strict upper bound of the true ground state energy, which, as exemplified by (34), is essentially an alternative size-extensive connected moment expansion of the ground state energy. Other polynomial expansion approaches, in particular connected moment expansions, will be discussed in details in Section III.3.

Alternatively, one can resort to time evolution approaches to target the energy of the Hamiltonian. Here the time evolution approaches include both real time evolution (RTE) approach and its imaginary analogue. In RTE, the dynamics of a quantum state $|\phi_{t}\rangle$ is governed by the time-dependent Schrödinger equation

In the Schrödinger picture the solution of Eq. (35), $|\phi_{t}\rangle$, can be expressed by acting a time-evolution operator $\hat{U}(t)=e^{-i\hat{H}t}$ (or its truncated Taylor expansion in practice) on an initial state $|\phi_{0}\rangle$,

The introduction of $\hat{\mathcal{H}}^{n}$ in (36) effectively builds a non-orthogonal Krylov subspace. Now consider an order-($r+1$) Krylov subspace $\mathcal{K}_{r+1}=\{\hat{\mathcal{H}}^{n}|\phi_{0}\rangle,n=0,1,\cdots,r\}$, the ground state can then be approximated as a linear combination of the Krylov basis

under the constraint $\langle\phi_{t}|\phi_{t}\rangle=1$ with $c_{n}$’s being the expansion coefficients to be determined. Plugging the expanded form (37) to Eq. (35) and projecting the approximate evolution on to the subspace $\mathcal{K}_{r+1}$, we end up with a set of time-dependent coupled equations in terms of the first $2r+1$ Hamiltonian moments, which in matrix form can be expressed as

where the vector $\mathbf{c}=(c_{0},c_{1},\cdots,c_{n})^{T}$ and $\dot{\mathbf{c}}$ is its time-derivative. Hamiltonian moments matrices $\mathbf{L}$ and $\mathbf{R}$ are defined as follows,

Eq. (38) is a first order, linear ordinary differential equation (ODE), whose solution depends on the eigenvalues of the matrix $\mathbf{L}^{-1}\mathbf{R}$ that essentially approximate some of the eigenvalues of $\hat{\mathcal{H}}$. In a recent study Guzman and Lacroix (2021), this approach is demonstrate to be able to to target both the ground and excited states for pairing Hamiltonian. Alternatively, as shown in Ref. 71, employing Rayleigh-Ritz technique, $\mathbf{L}$ and $\mathbf{R}$ constitute general eigenvalue problem

with the eigenpair ($\epsilon,\mathbf{v}$) approximate the energy and state vector of the true ground state. To solve Eq. (49) for ($\epsilon,\mathbf{v}$), employ canonical orthogonalization to do the following transformation

where $\mathbf{U}\mathbf{U}^{\dagger}=\mathbf{L}^{-1}$. Since $\mathbf{L}$ is a hermitian that can be diagonalized by a unitary matrix $\mathbf{V}$ with eigenvalues constituting a diagonal matrix $\mathbf{s}$, i.e. $\mathbf{L}=\mathbf{V}\mathbf{s}\mathbf{V}^{\dagger}$, $\mathbf{U}$ can be constructed through $\mathbf{U}=\mathbf{V}\mathbf{s}^{-1/2}$. Similar subspace diagonalization schemes have been adopted in many other hybrid quantum-classical algorithms such as quantum subspace expansion McClean et al. (2017); Colless et al. (2018), and some variants of VQE Parrish et al. (2019); Nakanishi, Mitarai, and Fujii (2019); Huggins et al. (2020).

The imaginary time evolution (ITE) approach has a long history of being a robust computational approach to solve the ground state of a many-body quantum system. To see how it works, let’s first assume the quantum state at time $t$, $|\phi_{t}\rangle$, is now expanded in the eigen-space of $\hat{\mathcal{H}}$, $\{|\psi_{n}\rangle,n=0,1,\cdots\}$, with $E_{n}$’s being the corresponding eigenvalues, then as we show in the previous section $|\phi_{t}\rangle$ can be expressed as

with $c_{n}$’s being the expansion coefficients at the initial time. Now, if we consider a replacement $\tau=it$ ($\tau$ is often called imaginary time) and confine $\tau\geq 0$, Eq. (51) becomes

Comparing Eqs. (51) and (52), we see that the wavefunction is driven from an “oscillating” superposition of the Hamiltonian eigenstates to an “exponential decaying” superposition of the eigenstates with the decay rate proportional to $E_{n}$. More importantly, in the limit of large $\tau$

the ground state $\psi_{0}(x)$ is “screened out” because its exponential decay rate is the smallest in the eigen spectrum of $\hat{\mathcal{H}}$. Briefly speaking, given a trial wavefunction $|\phi\rangle$ that has non-zero overlap with the true ground state wavefunction $\psi_{0}$, $\langle\phi|\psi_{0}\rangle\neq 0$, the ITE approach, in comparison to the RTE of the trial state, guarantees a monotonically decreasing energy functional in $\tau$ that converges to the ground state energy at $\tau\rightarrow\infty$, i.e.

Here the ITE of the trial wavefunction, $|\phi_{t}\rangle$, is defined as

and the first derivative of $E(\tau)$ with respect to $\tau$ is non-positive

where the equality sign holds if and only if $|\phi_{t}\rangle$ is an exact eigenstate of $\hat{\mathcal{H}}$ (e.g. $|\psi_{0}\rangle$).

The development and application of ITE approach targeting the ground state wave function and energy dates back to 1970s, when the similar random-walk imaginary-time technique were developed for diffusion Monte-Carlo methods Davies et al. (1980); Anderson (1975, 1979, 1980). Later on, in the pursuit of non-perturbative analytical tool for a wide variety of Hamiltonian systems that can be systematically improved, various polynomial expansions of Eq. (54) have been studied extensively. For example, a straightforward way is to utilize a finite number of Hamiltonian moments, $\hat{\mathcal{H}}^{n}$, to constitute a truncated Taylor expansion of $E(\tau)$ which we will have a detailed examination in Section IV.

Now if we recursively expand $\frac{1}{s+\hat{\mathcal{H}}}$ in (68) one more time, and introducing another parameter parameter $a_{2}\in\mathbb{R}$, and $a_{2}\neq-s$, we have

with a refined residual term $R_{2}(s,a_{1},a_{2})$ also being non-negative for $\Re(s)>-E_{0}$

Examining the first and second derivatives of $R_{2}(s,a_{1},a_{2})$ with respect to $a_{1}$, $a_{2}$ shows when

$R_{2}$ reaches its local minima.

Note that there is an important feature about $a_{i}$ ($i=1,2$), a little transform of (74) shows

from which it is straightforward to show

therefore all the terms in (70) are non-negative, and

(76) and (77) tell us that by keeping recursively refining $\frac{1}{s+\hat{\mathcal{H}}}$ using identity (68), we are able to get tighter upper bounds to $E_{0}$. This then constitutes the basic idea of PDS approach. For example, if we recursively expand (68) $K$ times ($K\in\mathbb{N},K\geq 2$) with each time introducing a new parameter $a_{i}$, $i=1,\ldots,K$), we will get a more refined non-negative residual term $R_{K}(s,a_{1},\ldots,a_{K})$

and

The $K$-th order PDS formalism, PDS($K$), is then associated with determining the introduced $K$ real parameters $(a_{1}^{(K)},\ldots,a_{K}^{(K)})$ that minimize the value of $R_{K}$ through

which can be alternatively translated to

with each $a_{i}^{(K)}$ ($i=1,\ldots,K$) being an strict upper bound to $E_{0}$, and can be expressed as