# QUBO-based density matrix electronic structure method

Christian F. A. Negre, Alejandro Lopez-Bezanilla, Yu Zhang, Prosper D.
Akrobotu, Susan M. Mniszewski, Sergei Tretiak, and Pavel A. Dub

###### Abstract

Density matrix electronic structure theory is used in many quantum chemistry
methods to “alleviate” the computational cost that arises from directly using
wave functions. Although density matrix based methods are computationally more
efficient than wave functions based methods, yet significant computational
effort is involved. Since the Schrödinger equation needs to be solved as an
eigenvalue problem, the time-to-solution scales cubically with the system
size, and is solved as many times in order to reach charge or field self-
consistency. We hereby propose and study a method to compute the density
matrix by using a quadratic unconstrained binary optimization (QUBO) solver.
This method could be useful to solve the problem with quantum computers, and
more specifically, quantum annealers. The method hereby proposed is based on a
direct construction of the density matrix using a QUBO eigensolver. We explore
the main parameters of the algorithm focusing on precision and efficiency. We
show that, while direct construction of the density matrix using a QUBO
formulation is possible, the efficiency and precision have room for
improvement. Moreover, calculations performing Quantum Annealing with the
D-Wave’s new Advantage quantum processing units is compared with classical
Simulated annealing, further highlighting some problems of the proposed
method. We also show some alternative methods that could lead to a better
performance of the density matrix construction.

††preprint: AIP/123-QED

##  I Introduction

The correct choice of the model or method used to represent the electronic
structure of a chemical system is crucial to get a good description of any
observable (forces, energies, etc.) of such a system. The quality of a model
can be estimated as the ratio between its predictive power and its complexity
(Quality = Predictive Power/Complexity), the predictive power being related to
the amount of information that can be extracted from a calculation, and the
complexity mainly referring to how the time-to-solution scales with the system
size Finnis (2004). Although this definition is vague, it gives us the general
understanding that, in order to optimize a method (increase its quality) it is
necessary to work on approaches that could decrease the complexity, and, at
the same time, avoid losing predictive power; In many cases the predictive
power translates into accuracy. Reducing the complexity is the typical
approach that has been taken given that computational resources have always
been limited due to the ever increasing need of addressing larger systems.
However, with the arrival of quantum computers (an alternative computational
paradigm) and more specifically, quantum annealers (QAs) dwa (2018), it is
possible that more complex formulations (in terms of traditional computation)
of a given problem end up having shorter or even instantaneous time-to-
solution without accuracy resignation. This has been the case of graph
problems including graph partitioning and community detection solved with QAs
that was made possible due to a quadratic unconstrained binary optimization
(QUBO) reformulation of the problem Ushijima-Mwesigwa, Negre, and Mniszewski
(2017); Negre, Ushijima-Mwesigwa, and Mniszewski (2019). For the case of
quantum chemistry methods, the long term vision to reach is the one in which
all the computational burden of the eigenvalue problem and self-consistent
process could be shifted towards an “extremely complex,” yet easily solvable
QUBO problem (See Figure 1). In this paper we take an initial step towards the
aforementioned vision by exploring how feasible it is to construct the density
matrix (DM) using linear algebra QUBO algorithms. Previous work on QAs used in
quantum chemistry has been recently applied to determine the ground and firsts
excited states on the Full Configuration Interaction (CI) method by
diagonalizing the Full CI (FCI) matrix Teplukhin  _et al._ (2020); the
downfolding of the FCI matrix to reduce the size of the effective matrix
Mniszewski  _et al._ (2021); and Time Dependent Density Functional Theory
(DFT) to compute excitations Teplukhin  _et al._ (2021).

$H=H^{0}+V$ if $|\rho-\rho_{old}|<tol$, STOP
$\rho=Cf_{F}(\epsilon)C^{\dagger}$$V=f_{C}(\rho)$$\mathbf{x}^{t}Q\mathbf{x}$
Paradigm shift Figure 1: Scheme showing the paradigm shift hereby proposed. On
the left we represent the traditional self-consistent charge (SCC) approach
with steps that include the correction of the Hamiltonian elements $H$; the
construction of the density matrix $\rho$ from a Fermi function $f_{F}$ of the
Hamiltonian; and the calculation of the effective potential $V$ as a Coulombic
function $f_{C}$ of $\rho$ to be used to reconstruct the Hamiltonian in the
subsequent SCC step. On the right, the alternative QUBO idea is proposed.

The density matrix (DM) formalism is ubiquitous to several methods in
computational chemistry that compute the electronic structure of a molecular
system. Among these methods, the Khon-Sham method within DFT is arguably the
most frequently used of the quantum chemistry methods, preferred over post
Hartree-Fock methods, for its performance in terms of predictive power
relative to its computational cost Vázquez-Mayagoitia  _et al._ (2014);
Fiolhais, Nogueira, and Marques (2003). The foundations of this theory will be
briefly stated below.

The many-body wave function of interacting electrons in a chemical system
resulting in the solution of the Schrödinger equation depends on 3$\times
N_{e}$ spatial coordinates, where $N_{e}$ is the number of electrons in the
system. DFT can be used to overcome this complexity, allowing for the
calculation of observables just by knowing the electronic density of the
system. This theory is based on using the electronic density $n(r)$ ($r$ being
the electron coordinates) as a variable of the total energy functional
$E[n(r)]$ of the system, and the variational principle to minimize the latter.
The theorems of P. Hohenberg and W. Kohn Hohenberg and Kohn (1964) ensure that
the energy, wave function, and all other electronic properties are uniquely
determined by the electron density $n(r)$. Moreover, for any external
potential $v_{ext}(r)$ generated by the nuclear configuration of a molecular
system, there is a bijective relation with a unique density function $n(r)$.
To solve for $n(r)$ the Kohn-Sham method relies on auxiliary wave functions
$\psi_{i}(r)$ that are used to transform the minimization of the energy
functional into the following eigenvalue problem:

$\displaystyle\left(-\frac{\nabla^{2}}{2}+v_{eff}\right)\psi_{i}(r)=\epsilon_{i}\psi_{i}(r)$
(1) $\displaystyle\hat{H}\psi_{i}(r)=\epsilon_{i}\psi_{i}(r)$ (2)

where:

$v_{eff}=v_{ext}(r)+\int\frac{n(r^{\prime})}{|r-r^{\prime}|}dr^{\prime}+\frac{\delta
E_{xc}[n(r)]}{\delta n(r)}$ (3)

Here, $v_{ext}(r)$ is the external potential given by the nuclei; It is the
part of the energy functional that characterizes the chemical system in
question. $E_{xc}$ is the functional that takes into account the energy of
exchange-correlation. It is written as the sum of two components: The kinetic
energy correlation component, which comes from correcting the error by
considering the kinetic energy of non-interacting electrons, and the inter-
electronic interaction containing the non-kinetic correlation and exchange.
The exact form of this functional is not known, but many approximations exists
nowadays with varying accuracies and computational costs Sousa, Fernandes, and
Ramos (2007). In this work we will use the most basic one where the functional
is represented by a local density approximation (LDA) of the form:
$E_{xc}[n(r)]=\int\epsilon_{xc}(n)n(r)d^{3}r$, where $\epsilon_{xc}(n)$ is the
correlation and exchange energy per electron for a uniform gas of density $n$.
Equation 2 is solved iterativelly using equation 3 until reaching self-
consistency between the density $n(r)$ and the effective potential
$v_{eff}(r)$.

When written using a basis set of functions $\\{\phi_{i}\\}_{i\in 1,..,N}$
with $S_{ij}=\int_{r}\phi^{*}_{i}(r)\phi_{j}(r)dr$ and
$H_{ij}=\int_{r}\phi^{*}_{i}(r)\hat{H}\phi_{j}(r)dr$ equation 2 turns into the
following generalized eigenvalue problem:

$H\mathbf{c}=\epsilon_{i}S\mathbf{c}$ (4)

where, $c_{k}$ is the k-th component of the expansion of the wavefunction
$\psi_{i}(r)$ in the chosen basis set. In matrix form equation 4 becomes:

$HC=\epsilon SC$ (5)

which is a generalized matrix eigenvalue problem, where $\epsilon$ is a
diagonal matrix with entries $\epsilon_{ij}=\delta_{ij}\epsilon_{i}$
containing all the eigenvalues of matrix $H$. Finally, $C$ is a matrix
containing all the expansion coefficients for each of the $N$ eigenfunctions
$\psi_{i}(r)$.

Given the expansion coefficients $C$, the DM which is our object of interest
here, is constructed as:

$\rho=Cf(\epsilon)C^{t}$ (6)

where we have used the Fermi distribution function $f_{n}$ for a system in
thermal equilibrium defined as:

$f_{n}(\epsilon)=\frac{2}{1+\exp(\frac{\epsilon-\epsilon_{F}}{k_{b}T})}$ (7)

where $k_{b}$ is the Boltzman constant, $\epsilon_{F}$ is the Fermi level (or
electrochemical potential) and $T$ is the system’s temperature. For practical
purposes, when temperature is low,
$f_{n}(\epsilon)=2\Theta(\epsilon-\epsilon_{F})$, leading to:

$\rho=\sum_{i}^{occ}2\mathbf{c}_{i}\mathbf{c}_{i}^{T}$ (8)

The DM can be hence directly constructed by summing the self outer products of
each eigenvector (in order of increasing energy $\epsilon_{i}$) up to the
number of occupied states ($occ=N_{e}/2$). Given the DM, many properties can
be computed since the expectation value of any operator $\hat{A}$ can be
calculated as $\langle A\rangle=tr(\rho A)$. For instance, the trace of the DM
gives the total number of electrons of the system, and the trace of $\rho H$
provides the total electronic energy also known as band energy of the system
Szabo and Ostlund (1996). Another important property is that we can easily
compute $n(r)$ as follows:

$\begin{split}n(r)=\sum_{k,l}\phi_{k}(r)\rho_{kl}\phi_{l}(r)\end{split}$

which can be reused in equation 3 to recompute $\hat{H}$. Other methods such
as Hartree-Fock, Tight-Binding, and Semiempirical are also based on
constructing the single-particle DM Elstner  _et al._ (1998); McWeeny (1956,
1960); Finnis (2004). From what was explained above, matrices $\epsilon$ and
$C$ arise from diagonalizing the matrix representation of $\hat{H}$ in the
aforementioned basis set. This diagonalization step in all the density matrix
based methods is typically the bottleneck of the whole calculation, and scales
as $\mathcal{O}(N^{3})$, where $N$ is the number of elements in the basis set
(which scales linearly with the system size). Over the years, computational
chemists have dedicated lots of effort to develop algorithms that could solve
the density matrix with $\mathcal{O}(N)$ scaling Goedecker and Scusseria
(2003); Bowler and Gillan (1999); Niklasson (2002); leading to the so called
“order n methods.” Some of them have also adapted algorithms to the ever
evolving classical computer architectures Mniszewski  _et al._ (2015);
Finkelstein  _et al._ (2021). The problem with all these order n approaches is
that they require a trade-off between accuracy and efficiency. The faster we
want the method to be, the lower will be the accuracy of our results. In this
paper we analyze the tractibility of quantum annealing as an alternative
computational paradigm to get accurate results at a lower computational cost
for DM based quantum chemistry methods.

##  II Methods

In this section we will describe the methods used, and in particular we will
focus on our recently developed high precision QUBO eigensolver Krakoff,
Mniszewski, and Negre (2021).

###  II.1 Quantum annealers

Quantum Annealers (QAs) are specific adiabatic quantum computing devices
particularly suitable to solve combinatorial optimization problems dwa (2020).
In these devices the exploration of the combinatorial space is given by the
delocalization and degree of superposition of the initial quantum state that
decays to the state of minimum energy during a process called quantum
annealing. The D-Wave QA is an architecture that uses this quantum annealing
technique to solve for the ground state of an Ising Model Hamiltonian that is
derived from a combinatorial optimization problem. In other words, any problem
to be solved with the D-Wave QA will need to be reformulated to be solved as
an Ising Model Hamiltonian whose energy is evaluated as follows:

$E(\mathbf{s})=\sum_{i}h_{i}s_{i}+\sum_{i<j}J_{i,j}s_{i}s_{j}$ (9)

for which $s_{i}\in\\{-1,1\\}$ are magnetic spin variables; $h_{i}$ are local
programmable fields; and $J_{ij}$ are programmable field coupling
interactions. An Ising model Hamiltonian can be easily transformed into its
boolean equivalent using the following transformation:
$\mathbf{s}=2\mathbf{q}-\mathbf{1}$, where $\mathbf{1}=(1,\ldots,1)$. The
boolean equivalent problem that is obtained is referred to as a QUBO problem
and it is expressed as follows:

$E(\mathbf{q})=\mathbf{q}^{T}Q\mathbf{q}=\sum_{i}q_{i}Q_{ii}+\sum_{i\neq
j}q_{i}q_{j}Q_{ij}.$ (10)

with $q_{i}\in\\{0,1\\}$.

###  II.2 High-precision QUBO-based eigensolver

The solution of the extremal eigenvectors with high accuracy (eigenvector
corresponding to minimal or maximal eigenvalue) was recently possible thanks
to a new QUBO formulation of the problem. Krakoff, Mniszewski, and Negre
(2021). This method builds on top of the quantum annealer eigensolver (QAE)
work of Teplukhin et. al. Teplukhin  _et al._ (2020). The main difference
between this high-precision QUBO eigensolver method and the previous QAE, is
that it solves the problem as a fixed point method instead of using binary
search Hoffman and Frankel (2018). The new method also involves a steepest
descent search together with a systematic adjustment of the length of the
precision vector to guarantee high accuracy convergence. The precision vector
defined as follows is used to convert numbers from binary to real
representation:

$\mathbf{v}=(-1,2^{-1},2^{-2},...,2^{-m})$

where $m$ is the number of bits used in the representation. In order to
transform a number from binary $q$ to real $x$ representation we need to
perform $x=\mathbf{v}^{T}\mathbf{q}$, where $\mathbf{q}$ is the vectorized
form of $q$. Note that with this representation we cover the $[-1,1)$ segment
where the binary representation $q$ of $x$ will satisfy:
$-1\equiv(1,0,0,0,...,0)\leq q<(0,1,1,1,...,1)\equiv 1$. A real field
multiplication of $x\times y$ has the binary equivalent
$q_{x}^{T}(\mathbf{v}\otimes\mathbf{v}^{T})q_{y}$. A particular example of
this multiplication is given in the Appendix.

With the help of this precision vector and the previously introduced
multiplication, any quadratic equation $\mathbf{x}^{T}B\mathbf{x}$ for which
we need to find an optimal $\mathbf{x}$ with $B\in\mathbb{R}^{N\times N}$,
will have binary equivalent
$\mathbf{q}^{T}(B\otimes(\mathbf{v}\otimes\mathbf{v}^{T}))\mathbf{q}$, where
$(B\otimes(\mathbf{v}\otimes\mathbf{v}^{T}))\in\mathcal{B}^{mN\times mN}$,
being $\mathcal{B}^{mN\times mN}$ the ${mN\times mN}$ binary field vector
space. Any linear term $\mathbf{r}^{T}\mathbf{x}$ with
$\mathbf{r}\in\mathbb{R}^{N}$ for which we need to find an optimal
$\mathbf{x}$, will have a binary equivalent
$\mathbf{q}^{T}\mathrm{Diag}(\mathbf{r}I\otimes\mathbf{v}))\mathbf{q}$, where
$I$ is the $N\times N$ identity matrix, and $\mathrm{Diag}()$ is an operator
that, given $\mathbf{a}$ in $\mathcal{B}^{mN}$ builds an operator $A$ in
$\mathcal{B}^{mN\times mN}$ with $A_{ij}=\delta_{ij}a_{i}$.

Given the Hamiltonian ($H$) and overlap ($S$) matrices, every eigenpair
$(\lambda,\mathbf{x})$ will satisfy: $(H\mathbf{x}-\lambda S\mathbf{x})=0$,
and, in particular, the lowest eigenpair ($\lambda_{min}$, $\mathbf{x}_{min}$)
will be an absolute minimum of:

$(H\mathbf{x}-\lambda S\mathbf{x})^{T}(H\mathbf{x}-\lambda S\mathbf{x})$ (11)

Since $H-\lambda S$ is positive semidefinite, we can instead minimize:

$\mathbf{x}^{T}(H-\lambda S)\mathbf{x}\quad\equiv\quad\mathbf{q}^{T}(H-\lambda
S)\otimes(\mathbf{v}\otimes\mathbf{v}^{T})\mathbf{q}$ (12)

Since $\lambda_{min}$ is not known, equation 12 needs to be solved
iterativelly as a fixed point. A good initial value for $\lambda$ happens to
be the average eigenvalue of $H$ computed as $\mathrm{tr}\left(H\right)/N$
Krakoff, Mniszewski, and Negre (2021). The method converges once the changes
in $\lambda_{min}$ are below a certain preset tolerance.

In order to achieve a higher precision in the calculation of both
$\mathbf{x}_{min}$ and $\lambda_{min}$, we have introduced a steepest descent
phase also solved as a QUBO problem. Setting the previously converged
$\mathbf{x}_{min}$ as an initial guess $\mathbf{x}^{0}_{min}$, we Taylor
expand $f(\mathbf{x})=\mathbf{x}^{T}(H-\lambda S)\mathbf{x}$ around
$\mathbf{x}^{0}_{min}$ to obtain:

$\displaystyle f(\mathbf{x})$ $\displaystyle=f(\mathbf{x}^{0}_{min})+\nabla
f\cdot\delta+\delta^{T}\frac{\mathcal{H}(f)}{2}\delta+o(||\delta||^{3})$
$\displaystyle\delta$ $\displaystyle:=\mathbf{x}-\mathbf{x}^{0}_{min}$

where $\nabla f$ and $\mathcal{H}(f)$ are the gradient and Hessian of $f$
respectivelly. For the case of $f(\mathbf{x})$, the gradient and Hessian
becomes $\nabla f=2(H-\lambda S)(\mathbf{x}^{0}_{min})^{T}$, and
$\mathcal{H}(f)=2(H-\lambda S)$ respectivelly. Hence, we can compute a high-
precision approximation to the lowest eigenpair as:

$\displaystyle\mathbf{x}$
$\displaystyle=\mathbf{x}^{0}_{min}+\underset{\mathbf{\delta}}{\text{argmin}}\left(\nabla
f\cdot\mathbf{\delta}+\mathbf{\delta}\cdot\mathcal{H}\mathbf{\delta}\right)$

In this case, $\mathbf{\delta}$ converges to a fixed point, and, once
$\mathbf{\delta}$ does not change anymore, the length of the precision vector
is shortened to 10% of its previous value. This procedure guarantees
convergence to any arbitrary desired precision.

For every “next” eigenpair we want to compute, the previous eigenvalue
$\mathbf{x}_{min}\equiv\mathbf{x}_{i-1}$ of $H$ needs to be “pushed out” of
the eigenspectrum by using the following transformation:

$H=H-wS\mathbf{x}_{i-1}\mathbf{x}^{T}_{i-1}S$ (13)

where $w$ is an estimation for the “eigenspectrum width.” A good estimation of
$w$ can be calculated as $w=h-l$ where $h$ and $l$ are the highest and lowest
Gershgorin bounds of $H$ respectively.

Using the previously explained procedure to compute eigenvectors, the density
matrix is constructed as explained before; using a straightforward summation
of eigenvector outer products up to the number of occupied states. This is:

$\rho=\sum_{k}^{nocc}\mathbf{x}_{k}\mathbf{x}^{T}_{k}$ (14)

A python like pseudocode for the full algorithm is detailed in Pseudocode 1.
The inputs are the Hamiltonian $H$, Overlap $S$, the occupation $nocc$, and
the target precision $prec_{T}$. The output is the DM $\rho$.

def get_rho($H,S,nocc,prec_{T}$):

#Eigenspectrum width estimation

$w$ = get_Gershgorin_width($H$)

for i in range(nocc):  #Loop over occ states

$\lambda$ = $\mathrm{tr}\left(H\right)/N$ #Initial estimate of $\lambda$

#Initial phase

while ($err>tol$):

$\mathbf{x}=$solve($\mathbf{x}^{T}(H-\lambda S)\mathbf{x}$)  #Send to Annealer

$\mathbf{x}\leftarrow\mathbf{x}/(\mathbf{x}^{T}S\mathbf{x})$ #Normalize

$\lambda=\mathbf{x}^{T}H\mathbf{x}$ #Compute new eigenvalue

$err=\mathbf{x}^{T}(H-\lambda S)\mathbf{x}$ #Compute error

#Descent phase

$prec=0.1$

while ($prec>prec_{T}$):

$\mathbf{\delta}=$solve($\nabla f\delta+\delta\mathcal{H}\delta$)  #Send to
Annealer

$\mathbf{x}\leftarrow\mathbf{x}+\mathbf{\delta}$

$\mathbf{x}\leftarrow\mathbf{x}/(\mathbf{x}^{T}S\mathbf{x})$ #Normalize

$\lambda=\mathbf{x}^{T}H\mathbf{x}$ #Compute new eigenvalue

$err=\mathbf{x}^{T}(H-\lambda S)\mathbf{x}$ #Compute error

if $(err>err_{old}):$

$\mathbf{v}\leftarrow 0.1\mathbf{v}$

$prec=0.1prec$

$err_{old}=err$

$H\leftarrow H-wS\mathbf{x}^{T}\mathbf{x}S$ #Shift eigenvalue up

$\rho\leftarrow\rho+2\mathbf{x}^{T}\mathbf{x}$ #Construct partial DM

return $\rho$

Pseudocode 1 Full DM construction algorithm using high-precision QUBO
eigensolver. Description of the variables follows. $H$: Input Hamiltonian;
$prec_{T}$: Target precision; $nooc$: Number of occupied states; $\lambda$:
Minimal eigenvalue being computed; $N$: Size of $H$; $err$: Convergence error;
$tol$: Tolerance for the initial phase; $S$: Overlap matrix; $\mathbf{x}$:
Minimal eigenvector being computed; $prec$: Convergence precision; $\nabla f$:
Gradient of $f$; $\mathcal{H}$: Hessian matrix; $\rho$: Output DM

Pseudocode 1 implements a straightforward way of computing the density matrix
using the high-precision QUBO eigensolver. From this pseudocode we can
immediately identify some important properties. We see that the matrix
operations in the inner loops are $\mathcal{O}(N^{2})$ unless they can be
solved using sparse matrix formats to get linear scaling. The other important
point is that the two inner loops (initial and descent phases) have to be
solved $nocc$ times, where $nocc$ is the number of occupied states that
depends on the number of electrons which scales linearly with the number of
orbitals $N$. From this, we can deduce that even if the solution of the
quadratic problem takes no time, the overall scaling of the algorithm gives
back the “undesired” $\mathcal{O}(N^{3})$ scaling. Provided the matrix
operations in the inner loop could be solved with $\mathcal{O}(N)$ scaling
methods, and the quadratic problem takes no time, we will have at best an
$\mathcal{O}(N^{2})$ scaling. Even though this seems to be an undesired
result, it is a first direct approach to solve the DM using a QUBO eigensolver
and deserves some further study.

###  II.3 DFT calculations

We used the DFT-based code SIESTA Ordejón, Artacho, and Soler (1996); Artacho
_et al._ (1999) to obtain all the molecules’ Hamiltonians in a tight-binding
like representation. SIESTA is a first-principles electronic structure code
for molecules and solids which represents wavefunctions as a combination of
atomic-like orbitalsDemkov  _et al._ (1995). The Kohn-Sham equations are
solved self-consistently in a linear combination of atomic orbitals basis set
and they are restricted to a fix-size cell condition. This framework allows us
to obtain a description of the operator associated with the system energy
employing a reduced number of orbitals in the basis set with complete
multiple-zeta and polarized bases, depending on the required accuracy. In this
study, all first-principles calculations were performed with a minimal
single-$\zeta$ basis for each atom. Description of the interaction between
atoms was conducted within the local density approximation approachKohn and
Sham (1965) for the exchange-correlation functional. The integration over the
Brillouin zone (BZ) was performed using a Monkhorst sampling in $\Gamma$
point. The radial extension of the orbitals had a finite range with a kinetic
energy cutoff of 50 meV. A separation of 20 Å in a cubic simulation box
prevents virtual periodic molecules from interacting. We used different
numbers of molecules to get the Hamiltonians ($H$) and overlap ($S$) matrices
needed to test our algorithm. We also computed $H$ and $S$ for benzene in
order to test for resiliency regarding degeneracy of eigenstates.

##  III Results and discussion

We first studied the propagation of the errors with the number of occupied
states. The error in the calculated DM is computed as follows:

$error=||\rho-\rho^{ex}||_{1}=\sum_{ij}|\rho_{ij}-\rho^{ex}_{ij}|$ (15)

where $\rho^{ex}$ is the DM obtained using regular diagonalization in double
precision which we hereby set as our standard. In order to reach chemical
accuracy, the error needs to satisfy $error\leq
1.0[kCal/mol]/\max_{lm}(H_{lm})$ as explained in the Appendix. Results showing
scaling of error and iteration with system size, occupation, and degeneracy
were performed using classical SA algorithm as implemented in the dwave-neal
code for the Ocean API dwa . A comparison with QA running on the D-Wave
Advantage 4.1 system is offered at the end of this section. From the
Pseudocode 1, we can deduce that the error committed on the calculations of
the very first computed eigenpairs is susceptible to be propagated towards the
higher eigenpairs, hence, making the calculations of DM with higher
occupations less accurate. Figure 2 shows the error as a function of the
occupation number for a single water molecule with six total atomic orbitals.
Three different precision values were used to show that this problem happens
regardless of the precision used. In general, we notice that the higher we set
the precision for the computation of the eigenvectors, the less error is
committed in the construction of the DM across different values of
occupations. The larger the number of occupied orbitals, the larger the error
that is committed in the computation of the DM.

Figure 2: Error as a function of the number of occupied states for a single
water molecule with $N=6$. Results for three different precision values are
shown.

The problem of error propagation worsens with system size. As it is shown in
Figure 3, the error gets larger when the total number of orbitals increases,
regardless of the precision and number of bits that are used. We notice that,
the error increases with the system size which is a sign that a higher
precision will be required for larger systems to reach chemical accuracy.

Figure 3: Error as a function of system size (number of orbitals).
Calculations using two and four bits of precision are shown.

Figure 4 shows the number of iterations (total QUBO solutions) that are needed
to construct the DM for two and four bits. The total number of iterations is
computed by adding up the number of iterations needed to achieve the desired
precision for the calculation of each eigenvector. From Figure 3 it does not
seem that much is gained by increasing the number of bits, however, as it is
shown in Figure 4, we notice that less iterations are needed to reach
convergence when using a larger number of bits.

Figure 4: Number of iterations as a function of the system size (number of
orbitals). Calculations using two and four bits of precision are shown.

Finally, we computed the number of iterations needed to achieve the desired
precision to compute every single state for the benzene molecule. This
particular case has been selected since, due to its high symmetry (see
molecular representation on the inset of figure 5), the electronic structure
of benzene is characterized by having several groups of degenerate molecular
orbitals (orbitals that are very close in energy). In general, the closer two
molecular orbitals are in energy, the more the degeneracy between them. Figure
5 shows the number of iterations for computing a single eigenvector as a
function of the proximity in energy to the next eigenvector. We notice that a
group of molecular orbitals that are close in energy to the next ($<$ 0.1 eV)
need significantly more iterations to be computed at the desired precision.

Figure 5: Number of required iterations to compute each eigenpair as a
function of the minimal energy difference from other orbitals. The benzene
molecule is shown in the inset with carbon and hydrogen atoms colored in cyan
and white respectivelly.

This point is a significantly limiting factor since the electronic structure
of several chemical systems including semiconductors and metals will have a
high degree of degeneracy. For semiconductors, degeneracy will most likely
happen around the band edges, while for metals it will happen around the Fermi
level.

##  IV Quantum vs Simulated Annealing

In this section we present a comparison between the results obtained using SA
with different sweeps and QA performed with the new D-Wave’s Advantage 4.1
machine Pau Farre (2021). The chain strength for QA was calculated as 30% of
the largest absolute value in the QUBO matrix D-Wave (2020). This helps ensure
minimal chain breaks contributing to erroneous solutions. The default anneal
time of 20 microseconds was used. The system DM calculation was performed over
samples Hamiltonians taken from different configurations over a 100 Molecular
Dynamics (MD) steps. We have extracted $H$ from 20 configurations evenly
spaced over the course of the MD trajectory. The DM was computed using an
eigenpair precision of $10^{-5}$. From Figure 6 we can see that, on average,
the results using QA have a larger error. Moreover, we have observed that more
iterations are needed to converge. In the case of SA a larger number of sweeps
tend to converge to results with lower errors. Means for the error between SA
and QA are statistically different even for small values of sweeps.

$0.4$$0.6$$0.8$$1$$1.2$$1.4$$1.6$$\cdot 10^{-2}$SA-10SA-100SA-1000DW5000QError
Figure 6: Figure showing different box plots for the error commited with
different methods. Quantum Annealer (red); and Simulating Annealing with
different number of sweeps (black)

One of the mayor issues is the fact that the QUBO matrices have a high dynamic
range which causes a “resolution problem.” When the number of bits is large
there are many small QUBO entries that could have large relative error at the
moment of embedding but that might have a large influence in the results. One
could scale the QUBO matrix Q so that the minimun entry has a value of 1.0,
but this introduces the need of having very strong chain strengths that could
introduce yet another source of error. Moreover, we have observed that some
chain length are too long due to the high connectivity or matrix density which
also depends on the number of bits. Other types of optimization problems have
also revealed no advantage using QA vs SA Koshikawa  _et al._ (2021).

##  V Alternative approaches

In this section we will briefly describe some alternative approaches that
could potentially overcome the issues that were found by using a direct method
to construct the DM using a QUBO eigensolver. Three different approaches are
hereby suggested together with a discussion of possible issues they might
have. 1) Using a QUBO based linear solver on the $\mathcal{O}(N)$ method
recently proposed by L. Anderson Andersson (2019). This algorithm shifts the
computational cost of solving $\rho$ to the solution of a linear system
($A\mathbf{x}=b$) using Conjugate Gradient method which can be replaced by:

$\mathbf{x}=\underset{\mathbf{x}}{\text{argmin}}\left((A\mathbf{x}-b)^{2}\right)$

and can be easily formulated as a QUBO problem and solved on a QA. The problem
that this method could have is the fact that there might be several iterations
(QUBO solves) needed to reach the desired precision of $\mathbf{x}$. This will
add yet another inner loop into the algorithm.

2) We can use Green Functions ($G(z)$) method and find $\rho$ by integrating
over the complex contour (with variable $z\in\mathbb{C}$) containing the
eigenspectrum of $H$. For this proposition we would need to have a method to
invert a complex matrix using a QUBO formulation. Provided we can compute the
inverse of a matrix using a QUBO formulation, $\rho$ would be determined as:

$\displaystyle\rho=\int_{C}G(z)dz$ (16)

where $G(z)=(H-zI)^{-1}$ and C is the semicircle above the real axis.
Inverting a matrix $A$ can be simply done by solving $N$ linear systems of the
form $A\mathbf{x}_{i}-\mathbf{e}_{i}=0$, with $\mathbf{e}_{i}$ being the i-th
canonical vector, and $\mathbf{x}_{j}$ the j-th column of $A^{-1}$. This can
be easily formulated as a QUBO provided both, the solution and the matrix to
be inverted are real. If this is not the case, as for the Green Function
method here proposed, the formulation becomes more complicated and, as far as
we know, no previous work has addressed this issue.

3) Finally, we could think about solving for $\rho$ “all at once” from a QUBO
problem where $\rho$ would be written as a long vector $\mathbf{x}$ containing
$N\times N$ entries, such that: $\mathbf{x}_{k}=\rho_{ij}$ where
$i=\left\lfloor{k/N}\right\rfloor+1$, and
$j=k-N\left\lfloor{k/N}\right\rfloor$; where $\left\lfloor{.}\right\rfloor$ is
the operator that takes the “floor” value of a real number. In this case, the
function to be minimized will be the total energy

$E=\mathrm{tr}\left(\rho H\right)\equiv\mathbf{1}^{T}(H\otimes I)\mathbf{x}$
(17)

where in this case $\mathbf{1}_{k}=\delta_{ij}$, where
$i=\left\lfloor{k/N}\right\rfloor+1$, and
$j=k-N\left\lfloor{k/N}\right\rfloor$. Although equation 17 can be easily
formulated as a QUBO problem, the energy needs to be minimized under certain
constraints given by the properties of $\rho$; some of which can be extremely
difficult to formulate as a QUBO. Moreover, the size of the QUBO problem to be
solved will scale as $N^{2}\times N^{2}$, probably opening up other
bottlenecks in the algorithm even if the annealing process require no time.

##  VI Conclusion

We have demonstrated that the DM can be computed by using a QUBO based
eigensolver and summing up the self outer products of the eigenvectors up to
the number of occupied states. Although the results obtained with this direct
approach are encouraging, many issues still needs to be addressed. We reported
on the problem of the $\mathcal{O}(N^{2})$ operation inside loops that compute
the egenpairs which leads back to an $\mathcal{O}(N^{3})$ scaling. By
experimenting with different occupation numbers, we have seen that the error
rapidly increases with the occupation number regardless of the precision at
which we compute each eigenvector. We have also seen that the error increases
with the system size regardless of the number of bits used to compute the
eigenpairs. We observed that by increasing the number of bits the number of
iterations to converge decreases significantly. We have identified a problem
where the number of iterations to reach precision increases for degenerate
states. Finally we showed that QA does not show any advantage when compared to
SA for computing DM within this QUBO formulation and more research needs to be
done to fully understand the QA vs SA comparison.

We have also proposed other alternative methods to compute the DM and analyze
some possible issues. These methods include using linear scaling algorithms
based on linear systems that could be solved using QUBO solvers; using Green
functions obtained from a QUBO solver that then could be used to compute the
DM by a close integral surrounding the eigenspectrum through the complex
plane; and solving the full density matrix “all at once” from a QUBO problem
using a vectorized form of the DM. Even if all these alternative methods have
challenges there is still an opportunity for further improvements. An in depth
study of these other methods will be the subject of future work.

##  VII Acknowledgments

Research presented in this article was supported by the Laboratory Directed
Research and De-velopment (LDRD) program of Los Alamos National Laboratory
(LANL) under project number 20200056DR. This research was also supported by
the U.S. Department of Energy (DOE) National Nuclear Security Administration
(NNSA) Advanced Simulation and Computing (ASC) program at LANL. We acknowledge
the ASC program for providing the support for accessing D-Wave’s Advantage 4.1
computing resource. LANL is operated by Triad National Security, LLC, for the
National Nuclear Security Administration of U.S. Department of Energy
(Contract No. 89233218NCA000001). Assigned: LA-UR-22-20271

##  VIII Appendix

###  VIII.1 Multiplication in binary representation

Multiplying $x$ and $y$ in the binary representation leads to the following
operation:

$x\times y=q^{T}_{x}(\mathbf{v}^{T}\otimes\mathbf{v})q_{y}$

An example with 3-bits follows. Given we have: $x=0.5\equiv(0,1,0)$, and
$y=0.25\equiv(0,0,1)$

$\displaystyle x\times y=0.125$
$\displaystyle=(0,1,0)\begin{pmatrix}-1&-0.5&-0.25\\\ -0.5&0.25&0.125\\\
-0.25&0.125&0.0625\end{pmatrix}\begin{pmatrix}0\\\ 0\\\ 1\end{pmatrix}=0.125$

###  VIII.2 Density matrix upper bound error

The error in $\rho$ is computed as: $\Delta\rho=||\rho-\rho_{ex}||$; where
$\rho_{ex}$ is the exact DM constructed using regular diagonalization methods.
To ensure chemical accuracy we need to ensure that the total energy error
remains below 1.0 kCal/mol. This means: $\Delta E<1.0kCal/mol$.

$\Delta E=\frac{\partial\mathrm{tr}\left(\rho
H\right)}{\partial\rho}\Delta\rho+\frac{\partial\mathrm{tr}\left(\rho
H\right)}{\partial H}\Delta H$

Assuming no error is commited in the calculation of the Hamiltonian, we get:

$\Delta
E=\sum_{lm}\partial\frac{\sum_{ik}\rho_{ik}H_{ki}}{\partial\rho_{lm}}\Delta\rho_{lm}$

$\Delta
E=\sum_{lm}H_{lm}\Delta\rho_{lm}\leq\max_{lm}(H_{lm})\sum_{lm}|\Delta\rho_{lm}|$

meaning that:

$\Delta
E\leq\max_{lm}(H_{lm})||\rho-\rho_{ex}||_{1}=\max_{lm}(H_{lm})\Delta\rho$

From the last inequality, we deduce that the maximun tolerated error in $\rho$
needs to satisfy:

$\Delta E\leq\max_{lm}(H_{lm})\Delta\rho\leq 1kCal/mol$

then,

$\Delta\rho\leq\frac{1kCal/mol}{\max_{lm}(H_{lm})}$

## References

  * Finnis (2004) M. W. Finnis,  _Interatomic forces in condesed matter_ , Oxford series on materials modelling (Oxford University Press, 2004).
  * dwa (2018) “D-wave systems,” (2018), https://www.dwavesys.com/home.
  * Ushijima-Mwesigwa, Negre, and Mniszewski (2017) H. Ushijima-Mwesigwa, C. F. A. Negre, and S. M. Mniszewski, “Graph partitioning using quantum annealing on the D-Wave system,” in  _Proceedings of the Second International Workshop on Post Moores Era Supercomputing_ , PMES’17 (Association for Computing Machinery, New York, NY, USA, 2017) pp. 22–29.
  * Negre, Ushijima-Mwesigwa, and Mniszewski (2019) C. F. A. Negre, H. Ushijima-Mwesigwa, and S. M. Mniszewski, “Detecting multiple communities using quantum annealing on the D-Wave system,” PLOS ONE (2019).
  * Teplukhin  _et al._ (2020) A. Teplukhin, B. K. Kendrick, S. Tretiak, and P. A. Dub, “Electronic structure with direct diagonalization on a dwave quantum annealer,” Sci. Rep. 10, 20753 (2020).
  * Mniszewski  _et al._ (2021) S. M. Mniszewski, P. A. Dub, S. Tretiak, P. M. Anisimov, Y. Zhang, and C. F. A. Negre, “Reduction of the molecular hamiltonian matrix using quantum community detection,” Scientific Reports 11 (2021).
  * Teplukhin  _et al._ (2021) A. Teplukhin, B. K. Kendrick, S. M. Mniszewski, Y. Zhang, A. Kumar, C. F. A. Negre, P. M. Anisimov, S. Tretiak, and P. A. Dub, “Computing molecular excited states on a D-Wave quantum annealer,” Sci. Rep. 11, 18796 (2021).
  * Vázquez-Mayagoitia  _et al._ (2014) Á. Vázquez-Mayagoitia, W. Scott Thornton, J. R. Hammond, and R. J. Harrison, “Quantum chemistry methods with multiwavelet bases on massive parallel computers,” (2014).
  * Fiolhais, Nogueira, and Marques (2003) C. Fiolhais, F. Nogueira, and M. Marques,  _A Primer in Density Functional Theory_ (Springer-Verlag Berlin Heidelberg, 2003).
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, “Inhomogeneous electron gas,” Phys. Rev. 136, 864 – 870 (1964).
  * Sousa, Fernandes, and Ramos (2007) S. F. Sousa, P. A. Fernandes, and M. J. Ramos, “General performance of density functionals,” J. Phys. Chem. A 111, 10439–10452 (2007).
  * Szabo and Ostlund (1996) A. Szabo and N. S. Ostlund,  _Modern Quantum Chemistry_ (Dover Publications Inc., 1996).
  * Elstner  _et al._ (1998) M. Elstner, D. Porezag, G. Jungnickel, M. Haugk, T. Frauenheim, S. Suhai, and G. Seifer, “Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties,” Phys. Rev. B 58, 7260 – 7267 (1998).
  * McWeeny (1956) R. McWeeny, “The density matrix in self-consistent field theory .1. iterative construction of the density matrix,” Proc. R. Soc. London Ser. A-Math 235, 496 (1956).
  * McWeeny (1960) R. McWeeny, “Rev mod phys,” Rev. Mod. Phys. 32, 335 (1960).
  * Goedecker and Scusseria (2003) S. Goedecker and G. E. Scusseria, “Linear scaling electronic structure methods in chemistry and physics,” Comput. Sci. Eng. 5, 14–21 (2003).
  * Bowler and Gillan (1999) D. R. Bowler and M. J. Gillan, “Density matrices in O(N) electronic structure calculations: theory and applications,” Comput. Phys. Commun. 120, 95–108 (1999).
  * Niklasson (2002) A. M. N. Niklasson, Phys. Rev. B 66, 155115 (2002).
  * Mniszewski  _et al._ (2015) S. M. Mniszewski, M. J. Cawkwell, M. E. Wall, J. Mohd-Yusof, N. Bock, T. C. Germann, and A. M. N. Niklasson, “Efficient parallel linear scaling construction of the density matrix for born-oppenheimer molecular dynamics,” J. Chem. Theory Comput. 11, 4644 (2015).
  * Finkelstein  _et al._ (2021) J. Finkelstein, J. S. Smith, S. M. Mniszewski, K. Barros, C. F. A. Negre, E. H. Rubensson, and A. M. N. Niklasson, “Mixed precision fermi-operator expansion on tensor cores from a machine learning perspective,” J. Chem. Theory Comput. (2021).
  * Krakoff, Mniszewski, and Negre (2021) B. Krakoff, S. M. Mniszewski, and C. F. A. Negre, “A QUBO algorithm to compute eigenvectors of symmetric matrices,” (2021), arXiv:2104.11311 [cs.ET] .
  * dwa (2020) “Welcome to d-wave,”  (2020), online; accessed February 2, 2021.
  * Hoffman and Frankel (2018) J. D. Hoffman and S. Frankel,  _Numerical Methods for Engineers and Scientists_ (CRC Press, 2018).
  * Ordejón, Artacho, and Soler (1996) P. Ordejón, E. Artacho, and J. M. Soler, “Self-consistent order-$n$ density-functional calculations for very large systems,” Physical Review B 53, R10441–R10444 (1996).
  * Artacho  _et al._ (1999) E. Artacho, D. Sánchez-Portal, P. Ordejón, A. García, and J. M. Soler, “Linear-scaling ab-initio calculations for large and complex systems,” physica status solidi (b) 215, 809–817 (1999).
  * Demkov  _et al._ (1995) A. A. Demkov, J. Ortega, O. F. Sankey, and M. P. Grumbach, “Electronic structure approach for complex silicas,” Physical Review B 52, 1618–1630 (1995).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, “Self-consistent equations including exchange and correlation effects,” Phys. Rev. 140, A1133–A1138 (1965).
  * (28) “dwave-neal: An implementation of a simulated annealing sampler for general ising model graphs in c++ with a dimod python wrapper,” .
  * Pau Farre (2021) C. M. Pau Farre, “The advantage system: Performance update,” Tech. Rep. 14-1054A-A (DWAVE, 2021).
  * D-Wave (2020) D-Wave, “Programming the D-Wave QPU: Setting the chain strength,” D-Wave Systems Whitepaper 2020-04-14 , 11 (2020).
  * Koshikawa  _et al._ (2021) A. S. Koshikawa, M. Ohzeki, T. Kadowaki, and K. Tanaka, “Benchmark test of black-box optimization using D-Wave quantum annealer,” (2021), arXiv:2103.12320 [cond-mat.stat-mech] .
  * Andersson (2019) L. Andersson, “Linear-scaling recursive expansion of the Fermi-Dirac operator,” (2019).


* Liu et al. [2015] Zhao Liu, Jianke Zhu, Jiajun Bu, and Chun Chen.  A survey of human pose estimation: The body parts parsing based methods.  _Journal of Visual Communication and Image Representation_ , 32:10–19, July 2015.
  * Mahmoud et al. [2013] Marwa Mahmoud, Louis-Philippe Morency, and Peter Robinson.  _Automatic multimodal descriptors of rhythmic body movement_.  ACM, December 2013.
  * [41] Mental Health Foundation.  https://www.mentalhealth.org.uk/, 2019.  Accessed 2019/06/04.
  * Ory et al. [2013] Marcia G Ory, SangNam Ahn, Luohua Jiang, Kate Lorig, Phillip Ritter, Diana D Laurent, Nancy Whitelaw, and Matthew Lee Smith.  National Study of Chronic Disease Self-Management: Six-Month Outcome Findings.  _Journal of Aging and Health_ , 25(7):1258–1274, September 2013.
  * Ozdas et al. [2000] A Ozdas, R G Shiavi, S E Silverman, M K Silverman, and D M Wilkes.  Analysis of fundamental frequency for near term suicidal risk assessment.  In _IEEE International Conference on Systems, Man, and Cybernetics_ , pages 1853–1858. IEEE, 2000.
  * [44] Peer2Peer Cambridge.  https://www.peer2peer-cambridge.org/, 2019.  Accessed 2019/06/01.
  * Peng et al. [2005] Hanchuan Peng, Fuhui Long, and C Ding.  Feature selection based on mutual information criteria of max-dependency, max-relevance, and min-redundancy.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 27(8):1226–1238, August 2005.
  * Scherer et al. [2014] Stefan Scherer, Giota Stratou, Gale Lucas, Marwa Mahmoud, Jill Boberg, Jonathan Gratch, Albert Skip Rizzo, and Louis-Philippe Morency.  Automatic audiovisual behavior descriptors for psychological disorder analysis.  _Image and Vision Computing_ , 32(10):648–658, October 2014.
  * Sharif Razavian et al. [2014] Ali Sharif Razavian, Hossein Azizpour, Josephine Sullivan, and Stefan Carlsson.  CNN Features Off-the-Shelf: An Astounding Baseline for Recognition.  pages 806–813, 2014.
  * Shukla et al. [2017] Parul Shukla, Kanad K Biswas, and Prem K Kalra.  Recurrent Neural Network Based Action Recognition from 3D Skeleton Data.  In _2017 13th International Conference on Signal-Image Technology & Internet-Based Systems (SITIS)_, pages 339–345. IEEE, 2017.
  * Simon et al. [2017] Tomas Simon, Hanbyul Joo, Iain A Matthews, and Yaser Sheikh.  Hand Keypoint Detection in Single Images Using Multiview Bootstrapping.  _CVPR_ , cs.CV, 2017.
  * Sobin and Sackeim [2019] Christina Sobin and Harold A Sackeim.  Psychomotor symptoms of depression.  _academia.edu_ , 2019.
  * Song et al. [2013] Yale Song, Louis-Philippe Morency, and Randall Davis.  Learning a sparse codebook of facial and body microexpressions for emotion recognition.  In _the 15th ACM_ , pages 237–244, New York, New York, USA, 2013\. ACM Press.
  * Spitzer et al. [2006] Robert L Spitzer, Kurt Kroenke, Janet B W Williams, and Bernd Löwe.  A Brief Measure for Assessing Generalized Anxiety Disorder.  _Archives of Internal Medicine_ , 166(10):1092–1097, May 2006.
  * Srivastava et al. [2003] Sanjay Srivastava, Oliver P John, Samuel D Gosling, and Jeff Potter.  Development of personality in early and middle adulthood: Set like plaster or persistent change?  _Journal of Personality and Social Psychology_ , 84(5):1041–1053, 2003.
  * Su et al. [2017] Benyue Su, Huang Wu, and Min Sheng.  Human action recognition method based on hierarchical framework via Kinect skeleton data.  In _2017 International Conference on Machine Learning and Cybernetics (ICMLC)_ , pages 83–90. IEEE, 2017.
  * Syed et al. [2017] Zafi Sherhan Syed, Kirill A Sidorov, and A David Marshall.  Depression Severity Prediction Based on Biomarkers of Psychomotor Retardation.  _AVEC@ACM Multimedia_ , pages 37–43, 2017.
  * Wei et al. [2017] Shenghua Wei, Yonghong Song, and Yuanlin Zhang.  Human skeleton tree recurrent neural network with joint relative motion feature for skeleton based action recognition.  In _2017 IEEE International Conference on Image Processing (ICIP)_ , pages 91–95. IEEE, 2017.
  * Wei et al. [2016] Shih-En Wei, Varun Ramakrishna, Takeo Kanade, and Yaser Sheikh.  Convolutional Pose Machines.  In _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 4724–4732. IEEE, 2016.
  * Yang et al. [2017a] Le Yang, Dongmei Jiang, Xiaohan Xia, Ercheng Pei, Meshia Cédric Oveneke, and Hichem Sahli.  Multimodal Measurement of Depression Using Deep Learning Models.  _AVEC@ACM Multimedia_ , pages 53–59, 2017a.
  * Yang et al. [2017b] Le Yang, Hichem Sahli, Xiaohan Xia, Ercheng Pei, Meshia Cédric Oveneke, and Dongmei Jiang.  Hybrid Depression Classification and Estimation from Audio Video and Text Information.  _AVEC@ACM Multimedia_ , pages 45–51, 2017b.
  * Zhang et al. [2016] Hong Bo Zhang, Qing Lei, Bi Neng Zhong, Ji Xiang Du, and Jia Lin Peng.  A Survey on Human Pose Estimation.  _Intelligent Automation and Soft Computing_ , 22(3):483–489, July 2016.