# A two-circuit approach to reducing quantum resources for
the quantum lattice Boltzmann method

Sriharsha Kocherla  School of Computational Science & Engineering, Georgia
Institute of Technology, Atlanta, GA 30332, USA  Austin Adams  School of
Computer Science, Georgia Institute of Technology, Atlanta, GA 30332, USA
Zhixin Song  School of Physics, Georgia Institute of Technology, Atlanta, GA
30332, USA  Alexander Alexeev  George W. Woodruff School of Mechanical
Engineering, Georgia Institute of Technology, Atlanta, GA 30332, USA  Spencer
H. Bryngelson<EMAIL_ADDRESS>School of Computational Science & Engineering,
Georgia Institute of Technology, Atlanta, GA 30332, USA  Daniel Guggenheim
School of Aerospace Engineering, Georgia Institute of Technology, Atlanta, GA
30332, USA

###### Abstract

Computational fluid dynamics (CFD) simulations often entail a large
computational burden on classical computers. At present, these simulations can
require up to trillions of grid points and millions of time steps. To reduce
costs, novel architectures like quantum computers may be intrinsically more
efficient at the appropriate computation. Current quantum algorithms for
solving CFD problems use a single quantum circuit and, in some cases, lattice-
based methods. We introduce the a novel multiple circuits algorithm that makes
use of a quantum lattice Boltzmann method (QLBM). The two-circuit algorithm we
form solves the Navier–Stokes equations with a marked reduction in CNOT gates
compared to existing QLBM circuits. The problem is cast as a stream
function–vorticity formulation of the 2D Navier–Stokes equations and verified
and tested on a 2D lid-driven cavity flow. We show that using separate
circuits for the stream function and vorticity lead to a marked CNOT
reduction: 35% in total CNOT count and 16% in combined gate depth. This
strategy has the additional benefit of the circuits being able to run
concurrently, further halving the seen gate depth. This work is intended as a
step towards practical quantum circuits for solving differential equation-
based problems of scientific interest.

00footnotetext:  Code available at: https://github.com/comp-physics/QLBM-
frugal

##  I Introduction

Quantum computers can achieve complexity improvements by leveraging quantum
mechanics principles such as superposition and entanglement. One example of
quantum computing overcoming the limitations of computation on classical
computers is Shor’s algorithm [1], which computes the prime factorization of
integers in logarithmic time, providing an exponential speedup. Efforts are
being made to explore the possibility of using quantum resources to speed up
other important tasks. Example tasks include linear systems solvers [2, 3],
Monte Carlo methods [4], and machine learning algorithms [5, 6]. With the
recent progress of quantum hardware, various quantum algorithms have been
purposed and examined for optimization problems [7], quantum chemistry [8],
and finance [9].

An emerging use case of quantum computers is solving computational fluid
dynamics (CFD) problems [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22,
23, 24]. For example, some techniques for solving the incompressible
Navier–Stokes equations rely on an efficient on a sparse linear solver. The
Harrow–Hassidim–Lloyd (HHL) algorithm [2] promises an exponential speedup on
such tasks, assuming the state preparation process and information readout are
efficient. Lapworth [25] proposes a hybrid CFD solver based on the SIMPLE
(Semi-Implicit Method for Pressure Linked Equations) algorithm and HHL.
However, the number of operations HHL requires exceeds current quantum
hardware capability. Instead, more proposals focus on using variational
quantum algorithms (VQA) as near-term replacements for HHL to solve linear
systems without any guaranteed complexity advantage. Demirdjian  _et al._ [26]
solves 1D advection-diffusion equation with Carleman linearization technique
and variational quantum linear solver [27].

Mesoscale methods for solving the Navier–Stokes equations operate above the
molecular level but below the continuum scale. These methods provide indirect
solutions by solving the Boltzmann equation, which characterizes the
statistical behavior of a system of particles. Two methods that can solve the
Boltzmann equation are lattice gas automata [28, 29] and the lattice Boltzmann
method (LBM) [30, 31]. Traditional LBM algorithms are suited for classical
algorithms, though some algorithms for quantum devices have been developed.
For example, previous art developed quantum algorithms for the lattice gas
[32, 33, 34]. Itani and Succi [35] demonstrate a Carleman linearization for
the collision term and further work in Itani  _et al._ [36] shows unitary
evolution for both collision and streaming operators. Further works have
pushed these methods forward [37, 38].

Budinski [39] introduced a quantum lattice Boltzmann method (QLBM) for the
advection–diffusion equation. This method simulates particles across a grid of
lattice sites and extracts macroscopic quantities from mesoscopic simulations.
Budinski [40] adapted the quantum circuit approximation to the Navier–Stokes
equations. Using the stream function–vorticity formulation of the
Navier–Stokes equations, the pressure term is removed, and the problem is
reduced to an advection–diffusion equation and a Poisson equation. The
strategy of these works involves a single quantum circuit for both equations.

This work formulates a two-circuit 2D QLBM to solve the incompressible Navier-
Stokes equations. The method is validated against the classical LBM solution
to the 2D lid-driven cavity flow problem. The quantum resources required for
this QLBM approach are compared to other approaches. Section II describes the
classical LBM theory for simulating the advection–diffusion equation. Section
III.1 shows a single-circuit quantum implementation of the LBM for solving the
advection–diffusion equation. In section IV, a description of the stream
function-vorticity formulation of the Navier–Stokes equations and work on a
two-circuit quantum implementation for solving the 2D Navier–Stokes equations
are provided. Section V presents verification against the classically-solved
LBM simulation for the lid-driven cavity problem and quantum resource
estimation. Section VI serves to clarify the contributions of the method
presented in this manuscript and its limitations.

##  II Classical lattice Boltzmann method

The LBM uses a grid of cells (shown in fig. 1) that contain particle
distributions. The particles move into neighbor cells during each time step
according to the Boltzmann equation [41, 42]:

$\displaystyle f_{\alpha}(\bm{r}+e_{\alpha}\Delta t,t+\Delta
t)=(1-\epsilon)f_{\alpha}(\bm{r},t)+\epsilon f_{\alpha}^{\text{(eq)}}+\Delta
w_{\alpha}S,$ (1)

where $e_{\alpha}$ is the velocity of a particle in link $\alpha$,
$f_{\alpha}$ is the particle distribution along each $\alpha$, $t$ is time,
$\bm{r}$ is the cell position, $S$ is the source term, $w_{\alpha}$ is the
proportion of particles streaming in link $\alpha$, and $\epsilon=\Delta
t/\tau$ where $\tau$ is a relaxation time. A D$n$Q$m$ scheme denotes an
$n$-dimensional lattice grid with $m$ links, where a link is a direction in
which particles can be propagated. The weights for the D2Q5 scheme (a 2D grid
with 5 links per cell) are $w_{0}=2/6$ and $w_{2,3,4}=1/6$ following the
directions shown in fig. 1 [43].

Figure 1:  Illustration of a D2Q5 LBM lattice with a streaming particle and
link weights $w_{\alpha}$.

We apply a new formulation to two PDE problems, the advection–diffusion and
incompressible Navier–Stokes equations. Using the stream function–vorticity
formulation, the Navier–Stokes equations reduce into two equations that
closely resemble the advection–diffusion equation, so the latter will be
examined first.

The advection–diffusion is

$\displaystyle\frac{\partial\phi}{\partial t}+c\frac{\partial\phi}{\partial
x}=D\frac{\partial^{2}\phi}{\partial x^{2}},$ (2)

which represents advection as $c({\partial\phi}/{\partial x})$ and diffusion
as $D({\partial^{2}\phi}/{\partial x^{2}})$. Here, $c$ and $D$ are advection
and diffusion coefficients and $\phi(t,x)$ is a scalar concentration field.

For the advection–diffusion equation of (2), the equilibrium distribution
function is

$\displaystyle
f_{\alpha}^{\text{(eq)}}(\bm{r},t)=w_{\alpha}\phi(\bm{r},t)\left(1+\frac{\bm{e_{\alpha}}\cdot\bm{c}}{c_{s}^{2}}\right),$
(3)

where $\phi(\bm{r},t)$ corresponds to the lattice site at lattice position
$\bm{r}$ = $(x_{i},y_{j})$, $\bm{c}$ is the advection velocity vector, and
$c_{s}=1/\sqrt{3}$ is the speed of sound for the D2Q5 scheme [44]. The
notation $x_{i}$ and $y_{j}$ denotes the $i$-th cell in the $x$ direction and
the $j$-th cell in the $y$ coordinate direction.

Following standard practice [30], the relaxation time $\tau$ relates to the
diffusion constant $D$ with

$\displaystyle D=c_{s}^{2}\left(\tau-\frac{\Delta t}{2}\right).$ (4)

Herein we use $D=1/6$, $c_{s}=1/\sqrt{3}$, and $\tau=\Delta t$. This gives
$\varepsilon=\tau/\Delta t=1$ and

$\displaystyle f_{\alpha}(\bm{r}+\bm{e_{\alpha}}\Delta t,t+\Delta
t)=f_{\alpha}^{\text{(eq)}}$ (5)

simplifies as

$f_{\alpha}(\bm{r}+\bm{e_{\alpha}}\Delta t,t+\Delta
t)=w_{\alpha}\phi(\bm{r},t)\left(1+\frac{\bm{e_{\alpha}}\cdot\bm{c}}{c_{s}^{2}}\right),$
(6)

which includes the collision (3) and the streaming steps (5). The
concentration field is thus

$\phi(\bm{r},t)=\sum_{\alpha=0}^{N-1}f_{\alpha}(\bm{r},t),$ (7)

summing the particle distributions across all link directions $\alpha$.

##  III Quantum lattice Boltzmann method

###  III.1 Advection–diffusion equation

A quantum circuit can implement the lattice Boltzmann method, which involves a
similar approach to classical computation. Specifically, the quantum circuit
performs the collision operator followed by the streaming of particles and
subsequent recalculation of macroscopic variables. Boundary conditions are
applied at the end of each time step. No special care is needed when boundary
conditions are periodic; the left and right shift gates $R$ and $L$
automatically propagate boundary conditions to the designated lattice site for
each link $\alpha$.

It is important to note that encoding input is a bottleneck in the algorithm.
For the current study, the amplitude encoding technique described in the work
of Shende  _et al._ [45] is employed. Amplitude encoding is part of the Qiskit
toolkit [46], though it is also a generic encoding algorithm with an expensive
gate count.

We organize the quantum circuit for the advection–diffusion equation using a
D2Q5 scheme with qubits organized into 4 registers, which are groupings of
qubits. Registers $r_{0}$ and $r_{1}$ contain $\log M$ qubits, where $M$ is
the number of lattice sites in each dimension. Register $d$ has
$\lceil\log_{2}Q\rceil$ qubits, where $Q$ is the number of LBM links $\alpha$.
Register $a$ holds a single ancilla qubit necessary for applying a non-unitary
collision operator.

####  III.1.1 Encoding input

At the start of each time step, given a distribution $\phi(\bm{r},t)$, define
$\phi(\alpha,\bm{r})$ to be the value of $\phi$ at time $t$ and position
$\bm{r}$ at link direction $\alpha$. Given a discretized concentration

$\displaystyle\phi(\alpha,\bm{r})=[\phi(0,0),\phi(0,1),\phi(0,2),\dots,\phi(4,M-2),\phi(4,M-1),\phi(4,M)],$
(8)

the initial statevector $\ket{\psi_{0}}$ is

$\displaystyle\ket{\psi_{0}}=\ket{0}_{a}\otimes\frac{1}{||\phi||}\sum_{i=0}^{2M-1}\phi(\alpha,\bm{r})\ket{i}.$
(9)

This normalizes $\phi(\alpha,\bm{r})$, and the initial qubit states follow
from amplitude encoding. In fig. 2, register $r_{0}$ stores the data at each
link direction $\alpha$, and it can be retrieved by specifying the link
$\alpha$ in register $d$.

Applying the collision operator to the initial statevector $\ket{\psi_{0}}$ is
equivalent to multiplying the $\phi(\alpha,\bm{r})$ with weight coefficients

$\displaystyle
w_{\alpha}\left(1+\frac{\bm{e_{\alpha}}\cdot\bm{c}}{c_{s}^{2}}\right),$ (10)

which follow from (3). This strategy is discussed further in the next
subsection.

####  III.1.2 Collision operator

Figure 2:  QLBM circuit for an advection–diffusion time step. The collision
operator $A=(C_{1}+C_{2})/2$, and shift operators $R$, $L$, $D$, and $U$
propagate particles in each link direction $\alpha$. In panel (c), a Hadamard
[47] gate is applied to the qubits in register $d$.

The collision step (fig. 2 (a)) computes the equilibrium distribution function
$f_{a}^{\text{(eq)}}$, which requires computing the proportion of the
distribution $\phi$ in each link $\alpha$. The collision operator entails
applying the coefficient matrix $A$ to the current statevector
$|\psi_{0}\rangle$.

The coefficient matrix $A$ is not unitary, so it cannot be directly translated
into a quantum gate. An important strategy was described by Budinski [40]: A
linear combination of unitary matrices split $A$ into two matrices $C_{1}$ and
$C_{2}$, related to the original matrix as

$\displaystyle C_{1,2}=A\pm i\sqrt{I-A^{2}}.$ (11)

As $A=(C_{1}+C_{2})/2$, an operation with $A$ is computed via block encoding
[48, 49], where $C_{1}$ and $C_{2}$ are unitary, but $A$ is, in general, not.

Figure 3:  A block encoding of $B$ = $(V+W)/2$. Here, $V$ and $W$ are unitary
matrices.

The circuit of fig. 3 evolves an input statevector $\ket{\psi}$ according to
the unitary

$\displaystyle U=\frac{1}{2}\begin{bmatrix}V+W&V-W\\\ V-W&V+W\end{bmatrix},$
(12)

where $B=(V+W)/2$. Thus, $B$ is a subblock of the block matrix $U$.

Post-selection selects quantum states for specific measurement outcomes. Here,
we use post-selection to measure the result after the collision operator for a
$0$ ancilla (after the block encoding).

The collision matrix is

$\displaystyle A=\begin{bmatrix}k_{1}I_{n}&0\\\
0&k_{2}I_{n}\end{bmatrix},\quad\text{where}\quad
k_{\alpha}=w_{\alpha}\left(1+\frac{\bm{e_{\alpha}}\cdot\bm{c}}{c_{s}^{2}}\right)$
(13)

are the link coefficients described by (10). Panel (a) of fig. 2 uses a linear
combination of unitaries to apply matrix $A$ to the input. The Hadamard gates
are used for the block encoding process along with $C_{1}$ and $C_{2}$
operations, which are derived from the unitary matrices of (11). The
coefficients representing the proportion of particles in each link $\alpha$
are $k_{1}$ and $k_{2}$. The collision matrices are

$\displaystyle C_{1,2}=\begin{bmatrix}\exp(\pm i\cos^{-1}(k_{1}))I_{n}&0\\\
0&\exp(\pm i\cos^{-1}(k_{2}))I_{n}\end{bmatrix}.$ (14)

The collision operator $A$ transforms statevector $\ket{\psi_{0}}$ via a
linear combination of $C_{1}$ and $C_{2}$, but requires an ancilla qubit $a$,
which stores orthogonal data $(C_{1}-C_{2})/2$ where the ancilla is $\ket{1}$.
The orthogonal data was ignored through post-selection. The result of this
linear combination is

$\displaystyle\ket{\psi_{1}}=\frac{1}{||\phi||}\sum_{i}a_{i,i}\phi_{i,i}\ket{i},$
(15)

which encodes the post-collision values for each link direction $\alpha$, for
which the ancilla is $\ket{0}$.

###  III.2 Particle streaming

The streaming step propagates particles in each link $\alpha$ to the
neighboring site. Figure 5 shows the shift operators $R_{n}$ and $L_{n}$,
which are controlled on link qubits $d$ and stream particles to neighboring
lattice sites. The shift matrices are permutation matrices as

$\displaystyle R_{n}=\begin{bmatrix}0&0&\cdots&0&1\\\ 1&0&\cdots&0&0\\\
0&1&\cdots&0&0\\\ \vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&\cdots&1&0\end{bmatrix}\quad\text{and}\quad
L_{n}=\begin{bmatrix}0&1&\cdots&0&0\\\ \vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&\cdots&1&0\\\ 0&0&\cdots&0&1\\\ 1&0&\cdots&0&0\end{bmatrix},$ (16)

and are both unitary matrices.

The resulting statevector $\ket{\psi_{2}}$ has a distribution shifted to a
neighbor lattice site corresponding to its link value $\alpha_{i}$. The right
and left shift gates $R$ and $L$, shown in fig. 4, are controlled on the link
qubits $d$ to act on the part of the statevector corresponding to a
distribution $\alpha$. Figure 5 shows the up and down operators $U$ and $D$,
which are implemented via the $L$ and $R$ operators applied via $r_{1}$.

(a) Left shift operator $L$. (b) Right shift operator $R$.

Figure 4:  Circuit decomposition for (a) left shift and (b) right shift
operators. These gates are applied to qubits in $r_{1}$ for left and right
shifts and qubits in $r_{2}$ for up- and down-shifts.  Figure 5:  Streaming
for a two-dimensional lattice grid. The streaming step uses right and left
shift operators $R$ and $L$ to shift distributions in their respective link
directions $\alpha_{1}$ and $\alpha_{2}$ by controlling the gates on qubits in
each link register $d$.

###  III.3 Macroscopic variable retrieval

We retrieve the distribution $\phi(\bm{r},t)$ by summing
$f_{\alpha}^{\mathrm{(eq)}}$ over the link directions $\alpha$. In terms of
quantum operations, this is accomplished via the application of Hadamard gates
to each of the qubits in link registers $d$, as shown in fig. 2 (c) [47]. A
Hadamard gate $H$ applied to a statevector $\ket{\psi}$, resulting in

$\displaystyle H|\psi\rangle=\frac{1}{\sqrt{2}}\begin{bmatrix}1&1\\\
1&-1\end{bmatrix}\begin{bmatrix}a\\\
b\end{bmatrix}=\frac{1}{\sqrt{2}}\begin{bmatrix}a+b\\\
a-b\end{bmatrix}=\frac{a+b}{\sqrt{2}}|0\rangle+\frac{a-b}{\sqrt{2}}|1\rangle.$
(17)

Thus, when the Hadamard gate is applied to a link qubit in $d$, it stores the
sum of the two amplitudes as the $\ket{0}$ amplitude and the difference as the
$\ket{1}$ amplitude. So, the sum can be post-selected by ignoring $\ket{1}$
measurements. This means that the Hadamard gates will sum the distributions
but introduce a factor of $1/\sqrt{2}$ per gate. This prefactor is post-
processed out of the computation by multiplication of a factor of
$\sqrt{2}^{\log Q}$ where $Q=|\alpha|$ is the number of link distributions. At
each time step, we retrieve the circuit result via state tomography, which
must be used to extract the first $M$ lattice site elements.

##  IV Quantum lattice Boltzmann method for the Navier–Stokes Equations

###  IV.1 Stream function–vorticity formulation

The incompressible 2D Navier–Stokes equations in Cartesian coordinates are

$\displaystyle\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial
x}+v\frac{\partial u}{\partial y}$ $\displaystyle=-\frac{\partial p}{\partial
x}+\frac{1}{\mbox{\text{Re}}}\left(\frac{\partial^{2}u}{\partial
x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}\right),$ (18)
$\displaystyle\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial
x}+v\frac{\partial v}{\partial y}$ $\displaystyle=-\frac{\partial p}{\partial
y}+\frac{1}{\mbox{\text{Re}}}\left(\frac{\partial^{2}v}{\partial
x^{2}}+\frac{\partial^{2}v}{\partial y^{2}}\right),$ (19)

where $u$ and $v$ are the velocity components in the $x$ and $y$ coordinate
directions, $p$ is the pressure, and Re is the Reynolds number, which is the
ratio of inertial to viscous effects [50].

Taking the curl of the above Navier–Stokes equations recasts them in the so-
called vorticity–stream function formulation, removing the pressure term $p$
and yielding

$\displaystyle\frac{\partial^{2}\psi}{\partial
x^{2}}+\frac{\partial^{2}\psi}{\partial y^{2}}$ $\displaystyle=-\omega,$ (20)
$\displaystyle\frac{\partial\omega}{\partial
t}+u\frac{\partial\omega}{\partial x}+v\frac{\partial\omega}{\partial y}$
$\displaystyle=\frac{1}{\mbox{\text{Re}}}\left(\frac{\partial^{2}\omega}{\partial
x^{2}}+\frac{\partial^{2}\omega}{\partial y^{2}}\right).$ (21)

In this formulation, (20) and (21) use vorticity $\omega$ and stream function
$\psi$ instead of directional speeds $u$ and $v$. The velocity vector is thus
$\bm{u}=\\{u,v\\}$. The stream function relates to the directional velocities
as

$\displaystyle\frac{\partial\psi}{\partial
x}=u\quad\text{and}\quad\frac{\partial\psi}{\partial y}=-v.$ (22)

So, (20) is a Poisson equation in stream function $\phi$ and (21) is a
advection–diffusion equation in vorticity $\omega$.

####  IV.1.1 Lattice-based representation

With the stream function–vorticity formulation, the collision, streaming, and
macro lattice stages follow as

$\displaystyle f_{\alpha}^{\text{(eq)}}(\bm{r},t)$
$\displaystyle=w_{\alpha}\omega(\bm{r},t)\left(1+\frac{\bm{e}_{\alpha}\cdot\bm{u}}{c_{s}^{2}}\right),$
(23) $\displaystyle f_{\alpha}(\bm{r}+\bm{e}_{\alpha}\Delta t,t+\Delta t)$
$\displaystyle=f_{\alpha}^{\text{(eq)}},$ (24) $\displaystyle\omega(\bm{r},t)$
$\displaystyle=\sum_{\alpha}f_{\alpha}(\bm{r},t),$ (25)

respectively. These stages, buttressed via the circuits of section III.1,
enable the vorticity $\omega$ computation.

The equilibrium distribution function for the Poisson equation
($\nabla^{2}\psi=-\omega$) is
$g_{\alpha}^{\text{(eq)}}(\bm{r},t)=w_{\alpha}\psi(\bm{r},t)$. The streaming
and macro steps match those of section III.2 and section III.3, but the source
term $S=-\omega$ is added during the collision step (a). Thus, the relaxation
operator is

$\displaystyle g_{\alpha}(\bm{r}+\bm{e}_{\alpha}\Delta t,t+\Delta
t)=g_{\alpha}^{\text{(eq)}}+\Delta w_{\alpha}S,$ (26)

and macro retrieval equation

$\displaystyle\psi(\bm{r},t)=\sum_{\alpha}g_{\alpha}(\bm{r},t).$ (27)

####  IV.1.2 Boundary conditions

The validation problem for the proposed method is a 2D lid-driven cavity flow.
The spatial domain is $\Omega\in x,y$ with lengths $L_{x}$ and $L_{y}$ and
boundary $\partial\Omega$. The stream function $\psi$ is constant along the
boundaries, with $\psi=0$ used here. For the lattice Boltzmann method,

$\displaystyle\psi=\sum_{\alpha=0}^{N_{l}-1}g_{\alpha}^{\text{(eq)}}=g_{0}+g_{1}+g_{2}+g_{3}+g_{4},\quad\text{so,}\quad
g_{\partial\Omega}=-\sum_{\alpha|\alpha\neq\partial\Omega}g_{\alpha}.$ (28)

Defining the vorticity expression (20) in terms of the stream function and
expanding it in its Taylor series gives

$\displaystyle\omega_{i,N}=-2\left(\frac{\psi}{\Delta y^{2}}+\frac{U}{\Delta
y}\right),$ (29)

along the boundaries $\partial\Omega$, where $U$ is wall-parallel velocity of
the top wall. For a stationary wall, $U=0$. The wall equilibrium distribution
in the direction of the wall is

$\displaystyle
g(x,y)_{\partial\Omega}=-{\sum_{\alpha|\alpha\neq\partial\Omega}g_{\alpha}}-2\left(\frac{-\psi}{\Delta
y^{2}}+\frac{U}{\Delta y}\right).$ (30)

To implement (30), the matrix

$\displaystyle B=\begin{bmatrix}0&\cdots&0\\\ \vdots&I_{N-2}&\vdots\\\
0&\cdots&0\end{bmatrix}$ (31)

is applied to the statevector $|\phi\rangle$, where $B$ is of size $N\times
N$, where $N=\log(n)$ in each dimension and $I_{n}$ denotes the size $n$
identity matrix. So, $g_{\partial\Omega}=0$, while retaining other
distribution values, enforces the boundary condition for the stream function
circuit $|\psi_{\text{b.c.}}\rangle=|0\rangle^{\otimes n}$.

Due to its non-unitary properties, the vorticity circuit boundary conditions
are applied via a linear combination of $B$. This linear unitary combination
is

$\displaystyle D_{1,2}=B\pm i\sqrt{I_{N}-B^{2}},$ (32)

following section III.1.2.

Figure 6:  A D2Q5 vorticity circuit, including boundary conditions. Compared
with the advection–diffusion circuit in fig. 2, this has an additional qubit
$b$ required to hold boundary conditions and an additional operation required
to compute them. However, if boundary conditions are computed classically, the
vorticity and advection–diffusion circuits are identical.

###  IV.2 Quantum circuits

####  IV.2.1 Vorticity circuit

Figure 7:  Collision portion of the stream function circuit. An additional
qubit $s$ stores the source term, and an additional Hadamard gate on qubit s
after the block encoding adds the source term to the link distributions
$\alpha_{i}$ before streaming.

The vorticity circuit computes the vorticity $\omega$ for the current time
step. Note the vorticity circuit in fig. 6 and advection–diffusion circuit in
fig. 2 match, except for boundary conditions. This occurs because the
vorticity equation (described in eq. 21) follows the same form as the
diffusion equation. The Navier–Stokes algorithm’s vorticity circuit is
identical to the advection–diffusion circuit for classically-computed boundary
conditions.

If boundary conditions are included, an extra qubit $b$ stores them. Boundary
conditions require additional computation, as enforcing such conditions is not
a unitary operation. For this, the linear combination of unitaries is used
[51], described further in section III.1.2, The input to this circuit for the
no-boundary version is the previous time step’s vorticity $\omega_{t-1}$. When
using boundary conditions, the circuit input includes pre-computed boundary
conditions in the boundary qubit $b$, computed following the descriptions in
section IV.1.2.

####  IV.2.2 Stream function circuit

The stream function circuit is the other quantum circuit used in this two-
circuit model. The D2Q5 stream function circuit shown in fig. 7 is referred to
for additional context. This circuit is similar in character to the
advection–diffusion circuit with classical boundary conditions but an
additional source. A key difference is the requirement of an extra qubit $s$
to store the source term, $S=-\omega$, as shown in fig. 7.

The $C_{1}$ and $C_{2}$ gates also operate on the source term, and the
Hadamard gate on qubit $s$ adds this source to the equilibrium stream function
distribution function $g_{\alpha}^{\mathrm{(eq)}}$. The same qubit-addition
process from section IV.2.1 applies. Boundary conditions require an additional
$b$ qubit store and a boundary gate $B$, which is not unitary and again is
implemented via a linear combination of matrices following (12). The previous
stream function, $\psi_{t-1}$, and source term, $S=-\omega_{t-1}$, serve as
the input to the stream function circuit. The boundary conditions are computed
according to section IV.1.2 for our quantum boundary condition variant.

##  V Simulations and results

###  V.1 Advection–diffusion equation

Example results for the QLBM circuit applied to the advection–diffusion
equation are shown for two exemplar cases. In the 1D case, D1Q2 and D1Q3
lattice schemes simulate a dense concentration of $\rho=0.2$ at source
$x_{i}=10$ undergoing advection and diffusion with uniform advection velocity
of $c=1/5$ and diffusion coefficient $D=1/6$. The validation problem for the
2D case follows a diffusing concentration $\rho=0.3$ at source
$(x_{i},y_{j})=(4,4)$ and $0.1$ elsewhere, solved via a D2Q5 scheme.

|
---|---
(a) D1Q2 and D1Q3 | (b) D2Q5
Figure 8: Quantum LBM (a) D1Q2, D1Q3, and (b) D2Q5 results to the
advection–diffusion equation. The initial condition in (a) is a point source
of $\rho=0.2$ at $x=10$ and $\rho=0.1$ in all other sites. The initial
condition in (b) is a point source of $\rho=0.3$ at $(x,y)=(4,4)$ and
$\rho=0.1$ elsewhere.

Figure 8 shows the results of these simulations, showing the checkerboard
pattern that arises in D1Q2 due to the distribution moving wholly into
neighboring areas, resulting in half of the lattice sites having zero
particles. The algorithmic deficiency is remedied via the D1Q3 lattice scheme
with an $w_{0}=2/3$ of the particles, which remain in their current lattice
site. Figure 8 (a) shows that the D1Q2 and D1Q3 schemes return similar
results, but D1Q3 resolves the checkerboarding problem.

Figure 8 (b) shows the results of the 2D test problem solved with the D2Q5
scheme. We observe a source with an advection velocity in the positive $x$ and
$y$ directions, advecting in the direction of the velocity and diffusing
outwards from an initial source at $(x_{i}=4,y_{j}=4)$. The solution behaves
as expected and agrees with classical LBM results.

###  V.2 Navier–Stokes equations

Figure 9 shows isocontours of the stream function for a lid-driven cavity
problem. The problem serves to verify the two-circuit QLBM against a classical
implementation of the lattice Boltzmann method.

Figure 9:  A 2D lid-driven cavity flow. Steady-state stream function
isocontours are shown as labeled. The initial conditions are $\psi=0$, $w=0$,
and $U=1$.

To define relative errors between the quantum and classical lattice
algorithms, we denote

$\displaystyle\psi_{i,j}=\psi_{(x_{i},y_{j})}\quad\text{and}\quad\omega_{i,j}=\omega_{(x_{i},y_{j})},$
(33)

where $x_{i}=i\Delta x$ and $y_{j}=j\Delta y$. The local $L_{1}$ relative
error between the classical and two-circuit QLBM Navier–Stokes solver is,
thus,

$\displaystyle\varepsilon_{\psi;i,j}=\frac{\psi_{i,j}^{\text{classic.}}-\psi_{i,j}^{\text{quant.}}}{\psi_{i,j}^{\text{classic.}}}\quad\text{and}\quad\varepsilon_{\omega;i,j}=\frac{\omega_{i,j}^{\text{classic.}}-\omega_{i,j}^{\text{quant.}}}{\omega_{i,j}^{\text{classic.}}}.$
(34)

(a) Relative error $\varepsilon_{\omega}$ for vorticity circuit (b) Relative
error $\varepsilon_{\psi}$ for stream function circuit

Figure 10:  Relative error between the classical LBM and two-circuit quantum
methods for (a) vorticity and (b) stream function.

The relative errors are shown for the cavity problem in fig. 10. The two-
circuit QLBM has a close agreement with the results of the classical LBM when
simulated with the statevector simulator in Qiskit’s SDK [46], as shown in
fig. 10. The state vector simulator perfects exact tomography but is generally
limited to relatively small simulations due to exponentially increasing memory
requirements with qubit numbers.

####  V.2.1 Quantum resource estimation and improvement

| CNOT Gates | Circuit Depth
---|---|---
Single-circuit QLBM | 25 | 58
Stream function | 4.3 | 9.4
Vorticity | 12 | 39
Stream function without boundaries | 4.2 | 15
Vorticity without boundaries | 3 | 6.5
Table 1: Quantum resource estimation (all counts are in units of $10^{4}$) for
a D2Q5 algorithm with lattice size $64\times 64$.

Implementing the two-circuit method for solving the Navier–Stokes equations
using quantum lattice-based algorithms shows quantum resource efficiency
advantages over the single-circuit method. The difference in accuracy between
the classical lattice Boltzmann method (LBM) and the two-circuit
implementation is negligible. The stream function and vorticity circuits can
be run concurrently for each time step, saving computational time. Other
improvements include circuits with classical boundary conditions, decreasing
circuit depth, and the number of CNOT gates compared to a single-circuit
strategy. CNOT gates are slower than single qubit gates such as Hadamard
gates, so CNOT gate count is used here to compare quantum resource costs.

The Qiskit transpiler translates the circuits of previous sections into single
qubit and CNOT gates. The results of table 1 show that pre-computing the
boundary conditions using classical methods reduces the gate count by at least
$35\%$. Even when these circuits are computed in serial, the combined quantum
resources for stream function and vorticity circuits remain notably lower than
those of a single-circuit QLBM. Of course, parallel implementations are
natural with the two-circuit approach.

##  VI Conclusion

This work presents improvements on the quantum lattice Boltzmann method for
solving the two-dimensional Navier–Stokes equations. Testing on the validation
problem shows that the error between the quantum and classical LBM methods is
negligible. Runtime and quantum resources markedly decrease from a previous
single-circuit implementation of an otherwise similar algorithm [40]. The
number of qubits required scales as $O(\log M)$ where $M$ is the size of the
grid, as opposed to the linear scaling of previous related works such as Yepez
[34], which makes the prospect of near-term implementations on quantum
hardware potentially viable, though not tested here.

The CNOT count of the vorticity and stream function circuits scales as
$O(\alpha M)$. However, the encoding step is an efficiency bottleneck,
occurring after each time step. Readout of the statevector and subsequent re-
encoding decreases the efficiency of this method and is an important
consideration for future work. Using the two-circuit quantum lattice Boltzmann
method simplifies the quantum algorithm implementation and reduces the
required quantum resources.

## Acknowledgements

We thank Professors J. Yepez and L. Budinski for fruitful discussions of this
work. This work was supported in part by DARPA grant number HR0011-23-3-0006,
the Georgia Tech Quantum Alliance, and a Georgia Tech Seed Grant.

## References

  * Shor [1997] P. W. Shor, Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer, SIAM Journal on Computing 26, 1484 (1997).
  * Harrow  _et al._ [2009] A. W. Harrow, A. Hassidim, and S. Lloyd, Quantum algorithm for solving linear systems of equations, Physical Review Letters 103, 150502 (2009).
  * Clader  _et al._ [2013] B. D. Clader, B. C. Jacobs, and C. R. Sprouse, Preconditioned quantum linear system algorithm, Physical Review Letters 110, 250504 (2013).
  * Montanaro [2015] A. Montanaro, Quantum speedup of Monte Carlo methods, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 471, 20150301 (2015).
  * Havlíček  _et al._ [2019] V. Havlíček, A. D. Córcoles, K. Temme, A. W. Harrow, A. Kandala, J. M. Chow, and J. M. Gambetta, Supervised learning with quantum-enhanced feature spaces, Nature 567, 209 (2019).
  * Liu  _et al._ [2021] Y. Liu, S. Arunachalam, and K. Temme, A rigorous and robust quantum speed-up in supervised machine learning, Nature Physics 17, 1013 (2021).
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, A quantum approximate optimization algorithm, arXiv preprint arXiv:1411.4028 (2014).
  * Cao  _et al._ [2019] Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. D. Sawaya, S. Sim, L. Veis, and A. Aspuru-Guzik, Quantum chemistry in the age of quantum computing, Chemical Reviews 119, 10856 (2019).
  * Herman  _et al._ [2023] D. Herman, C. Googin, X. Liu, Y. Sun, A. Galda, I. Safro, M. Pistoia, and Y. Alexeev, Quantum computing for finance, Nature Reviews Physics 5, 450 (2023).
  * Succi  _et al._ [2023] S. Succi, W. Itani, K. Sreenivasan, and R. Steijl, Quantum computing for fluids: Where do we stand?, Europhysics Letters 144, 10001 (2023).
  * Succi  _et al._ [2024] S. Succi, W. Itani, C. Sanavio, K. R. Sreenivasan, and R. Steijl, Ensemble fluid simulations on quantum computers, Computers & Fluids 270, 106148 (2024).
  * Steijl [2019] R. Steijl, Quantum algorithms for fluid simulations, Advances in Quantum Communication and Information , 31 (2019).
  * Moawad  _et al._ [2022] Y. Moawad, W. Vanderbauwhede, and R. Steijl, Investigating hardware acceleration for simulation of CFD quantum circuits, Frontiers in Mechanical Engineering 8, 925637 (2022).
  * Steijl and Barakos [2018] R. Steijl and G. N. Barakos, Parallel evaluation of quantum algorithms for computational fluid dynamics, Computers & Fluids 173, 22 (2018).
  * Gaitan [2020] F. Gaitan, Finding flows of a Navier–Stokes fluid through quantum computing, npj Quantum Information 6, 61 (2020).
  * Bharadwaj and Sreenivasan [2020] S. S. Bharadwaj and K. R. Sreenivasan, Quantum computation of fluid dynamics, arXiv preprint arXiv:2007.09147 (2020).
  * Bharadwaj and Sreenivasan [2023] S. S. Bharadwaj and K. R. Sreenivasan, Hybrid quantum algorithms for flow problems, Proceedings of the National Academy of Sciences 120, e2311014120 (2023).
  * Li  _et al._ [2023] X. Li, X. Yin, N. Wiebe, J. Chun, G. K. Schenter, M. S. Cheung, and J. Mülmenstädt, Potential quantum advantage for simulation of fluid dynamics, arXiv preprint arXiv:2303.16550 (2023).
  * Oz  _et al._ [2022] F. Oz, R. K. Vuppala, K. Kara, and F. Gaitan, Solving Burgers’ equation with quantum computing, Quantum Information Processing 21, 1 (2022).
  * Jaksch  _et al._ [2023] D. Jaksch, P. Givi, A. J. Daley, and T. Rung, Variational quantum algorithms for computational fluid dynamics, AIAA Journal 61, 1885 (2023).
  * Givi  _et al._ [2020] P. Givi, A. J. Daley, D. Mavriplis, and M. Malik, Quantum speedup for aeroscience and engineering, AIAA Journal 58, 3715 (2020).
  * Xu  _et al._ [2018] G. Xu, A. J. Daley, P. Givi, and R. D. Somma, Turbulent mixing simulation via a quantum algorithm, AIAA Journal 56, 687 (2018).
  * Xu  _et al._ [2019] G. Xu, A. J. Daley, P. Givi, and R. D. Somma, Quantum algorithm for the computation of the reactant conversion rate in homogeneous turbulence, Combustion Theory and Modelling 23, 1090 (2019).
  * Sammak  _et al._ [2015] S. Sammak, A. Nouri, N. Ansari, and P. Givi, Quantum computing and its potential for turbulence simulations, in  _Mathematical Modeling of Technological Processes: 8th International Conference, CITech 2015, Almaty, Kazakhstan, September 24-27, 2015, Proceedings 8_ (Springer, 2015) pp. 124–132.
  * Lapworth [2022] L. Lapworth, A hybrid quantum-classical CFD methodology with benchmark HHL solutions, arXiv preprint arXiv:2206.00419 (2022).
  * Demirdjian  _et al._ [2022] R. Demirdjian, D. Gunlycke, C. A. Reynolds, J. D. Doyle, and S. Tafur, Variational quantum solutions to the advection–diffusion equation for applications in fluid dynamics, Quantum Information Processing 21, 322 (2022).
  * Bravo-Prieto  _et al._ [2020] C. Bravo-Prieto, R. LaRose, M. Cerezo, Y. Subasi, L. Cincio, and P. J. Coles, Variational quantum linear solver, arXiv preprint arXiv:1909.05820 (2020).
  * Frisch  _et al._ [1986a] U. Frisch, D. d’Humieres, B. Hasslacher, P. Lallemand, Y. Pomeau, and J. P. Rivet, Lattice gas hydrodynamics in two and three dimensions, in  _Modern Approaches to Large Nonlinear Physical Systems Workshop_ (1986).
  * Frisch  _et al._ [1986b] U. Frisch, B. Hasslacher, and Y. Pomeau, Lattice-gas automata for the Navier–Stokes equation, Physical Review Letters 56, 1505 (1986b).
  * Krüger  _et al._ [2016] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. Viggen,  _The Lattice Boltzmann Method: Principles and Practice_ , Graduate Texts in Physics (Springer International Publishing, 2016).
  * Chen and Doolen [1998] S. Chen and G. D. Doolen, Lattice Boltzmann method for fluid flows, Annual Review of Fluid Mechanics 30, 329 (1998).
  * Yepez [1998] J. Yepez, Lattice-gas quantum computation, International Journal of Modern Physics C 09, 1587 (1998).
  * Yepez [2001] J. Yepez, Quantum lattice-gas model for computational fluid dynamics, Physical Review E 63, 046702 (2001).
  * Yepez [2002] J. Yepez, Quantum lattice-gas model for the Burgers equation, Journal of Statistical Physics 107, 203 (2002).
  * Itani and Succi [2022] W. Itani and S. Succi, Analysis of Carleman linearization of lattice Boltzmann, Fluids 7 (2022).
  * Itani  _et al._ [2023] W. Itani, K. R. Sreenivasan, and S. Succi, Quantum algorithm for lattice Boltzmann (QALB) simulation of incompressible fluids with a nonlinear collision term, arXiv preprint arXiv:2304.05915 (2023).
  * Sanavio and Succi [2023] C. Sanavio and S. Succi, Quantum lattice Boltzmann–Carleman algorithm, arXiv preprint arXiv:2310.17973 (2023).
  * Todorova and Steijl [2020] B. N. Todorova and R. Steijl, Quantum algorithm for the collisionless Boltzmann equation, Journal of Computational Physics 409, 109347 (2020).
  * Budinski [2021a] L. Budinski, Quantum algorithm for the advection–diffusion equation simulated with the lattice Boltzmann method, Quantum Information Processing 20, 57 (2021a).
  * Budinski [2021b] L. Budinski, Quantum algorithm for the Navier–Stokes equations, arXiv preprint arXiv:2103.03804 (2021b).
  * Rothman and Zaleski [1997] D. H. Rothman and S. Zaleski,  _Lattice-gas cellular automata: Simple models of complex hydrodynamics_ (Cambridge University Press, 1997).
  * Rivet and Boon [2001] J.-P. Rivet and J. P. Boon,  _Lattice Gas Hydrodynamics_ (Cambridge University Press, 2001).
  * Ginzburg [2005] I. Ginzburg, Equilibrium-type and link-type lattice Boltzmann models for generic advection and anisotropic-dispersion equation, Advances in Water Resources 28, 1171 (2005).
  * Servan-Camas and Tsai [2009] B. Servan-Camas and F. T.-C. Tsai, Non-negativity and stability analyses of lattice boltzmann method for advection–diffusion equation, Journal of Computational Physics 228, 236 (2009).
  * Shende  _et al._ [2006] V. Shende, S. Bullock, and I. Markov, Synthesis of quantum-logic circuits, IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 25, 1000 (2006).
  * Qiskit Contributors [2023] Qiskit Contributors, Qiskit: An open-source framework for quantum computing (2023).
  * Nielsen and Chuang [2000] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information_ (Cambridge University Press, 2000).
  * Chakraborty  _et al._ [2019] S. Chakraborty, A. Gilyén, and S. Jeffery, The power of block-encoded matrix powers: Improved regression techniques via faster Hamiltonian simulation, arXiv preprint arXiv:1804.01973 (2019).
  * Low and Chuang [2019] G. H. Low and I. L. Chuang, Hamiltonian simulation by qubitization, Quantum 3, 163 (2019).
  * Batchelor [1967] G. K. Batchelor,  _An introduction to fluid dynamics_ (Cambridge University Press, 1967).
  * Meyer [1997] D. A. Meyer, Quantum mechanics of lattice gas automata: One-particle plane waves and potentials, Physical Review E 55, 5261 (1997).


assume 100-$M_{\odot}$ PBHs comprising all the dark matter, but the
qualitative trends are general for all PBH masses. The suppressed amplitude in
the auto-power-spectrum of $\Theta^{(1)}_{\rm inh}$ (purple) compared to
$\Theta^{(1)}_{\rm hom}$ (red) and the poor correlation explains the large
difference in amplitude for the computed power spectra in Sec. IV.2.

##  Appendix F Redshift dependence of the temperature trispectrum induced by
accreting PBHs

In this appendix we inspect the redshift dependence of the temperature
trispectrum from accreting PBHs, by reproducing the forecast analysis of Sec.
V, but artificially imposing that the free-electron fraction perturbation
vanishes outside of redshift bins of size $\Delta z=50$. In a given redshift
bin, we compute the signal-to-noise, $S/N$, assuming a Planck-like experiment
for both the temperature-only trispectrum and power spectrum. That is, for the
trispectrum we compute $(S/N)_{\rm tri}\equiv 1/\sigma_{f_{\rm pbh}}$ from Eq.
(V.2). For the power spectrum we compute the similar forecasted quantity,

$\displaystyle(S/N)_{\rm ps}\equiv\left[\frac{f_{\rm
sky}}{2}\sum_{\ell}(2\ell+1)\left(\frac{C_{\ell}^{(1)}}{C_{\ell}^{\prime}}\right)^{2}\right]^{1/2},$
(175)

where $C_{\ell}^{(1)}=C_{\ell,\rm hom}^{(1)}+C_{\ell,\rm inh}^{(1)}$ is the
total (c.f. Fig. 6) perturbed $TT$ power spectrum due to accreting PBHs
(considering only the “direct” term discussed in Sec. III.4). Note that a
rigorous treatment would properly account for correlations between different
redshift bins, and involve a principal component analysis. Still, the simple
estimation of $S/N$ should give us a reasonable qualitative understanding of
the redshift dependence of the signal.

We compare the two $S/N$ as a function of redshift in Fig. 10 for
100-$M_{\odot}$ PBHs. We see that the temperature trispectrum $S/N$ is rather
sharply peaked around $z\sim 900-1000$, in contrast with the temperature power
spectrum signal, which receives comparable contributions from a broad range of
redshifts $200\lesssim z\lesssim 1200$.

This is consistent with the the following observations. First, by inspecting
the on-the-spot energy deposition limit discussed in Sec. II.3, we found that
the trispectrum is negligibly affected by photon propagation that is more
suppressive at late times. Namely the strongest spatial fluctuations for the
accreting PBHs due to relative velocities occur at a few $10$’s of Mpc scales
as shown in Paper I Fig. 13, but this is not noticeably suppressed until
$z\approx 800$ when inspecting Fig. 2. Second, we find that the trispectrum
constraints converge much more quickly than compared to the power spectrum
when varying the max multipole on the zeroth-order collision term present in
the line-of-sight source. Namely, the trispectrum is unaffected by higher
order multipoles of zeroth-order temperature anisotropy which are induced at
later times. Thirdly and more subtly, the $f_{\rm pbh}$-$M_{\rm pbh}$ powerlaw
dependence is weaker for the trispectrum constraints than it is for the power
spectrum constraints. As found in AK17, the luminosity of a spherically
accreting PBH is proportional to $M^{3}$ at all times when excluding their
radiative efficiency. The radiative efficiency, however, turns out to have an
inverse dependence on black hole mass whose power depends on redshift. This
power converges to zero at late times, and implies that the mass dependence on
$f_{\rm pbh}$ constraints is weaker if the signal receives support from
earlier redshifts. This can be seen directly in Fig. 8 of AK17 where they plot
the mean luminosity as a function of redshift for various $M_{\rm pbh}$.

Figure 10:  Forecasted Planck signal-to-noise ratio for the $TT$-only power
spectrum and $TTTT$ trispectrum induced by accreting PBHs. For ease of
comparison we normalize the curves such that they integrate to unity over
redshift. Each point is computed assuming the perturbed free-electron fraction
is only nonzero in redshift bins of size $\Delta z=50$.

## References

  * Bean and Magueijo (2002) R. Bean and J. Magueijo, Phys. Rev. D 66, 063505 (2002), arXiv:astro-ph/0204486 .
  * Bird  _et al._ (2016) S. Bird _et al._ , Phys. Rev. Lett. 116, 201301 (2016), arXiv:1603.00464 .
  * Carr  _et al._ (2021) B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Reports on Progress in Physics 84, 116902 (2021), arXiv:2002.12778 .
  * Ricotti  _et al._ (2008) M. Ricotti, J. P. Ostriker, and K. J. Mack, ApJ 680, 829 (2008), arXiv:0709.0524 .
  * Ali-Haïmoud and Kamionkowski (2017) Y. Ali-Haïmoud and M. Kamionkowski, Phys. Rev. D 95, 043534 (2017), arXiv:1612.05644 .
  * Poulin  _et al._ (2017) V. Poulin, P. D. Serpico, F. Calore, S. Clesse, and K. Kohri, Phys. Rev. D 96, 083524 (2017), arXiv:1707.04206 .
  * Senatore  _et al._ (2009) L. Senatore, S. Tassev, and M. Zaldarriaga, J. Cosmology Astropart. Phys 2009, 038 (2009), arXiv:0812.3658 .
  * Khatri and Wandelt (2009) R. Khatri and B. D. Wandelt, Phys. Rev. D 79, 023501 (2009), arXiv:0810.4370 .
  * Dvorkin  _et al._ (2013) C. Dvorkin, K. Blum, and M. Zaldarriaga, Phys. Rev. D 87, 103522 (2013), arXiv:1302.4753 .
  * Jensen and Ali-Haïmoud (2021) T. W. Jensen and Y. Ali-Haïmoud, Phys. Rev. D 104, 063534 (2021), arXiv:2106.10266 .
  * Tseliakhovich and Hirata (2010) D. Tseliakhovich and C. Hirata, Phys. Rev. D 82, 083520 (2010), arXiv:1005.2416 .
  * Novosyadlyj (2006) B. Novosyadlyj, MNRAS 370, 1771 (2006), arXiv:astro-ph/0603674 .
  * Huang and Vernizzi (2013) Z. Huang and F. Vernizzi, Phys. Rev. Lett. 110, 101303 (2013), arXiv:1212.3573 .
  * Planck Collaboration (2020a) Planck Collaboration, A&A 641, A9 (2020a), arXiv:1905.05697 .
  * Seljak and Zaldarriaga (1996) U. Seljak and M. Zaldarriaga, ApJ 469, 437 (1996), arXiv:astro-ph/9603033 .
  * Khatri and Wandelt (2010) R. Khatri and B. D. Wandelt, Phys. Rev. D 81, 103518 (2010), arXiv:0903.0871 .
  * Smith  _et al._ (2015) K. M. Smith, L. Senatore, and M. Zaldarriaga, arXiv e-prints (2015), arXiv:1502.00635 .
  * Lee and Ali-Haïmoud (2020) N. Lee and Y. Ali-Haïmoud, Phys. Rev. D 102, 083517 (2020), arXiv:2007.14114 .
  * Ali-Haïmoud and Hirata (2011) Y. Ali-Haïmoud and C. M. Hirata, Phys. Rev. D 83, 043513 (2011), arXiv:1011.3758 .
  * Ali-Haïmoud and Hirata (2010) Y. Ali-Haïmoud and C. M. Hirata, Phys. Rev. D 82, 063521 (2010), arXiv:1006.1355 .
  * Ma and Bertschinger (1995) C.-P. Ma and E. Bertschinger, ApJ 455, 7 (1995), arXiv:astro-ph/9506072 .
  * Blas  _et al._ (2011) D. Blas, J. Lesgourgues, and T. Tram, J. Cosmology Astropart. Phys 2011, 034 (2011), arXiv:1104.2933 .
  * Regan  _et al._ (2010) D. M. Regan, E. P. S. Shellard, and J. R. Fergusson, Phys. Rev. D 82, 023520 (2010), arXiv:1004.2915 .
  * Komatsu  _et al._ (2005) E. Komatsu, D. N. Spergel, and B. D. Wandelt, ApJ 634, 14 (2005), arXiv:astro-ph/0305189 .
  * Planck Collaboration (2020b) Planck Collaboration, A&A 641, A5 (2020b), arXiv:1907.12875 .
  * Planck Collaboration (2020c) Planck Collaboration, A&A 641, A6 (2020c), arXiv:1807.06209 .
  * Meerburg  _et al._ (2016) P. D. Meerburg, J. Meyers, A. van Engelen, and Y. Ali-Haïmoud, Phys. Rev. D 93, 123511 (2016), arXiv:1603.02243 .
  * Inman and Ali-Haïmoud (2019) D. Inman and Y. Ali-Haïmoud, Phys. Rev. D 100, 083528 (2019), arXiv:1907.08129 .
  * Varshalovich  _et al._ (1988) D. A. Varshalovich, A. N. Moskalev, and V. K. Khersonskii,  _Quantum Theory of Angular Momentum_ (WORLD SCIENTIFIC, 1988).


11institutetext: The Robotics Institute, Carnegie Mellon University
11email: {ktracy<EMAIL_ADDRESS>

# On the Differentiability of the
Primal-Dual Interior-Point Method

Kevin Tracy  Zachary Manchester

###### Abstract

Primal-Dual Interior-Point methods are capable of solving constrained convex
optimization problems to tight tolerances in a fast and robust manner. The
derivatives of the primal-dual solution with respect to the problem matrices
can be computed using the implicit function theorem, enabling efficient
differentiation of these optimizers for a fraction of the cost of the total
solution time. In the presence of active inequality constraints, this
technique is only capable of providing discontinuous subgradients that present
a challenge to algorithms that rely on the smoothness of these derivatives.
This paper presents a technique for relaxing primal-dual solutions with a
logarithmic barrier to provide smooth derivatives near active inequality
constraints, with the ability to specify a uniform and consistent amount of
smoothing. We pair this with an efficient primal-dual interior-point algorithm
for solving an always-feasible $\ell_{1}$-penalized variant of a convex
quadratic program, eliminating the issues surrounding learning potentially
infeasible problems. This parallelizable and smoothly differentiable solver is
demonstrated on a range of robotics tasks where smoothing is important. An
open source implementation in JAX is available at www.github.com/kevin-
tracy/qpax.

###### Keywords:

Differentiable Optimization Interior-Point Methods Convex Optimization.

##  1 Introduction

Convex optimization has seen widespread use in modern robotics, where the
guarantees of global optimality and polynomial time complexity have enabled
algorithms that span control [15, 6], state estimation [24, 25], actuator
allocation [22, 14], collision detection [10, 21], and simulation [5, 19].

For years, practitioners have exploited domain-specific knowledge to craft
convex optimization problems that enjoy fast and reliable convergence for both
offline and online use. In the era of data-driven robotics, differentiable
optimization has enabled automatic tuning or “learning” of optimization
problems directly from data. Although this may not replace domain-specific
knowledge, the ability to build complex tuneable functions with embedded
optimization problems is well suited for a variety of tasks.

The sensitivity analysis of linear programs has been studied for decades [7],
but recent advances in differentiable convex optimization and parametrized
convex layers in deep networks are gaining traction [4, 1, 2]. This has
resulted in convex modeling tools such as CVXPY layers that enable easy
incorporation of differentiable convex optimization into common workflows [8].

There are two main issues preventing the widespread adoption of differentiable
optimization in robotics as it exists today: the subgradient problem and the
infeasibility problem. For active inequality constraints in optimization
problems, the derivatives propagated through these solvers are restricted to
subgradients when there is no uniquely defined gradient. This results in
routines getting “stuck” near these inequalities due to their nonsmooth
nature. The second issue arises when infeasible problem instances are created
during auto-tuning. While infeasibility detection is commonplace in convex
optimization, detecting an infeasible problem does not address the root cause
of the infeasibility, nor provide informative gradient information to
encourage feasibility.

In this paper, we propose solutions to both of these problems. Our
contributions include:

  * •

A rigorous method for returning unique and smoothed derivatives of convex
optimization problems through the use of a relaxed logarithmic barrier

  * •

A primal-dual interior-point algorithm for solving an always-feasible
“elastic” quadratic program with minimal computational overhead

Together, these two advances are presented in an open source software package
written in JAX and demonstrated on the relevant robotics tasks of contact
dynamics and collision detection.

##  2 Background

This section introduces the notation for a standard form of the convex
Quadratic Program (QP), common solution methods, and extensions to
differentiable optimization.

###  2.1 Quadratic Programming

In this paper, we focus on the convex QP as a fundamental problem
specification. Modern algorithms such as the Primal-Dual Interior-Point (PDIP)
method solve these problems globally in a fast and efficient manner [16]. A
standard form quadratic program and its equivalent with a slack variable are
as follows:

x 12x^TQx + q^Tx Ax=bGx≤h,  |  $\Rightarrow$ x, s 12x^TQx + q^Tx Ax=bGx + s= hs≥0,
---|---

with a primal variable $x\in\mathbf{R}^{n}$, cost terms
$Q\in\mathbf{S}_{+}^{n}$ and $q\in\mathbf{R}^{n}$, equality constraints
described with $A\in\mathbf{R}^{m\times n}$ and $b\in\mathbf{R}^{m}$, and
inequality constraints with $G\in\mathbf{R}^{p\times n}$ and
$h\in\mathbf{R}^{p}$. Dual variables are introduced to enforce the
constraints, with $y\in\mathbf{R}^{m}$ associated with the equality
constraint, and $z\in\mathbf{R}^{p}$ with the inequality constraints [7]. The
slack variable $s\in\mathbf{R}^{p}$ is introduced for algorithmic simplicity.
The Lagrangian for this problem is then

$\displaystyle\mathcal{L}(x,s,z,y)$
$\displaystyle=\frac{1}{2}x^{T}Qx+q^{T}x+y^{T}(Ax-b)+z^{T}(Gx-h),$ (1)

resulting in the following KKT conditions for optimality:

$\displaystyle Qx+q+G^{T}z+A^{T}y$ $\displaystyle=0,$ (2) $\displaystyle
z\odot s$ $\displaystyle=0,$ (3) $\displaystyle Gx+s$ $\displaystyle=h,$ (4)
$\displaystyle Ax$ $\displaystyle=b,$ (5) $\displaystyle s$ $\displaystyle\geq
0,$ (6) $\displaystyle z$ $\displaystyle\geq 0,$ (7)

where $\odot$ denotes elementwise multiplication. A primal-dual solution
$(x^{*},s^{*},y^{*},z^{*})$ is globally optimal if it satisfies (2)-(7).

###  2.2 Primal-Dual Interior-Point Methods

PDIP methods solve (2.1) by treating a modified version of the system of
equations in (2)-(5) as a root finding problem, and then using Newton’s method
to find a solution while restricting $(s,z)>0$. As a result, the majority of
the computation time is spent solving linear systems of the following form:

$\displaystyle\begin{bmatrix}Q&0&G^{T}&A^{T}\\\ 0&D(z)&D(s)&0\\\ G&I&0&0\\\
A&0&0&0\end{bmatrix}\begin{bmatrix}\Delta x\\\ \Delta s\\\ \Delta z\\\ \Delta
y\end{bmatrix}=\begin{bmatrix}u_{1}\\\ u_{2}\\\ u_{3}\\\ u_{4}\end{bmatrix},$
(8)

where $D(\cdot)$ denotes the diagonal matrix constructor from a vector. Using
block reduction techniques, the linear system in (8) can be efficiently solved
with Alg. (1). This technique for solving linear systems of this form is
useful for both the solving and differentiation of the PDIP method.

After the step directions are computed with the linear system, a linesearch is
used to ensure the nonnegativity of $(s,z)$ . For an arbitrary variable $v$
and step direction $\Delta v$, a linesearch that solves for the largest
$\alpha\leq 1$ that keeps $v+\Delta v\geq 0$ is solved in closed form:

$\displaystyle\operatorname{linesearch}(v,\Delta v)$
$\displaystyle=\min\bigg{(}1,\,\min_{i:\Delta v_{i}<0}-\frac{v_{i}}{\Delta
v_{i}}\bigg{)}.$ (9)

This linesearch is performed for both $s$ and $z$, and the step length is
simply the minimum of the two. For a complete algorithmic specification and
implementation details, we refer the reader to [16] for a PDIP method for
which Alg. (1) can be used to solve the linear systems.

Algorithm 1 PDIP Linear System Solver

1:function $\operatorname{solve\\_kkt}(u_{1},u_{2},u_{3},u_{4})$
$\triangleright$ assume access to $Q,G,A,s,z$

2:$P\leftarrow D(s\oslash z)$

3:$H\leftarrow Q+G^{T}P^{-1}G$ $\triangleright$ cacheable Cholesky
decomposition

4:$F\leftarrow AH^{-1}A^{T}$ $\triangleright$ cacheable Cholesky decomposition

5:$r_{2}\leftarrow u_{3}-u_{2}\oslash z$

6:$p_{1}\leftarrow u_{1}+G^{T}P^{-1}r_{2}$

7:$\Delta y\leftarrow F^{-1}(AH^{-1}p_{1}-u_{4})$

8:$\Delta x\leftarrow H^{-1}(p_{1}-A^{T}\Delta y)$

9:$\Delta s\leftarrow u_{3}-G\Delta x$

10:$\Delta z\leftarrow(u_{2}-z\odot\Delta s)\oslash s$

11:return $(\Delta x,\Delta s,\Delta z,\Delta y)$

###  2.3 Differentiable Optimization

Modern automatic differentiation tools have made forming derivatives through
complex functions easier than ever. When one of these functions includes an
iterative routine, unrolling the iterations into a sequential computational
graph and proceeding with differentiation is not always possible. The first
concern is nonsmooth operations or logical branching in the iterations that
can stop the “flow” of the derivatives through the routine, and the second
concern is the rapidly decaying numerical precision inherent in
differentiating an iterative process.

To address both of these shortcomings, iterative routines such as numerical
optimizers are differentiated with methods that do not require propagating
derivatives through the iterations themselves. Instead, the optimization
problem is solved as normal, and the derivatives are then constructed directly
from the solution to the problem.

####  2.3.1 Implicit Function Theorem

When an iterative routine can be interpreted as finding a root or an
equilibrium point to an implicit function, the implicit function theorem can
be used to form the derivatives of interest. Given variables
$w\in\mathbf{R}^{a}$ and parameters $\theta\in\mathbf{R}^{b}$, an implicit
function is defined as

$\displaystyle r(w^{*},\theta)=0,$ (10)

at an equilibrium point $w^{*}$. This implicit function can be linearized
about this point resulting in the following first-order Taylor series:

$\displaystyle\frac{\partial{r}}{\partial{w}}\delta
w+\frac{\partial{r}}{\partial{\theta}}\delta\theta=0,$ (11)

which can simply be re-arranged to solve for

$\displaystyle\frac{\partial{w}}{\partial{\theta}}=-\bigg{(}\frac{\partial{r}}{\partial{w}}\bigg{)}^{-1}\frac{\partial{r}}{\partial{\theta}}.$
(12)

The implicit function theorem enables the differentiation of routines that
solve for equilibrium points without the need to unroll and differentiate the
iterations themselves.

By treating the KKT conditions from (2)-(5) as a residual function of the
primal-dual solution $(x,s,z,y)$ and problem parameters $\theta$, the implicit
function theorem is used to form the derivatives of the optimizer without
unrolling the iterations. The linear system from (12) applied to this residual
function is of the same form as that in Alg. (1), enabling fast and easy
computation of derivatives using a framework already available in the solver.

####  2.3.2 Efficient Computation of Gradients

Using the implicit function theorem to compute the Jacobians of the primal-
dual solution with respect to the problem parameters results in the need to
directly form potentially large Jacobians. In many cases, it is not these
specific Jacobians that are of interest, but rather the left matrix-vector
product with a backward pass vector.

For a loss function of the primal variable $\ell(x)$ that takes as input the
optimal primal variable from (2.1), reverse-mode automatic differentiation
will have constructed $\nabla_{x}\ell$ by the time it comes to the QP solver
in the backward pass of the computational graph. Instead of forming the
Jacobians of the primal variable with respect to each of the problem matrices
directly, it is instead desirable to form the left matrix-vector products
$\frac{\partial{\ell}}{\partial{x}}\frac{\partial{x}}{\partial{\square}}$,
where $\square$ simply denotes any of the problem matrices [3].

As shown in Alg. (2), the gradients of this loss function
$\nabla_{\square}\ell$ with respect to the problem parameters of the QP can be
computed by once again utilizing our interior-point linear system solver.

1:function
$\operatorname{compute\\_qp\\_grads}(Q,q,A,b,G,h,x,s,z,y,\nabla_{x}\ell)$

2:$dx,ds,d\tilde{z},dy\leftarrow\operatorname{solve\\_kkt}(-\nabla_{x}\ell,0,0,0)$
$\triangleright$ compute differentials with kkt system

3:$dz=d\tilde{z}\oslash z$

4:$\nabla_{Q}\ell\leftarrow\big{(}(dx)x^{T}+x(dx)^{T})\big{)}/2$

5:$\nabla_{q}\ell\leftarrow dx$

6:$\nabla_{A}\ell\leftarrow(dy)x^{T}+y(dx)^{T}$

7:$\nabla_{b}\ell\leftarrow-dy$

8:$\nabla_{G}\ell\leftarrow z\odot\big{(}(dz)x^{T}+z(dx)^{T}\big{)}$

9:$\nabla_{h}\ell\leftarrow-z\odot dz$

10:return
$\nabla_{Q}\ell,\nabla_{q}\ell,\nabla_{A}\ell,\nabla_{b}\ell,\nabla_{G}\ell,\nabla_{h}\ell$

Algorithm 2 Computing Gradients Through a QP

##  3 Logarithmic Barrier Smoothing

Figure 1: A sharp corner in a square as smoothed with the logarithmic barrier
at varying central path parameters $\kappa$. As $\kappa\rightarrow 0$, the
corner becomes more pronounced until it assumes a true $90^{\circ}$ corner at
$\kappa=0$. The logarithmic barrier effectively smooths out any sharp corners
of the feasible set enabling smooth differentiation in the presence of such
discontinuities.

Primal interior-point methods work by replacing constraints of the form
$h(x)\leq 0$ with a logarithmic barrier penalty in the form of
$\phi(x)=-\kappa\sum\log(-h(x))$, where $\kappa\in\mathbf{R}^{+}$ is referred
to as the central path parameter. This barrier function is a smooth
approximation of the indicator function, where feasible values of $x$ result
in no penalty, and infeasibility results in an infinite penalty [7]. By
solving a sequence of unconstrained problems as $\kappa\rightarrow 0$, the
logarithmic barrier becomes a closer and closer approximation of the indicator
function until acceptable convergence is achieved. This barrier function is
only defined for feasible values of $x$, hence the name “interior-point”.

For a quadratic program in standard form (2.1), the optimality conditions for
a barrier subproblem given a central path parameter $\kappa$ are almost
identical to the original KKT conditions (2)-(7) with the exception of the
complementarity condition (3) replaced with $z\odot s-\kappa=0$. This relaxed
complementarity condition allows for a certain amount of smoothing over the
feasible set, where larger values of $\kappa$ have a stronger smoothing
effect. The optimality conditions for this barrier subproblem are referred to
as the perturbed or relaxed KKT conditions.

###  3.1 Relaxing Primal-Dual Solutions

Algorithm 3 Relaxing a Quadratic Program

1:function $\operatorname{relax\\_qp}(Q,q,A,b,G,h,x,s,z,y,\kappa)$

2:for $i\leftarrow 1:\texttt{max\\_iters}$ do

3: /* evaluate relaxed KKT conditions and check convergence*/

4: $r_{1}\leftarrow Qx+q+A^{T}y+G^{T}z$

5: $r_{2}\leftarrow s\odot z-\kappa$ $\triangleright$ relaxed complementarity

6: $r_{3}\leftarrow Gx+s-h$

7: $r_{4}\leftarrow Ax-b$

8: if $\|(r_{1},r_{2},r_{3},r_{4},r_{5},r_{6})\|_{\infty}<\texttt{tol}$ then

9: return: $x,t,s_{1},s_{2},z_{1},z_{2}$

10: end if

11:

12: /* calculate and take Newton step */

13: $\Delta x,\Delta s,\Delta z,\Delta
y\leftarrow\operatorname{solve\\_kkt}(-r_{1},-r_{2},-r_{3},-r_{4})$
$\triangleright$ Alg. (1)

14: $\alpha\leftarrow 0.98\min\big{(}\operatorname{linesearch}(s,\Delta
s),\,\operatorname{linesearch}(z,\Delta z)\big{)}$ $\triangleright$ (9)

15: $(x,s,z,y)\leftarrow(x,s,z,y)+\alpha(\Delta x,\Delta s,\Delta z,\Delta y)$

16:end for

When differentiating a quadratic program near an active inequality constraint,
the implicit function theorem can produce subgradients that are potentially
uninformative [1]. In order to ensure smooth and continuous gradients in the
presence of sharp corners in the feasible set, proposed here is a relaxation
method that exploits the smoothing of the logarithmic barrier to effectively
round out any nonsmoothness. Specifically, the idea is to take a primal-dual
solution that is optimal to some low $\kappa_{\text{low}}$, and relax it to a
specified $\kappa_{\text{high}}$.

In PDIP methods, much care is taken to solve the perturbed KKT conditions for
decreasing $\kappa\rightarrow 0$. The most common strategy in PDIP methods is
a Mehrotra predictor-corrector method that adaptively updates the target
$\kappa$ until it is below the convergence criteria [17]. While solving this
system for a sequence where $\kappa\rightarrow 0$ is challenging, going the
other direction from $\kappa_{\text{low}}\rightarrow\kappa_{\text{high}}$ is
actually quite trivial. For this case, standard Newton steps on the perturbed
KKT conditions with a linesearch to ensure $(s,z)>0$ is able to converge to
$\kappa_{\text{high}}$ in only a few steps. This algorithm is shown in Alg.
(3), where once again the linear system solver from Alg. (1) is used.

Once the primal-dual solution to the quadratic program has been relaxed, Alg.
(2) is used to calculate gradients that benefit from the logarithmic barrier
smoothing. The full sequence for the solving, relaxation, and differentiation
of a quadratic problem is as follows:

  1. 1.

Solve the quadratic program to a specified tolerance and return the solution

  2. 2.

Relax the primal-dual solution to a target $\kappa$

  3. 3.

Form the derivatives of interest at the relaxed primal-dual solution

This means that the solver itself can return high quality solutions to tight
tolerances while still returning smooth gradients evaluated at the relaxed
solution. In modern automatic differentiation frameworks, this sequence is
written into a custom forward and backward pass through the function with
limited overhead.

The choice of a target $\kappa$ is left to the specifics of the problem at
hand. In some scenarios, only a little bit of smoothing is required, making a
lower $\kappa$ appropriate. For very sharp corners in the feasible set (like
the tip of a triangle), a larger value of $\kappa$ can be used to provide even
more smoothing. In either case, computing the derivatives of the solver with
this technique allows both the tolerance of the solver and the relaxed
$\kappa$ to be specified exactly and independently.

##  4 Elastic Quadratic Program Solver

While the convexity of a quadratic program guarantees a globally optimal
solution when one is available, there is generally no guarantee of feasibilty.
This infeasibility occurs when the quadratic program has a set of constraints
that are impossible to satisfy, something that can easily happen in practice
if the constraint matrices are learned. In the event of an infeasible problem,
standard PDIP methods are unable to return a useful solution.

Traditionally, infeasibility in convex optimization problems can be handled
with a homogenous self-dual embedding that allows for the computation of a
certificate of infeasibility [9, 23, 20]. This ensures that, even in the event
of an infeasible problem, the solver can “gracefully” fail in a way that
simply returns certification that there are no solutions that satisfy the
constraints. While this approach is useful for determining if a given problem
is feasible, ideally we would set up problems such that infeasibility is not
possible.

Given a convex quadratic program in a standard inequality-only form, x 12x^TQx
+ q^Tx Gx≤h, it is possible there is no $x$ that satisfies $Gx\leq h$. In
order to convert the optimization problem in (4) into one in which there is
always a solution, the hard constraints are converted into penalties. An
$\ell_{1}$-penalty on the constraint violation is chosen because the
$\ell_{1}$-norm encourages sparsity in the constraint violation, translating
into a penalty that encourages the solver to satisfy as many of the
constraints as possible. This is a common technique in nonlinear programming
[18] for the handling of infeasible subproblems. This “elastic” mode, as it is
referred to in the solver SNOPT [12], is a highly effective method to
guarantee that a problem always has a solution without sacrificing the quality
and utility of the solution. The problem from (4) is converted into its
elastic form as follows: x 12x^TQx + q^Tx + ∥ρ⊙max(0, Gx-h)∥_1, such that
feasible values of $x$ do not contribute to the cost function but
infeasibility is penalized with $\rho\in\mathbf{R}^{p}$. While (4) is an
unconstrained convex optimization problem, it is nonsmooth and does not pair
well with PDIP methods. We therefore reformulate it as: x, t 12x^TQx + q^Tx +
ρ^T t Gx - h≤tt≥0, where $t\in\mathbf{R}^{p}$ is a slack variable containing
the constraint violation, and a simple linear cost term is used to recover the
$\ell_{1}$-penalty from (4).

Solving the elastic quadratic program in (4) with a PDIP method normally comes
at an added computational expense since we have increased the number of
primal-dual variables from $n+2p$ to $n+5p$. To avoid the cubic complexity in
the increase primal-dual dimension, we introduce a custom algorithm for
solving these problems that exploits the sparsity of the constraints such that
the time to solve the elastic version of the problem is only a slight (5–20%)
increase compared to the time to solve the original problem.

As before, the PDIP method for solving (4) is dominated by the factorization
and solving of linear systems in the following form:

$\displaystyle\begin{bmatrix}Q&0&0&0&0&G^{T}\\\ 0&0&0&0&-I&-I\\\
0&0&Z_{1}&0&S_{1}&0\\\ 0&0&0&Z_{2}&0&S_{2}\\\ 0&-I&I&0&0&0\\\
G&-I&0&I&0&0\end{bmatrix}\begin{bmatrix}\Delta x\\\ \Delta t\\\ \Delta
s_{1}\\\ \Delta s_{2}\\\ \Delta z_{1}\\\ \Delta z_{2}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}r_{1}\\\ r_{2}\\\ r_{3}\\\ r_{4}\\\ r_{5}\\\
r_{6}\end{bmatrix}.$ (13)

As shown in Alg. (5), this linear system can be solved with block-wise
elimination where the only matrix factorization required is that of a positive
definite matrix the size of the primal variable $x$. This routine is used in
the full PDIP algorithm for the elastic QP as shown in Alg. (4).

The elastic QP is fully differentiable in the same way the original QP is. To
do this, Alg. (6) is used to take a primal-dual solution and relax it to a
specified $\kappa$, and Alg. (5) for solving the linear systems is re-used for
fast and efficient relaxation of the elastic problem. From this, the gradients
of a downstream loss function with respect to the problem parameters can be
constructed with Alg. (7), again using the same linear system solver.

The solving, relaxation, and differentiation of the elastic mode QP are all
only a slight increase in computational complexity compared to the original
QP, with the benefit of guaranteed feasibility. This enables the inclusion of
always-feasible quadratic programming in learned pipelines where feasibility
cannot be guaranteed by construction.

Algorithm 4 PDIP Method for Elastic Quadratic Programs

1:function $\operatorname{solve\\_qp\\_elastic}(Q,q,G,h,\rho)$

2:$x,s_{1},s_{2},z_{1},z_{2}\leftarrow\operatorname{initialize}(Q,q,G,h,\rho)$
$\triangleright$ 4.1

3:for $i\leftarrow 1:\texttt{max\\_iters}$ do

4: /* evaluate KKT conditions and check convergence*/

5: $r_{1}\leftarrow Qx+q+G^{T}z_{2}$

6: $r_{2}\leftarrow-z_{1}-z_{2}+\rho$

7: $r_{3}\leftarrow s_{1}\odot z_{1}$

8: $r_{4}\leftarrow s_{2}\odot z_{2}$

9: $r_{5}\leftarrow-t+s_{1}$

10: $r_{6}\leftarrow Gx-t+s_{2}-h$

11: if $\|(r_{1},r_{2},r_{3},r_{4},r_{5},r_{6})\|_{\infty}<\texttt{tol}$ then

12: return: $x,t,s_{1},s_{2},z_{1},z_{2}$

13: end if

14:

15: /* calculate affine step direction Alg. (5)*/

16: $\Delta x^{\text{a}},\Delta t^{\text{a}},\Delta s_{1}^{\text{a}},\Delta
s_{2}^{\text{a}},\Delta z_{1}^{\text{a}},\Delta
z_{2}^{\text{a}}\leftarrow\operatorname{elastic\\_kkt}(-r_{1},-r_{2},-r_{3},-r_{4},-r_{5},-r_{6})$

17: $\Delta s^{\text{a}},\,\Delta z^{\text{a}}\leftarrow(\Delta
s_{1}^{\text{a}},\Delta s_{2}^{\text{a}}),\,(\Delta z_{1}^{\text{a}},\Delta
z_{2}^{\text{a}})$

18: $s,\,z\leftarrow(s_{1},s_{2}),\,(z_{1},z_{2})$

19: /* calculate centering-plus-corrector step direction */

20: $\alpha^{\text{a}}=\min(\operatorname{linesearch}(s,\Delta
s^{\text{a}}),\,\operatorname{linesearch}(z,\Delta z^{\text{a}}))$
$\triangleright$ (9)

21: $\mu\leftarrow s^{T}z/\operatorname{len}(s)$

22: $\sigma\leftarrow[(s+\alpha^{\text{a}}\Delta
s^{\text{a}})^{T}(z+\alpha^{\text{a}}\Delta z^{\text{a}})/(s^{T}z)]^{3}$

23: $r_{3}\leftarrow r_{3}-\sigma\mu\mathbf{1}+\Delta
s_{1}^{\text{a}}\odot\Delta z_{1}^{\text{a}}$

24: $r_{4}\leftarrow r_{4}-\sigma\mu\mathbf{1}+\Delta
s_{2}^{\text{a}}\odot\Delta z_{2}^{\text{a}}$

25: $\Delta x,\Delta t,\Delta s_{1},\Delta s_{2},\Delta z_{1},\Delta
z_{2}\leftarrow\operatorname{elastic\\_kkt}(-r_{1},-r_{2},-r_{3},-r_{4},-r_{5},-r_{6})$

26:

27: /* update with linesearch */

28: $\alpha\leftarrow 0.98\min(\operatorname{linesearch}(s,\Delta
s),\,\operatorname{linesearch}(z,\Delta z))$

29:
$(x,t,s_{1},s_{2},z_{1},z_{2})\leftarrow(x,t,s_{1},s_{2},z_{1},z_{2})+\alpha(\Delta
x,\Delta t,\Delta s_{1},\Delta s_{2},\Delta z_{1},\Delta z_{2})$

30:end for

1:function $\operatorname{elastic\\_kkt}(r_{1},r_{2},r_{3},r_{4},r_{5},r_{6})$
$\triangleright$ assume access to $Q,G,A,s_{1},s_{2},z_{1},z_{2}$

2:$w_{1}\leftarrow r_{3}\oslash z_{1}$

3:$w_{2}\leftarrow r_{4}\oslash z_{2}$

4:$p_{1}\leftarrow r_{5}-r_{6}+w_{2}-w_{1}-(s_{1}\odot r_{2})\oslash z_{1}$

5:$A_{3}\leftarrow\operatorname{diag}(a_{1}+a_{2})$

6:$\Delta x\leftarrow(Q+G^{T}A_{3}^{-1}G)^{-1}(r_{1}-G^{T}A_{3}^{-1}p_{1})$
$\triangleright$ Cholesky factorization (cacheable)

7:$\Delta z_{2}\leftarrow A_{3}^{-1}(p_{1}+G\Delta x)$

8:$\Delta z_{1}\leftarrow-r2-\Delta z_{2}$

9:$\Delta s_{1}\leftarrow(r_{3}-s_{1}\odot\Delta z_{1})/z_{1}$

10:$\Delta s_{2}\leftarrow(r_{4}-s_{2}\odot\Delta z_{2})/z_{2}$

11:$\Delta t\leftarrow\Delta s_{1}-r_{5}$

12:return $\Delta x,\,\Delta s_{1},\,\Delta s_{2},\,\Delta z_{1},\,\Delta
z_{2}$

Algorithm 5 Solve Elastic KKT Linear System

###  4.1 Elastic Initialization

In order to initialize the PDIP method shown in Alg. (4), the only requirement
is that $s,z>0$. In practice, a more advanced initialization technique can
both reduce the number of iterations required for convergence and dramatically
improve the robustness of the solver. The initialization from [23] and [16] is
adapted for use in the elastic case.

First, the solution to the following problem is computed analytically: x, t,
s_1, s_2 12x^TQx + q^Tx + 12(s_1^Ts_1 + s_2^Ts_2) + ρ^T t s_1 - t=0Gx - t +
s_2= h, where the primal and dual solutions are the solution to the linear
system,

$\displaystyle\begin{bmatrix}Q&0&0&G^{T}\\\ 0&0&-I&-I\\\ 0&-I&-I&0\\\
G&-I&0&-I\end{bmatrix}\begin{bmatrix}\Delta x\\\ \Delta t\\\ \Delta z_{1}\\\
\Delta z_{2}\end{bmatrix}$ $\displaystyle=\begin{bmatrix}-q\\\ \rho\\\ 0\\\
h\end{bmatrix},$ (14)

which can be solved with a dense block reduction:

$\displaystyle x$
$\displaystyle=(Q+\frac{1}{2}G^{T}G)^{-1}\big{(}-q-\frac{1}{2}G^{T}(\rho-h)\big{)},$
(15) $\displaystyle z_{2}$ $\displaystyle=\frac{1}{2}(Gx+\rho-h),$ (16)
$\displaystyle z_{1}$ $\displaystyle=\rho-z_{2},$ (17) $\displaystyle t$
$\displaystyle=-z_{1}.$ (18)

This method only requires one positive definite matrix the size of $x$. From
here, we stack $z=(z_{1},z_{2})$ and initialize $s=(s_{1},s_{2})$ with

$\displaystyle\alpha_{p}$ $\displaystyle=-\min(-z),$ (19) $\displaystyle s$
$\displaystyle=\begin{cases}-z,&\alpha_{p}<0\\\
-z+(1+\alpha_{p}\mathbf{1}),&\alpha_{p}\geq 0\end{cases}.$ (20)

The dual variable for the inequality constraint is then initialized with

$\displaystyle\alpha_{d}$ $\displaystyle=-\min(-z),$ (21) $\displaystyle z$
$\displaystyle=\begin{cases}z+(1+\alpha_{d}\mathbf{1}),&\alpha_{d}\geq 0\\\
z,&\alpha_{d}<0\end{cases},$ (22)

finishing the initialization of the primal and dual variables.

Algorithm 6 Relaxing an Elastic Quadratic Program

1:function
$\operatorname{relax\\_qp\\_elastic}(Q,q,A,b,G,h,\rho,x,t,s_{1},s_{2},z_{1},z_{2},\kappa)$

2:for $i\leftarrow 1:\texttt{max\\_iters}$ do

3: /* evaluate KKT conditions and check convergence*/

4: $r_{1}\leftarrow Qx+q+G^{T}z_{2}$

5: $r_{2}\leftarrow-z_{1}-z_{2}+\rho$

6: $r_{3}\leftarrow s_{1}\odot z_{1}-\kappa$ $\triangleright$ relaxed
complementarity

7: $r_{4}\leftarrow s_{2}\odot z_{2}-\kappa$ $\triangleright$ relaxed
complementarity

8: $r_{5}\leftarrow-t+s_{1}$

9: $r_{6}\leftarrow Gx-t+s_{2}-h$

10: if $\|(r_{1},r_{2},r_{3},r_{4},r_{5},r_{6})\|_{\infty}<\texttt{tol}$ then

11: return: $x,t,s_{1},s_{2},z_{1},z_{2}$

12: end if

13:

14: /* calculate and take Newton step*/

15: $\Delta x,\Delta t,\Delta s_{1},\Delta s_{2},\Delta z_{1},\Delta
z_{2}\leftarrow\operatorname{elastic\\_kkt}(-r_{1},-r_{2},-r_{3},-r_{4},-r_{5},-r_{6})$

16: $\Delta s,\,\Delta z\leftarrow(\Delta s_{1},\Delta s_{2}),\,(\Delta
z_{1},\Delta z_{2})$

17: $s,\,z\leftarrow(s_{1},s_{2}),\,(z_{1},z_{2})$

18: $\alpha^{\text{a}}=0.98\min(\operatorname{linesearch}(s,\Delta
s),\,\operatorname{linesearch}(z,\Delta z))$ $\triangleright$ (9)

19:
$(x,t,s_{1},s_{2},z_{1},z_{2})\leftarrow(x,t,s_{1},s_{2},z_{1},z_{2})+\alpha(\Delta
x,\Delta t,\Delta s_{1},\Delta s_{2},\Delta z_{1},\Delta z_{2})$

20:end for

1:function
$\operatorname{compute\\_elastic\\_qp\\_grads}(Q,q,G,h,x,t,s_{1},s_{2},z_{1},z_{2},\nabla_{x}\ell)$

2:/* compute differentials with same linear system */

3:$dx,dt,ds_{1},ds_{2},d\tilde{z}_{1},d\tilde{z}_{2}\leftarrow\operatorname{elastic\\_kkt}(-\nabla_{x}\ell,0,0,0,0,0)$

4:$p,s,z\leftarrow(x,t),(s_{1},s_{2}),(z_{1},z_{2})$

5:$dp,ds,dz\leftarrow(dx,dt),(ds_{1},ds_{2}),(d\tilde{z}_{1},d\tilde{z}_{2})\oslash(z_{1},z_{2})$

6:/* create indices (one-based) and form gradients from differentials */

7:$ix\leftarrow 1:\operatorname{len}(q)$

8:$is\leftarrow(\operatorname{len}(q)+1):(\operatorname{len}(q)+\operatorname{len}(h))$

9:$\nabla_{Q}\ell\leftarrow\big{[}\big{(}(dp)p^{T}+p(dp)^{T})\big{)}/2\big{]}_{ix,ix}$

10:$\nabla_{q}\ell\leftarrow[dp]_{ix}$

11:$\nabla_{G}\ell\leftarrow\big{[}z\odot\big{(}(dz)p^{T}+z(dp)^{T}\big{)}\big{]}_{is,ix}$

12:$\nabla_{h}\ell\leftarrow[-z\odot dz]_{is}$

13:return
$\nabla_{Q}\ell,\nabla_{q}\ell,\nabla_{A}\ell,\nabla_{b}\ell,\nabla_{G}\ell,\nabla_{h}\ell$

Algorithm 7 Computing Gradients Through an Elastic QP

##  5 Numerical Experiments

The utility of the proposed approaches are demonstrated on the common
optimization-based robotics tasks of contact mechanics simulation and
collision detection. In each of these tasks, a QP makes up a core part of the
algorithm and smooth differentiation through the inherent nonsmoothness proves
useful.

###  5.1 Contact Mechanics

Figure 2: Contact dynamics for a two-dimensional block as modeled with a
quadratic program. When a horizontal force $f_{x}$ is applied to the block, it
must overcome the friction with the ground before it moves, and when a
vertical force $f_{y}$ is applied, it must overcome gravity. Despite these
discontinuities in the dynamics, the relaxed gradients from the differentiable
quadratic program solver are able to provide smooth and continuous derivative
information before and after the block begins to move. Figure 3: Using the
gradients from the block pushing example in Fig. 2, a vertical force is
optimized to accelerate the block to a set value from multiple different
initial force values. When the force is below the threshold for the block to
move and $\kappa=0$, the gradient is zero and the optimizer fails to make
progress. In the case of $\kappa=0.01$, even before the block moves there is
gradient information that pushes the optimizer to converge on the optimal
force regardless of the initial force value.

For a block at rest on a table with both gravity and friction, the nonsmooth
contact dynamics can be represented as the solution to a convex quadratic
program. For a full treatment of optimization-based dynamics, the reader is
referred to [5] and [13].

In the simplest case, a block is stationary until it is acted upon by a force
that exceeds the static frictional force in the horizontal direction, or the
gravitational force in the vertical direction. Until the applied forces exceed
these two thresholds, the block does not move. This exercise is demonstrated
in Fig. 2, where forces are applied in each of the two directions and the true
derivatives of these dynamics at $\kappa=0$ have a discontinuity as soon as
the block begins to move. When these derivatives are taken with a relaxed
$\kappa>0$, the discontinuous derivative is smoothed, allowing for an
informative gradient about the impending motion of the block before it moves.

This derivative information is used in an optimization routine in Fig. 3,
where an optimizer attempts to solve for an applied force that produces the
desired motion from the block. This optimizer is initialized at multiple
different force values, and uses either exact gradients with $\kappa=0$, or
smoothed gradients with $\kappa=0.01$. In the case of exact gradients, for
initial forces where the block does not move, the gradient is zero and the
solver is unable to find a descent direction. When $\kappa$ is relaxed, the
smooth gradients provide information about the motion of the block even before
the block itself begins to move, allowing the optimizer to converge on the
true solution for each initial force. This is a simple yet expressive
demonstration of the impact these smooth gradients have on optimization
routines in the presence of discontinuous subgradients.

###  5.2 Collision Detection

$\kappa=0$$\kappa=0.01$

Figure 4: Contact normal vectors from an optimization-based differentiable
collision detection routine with and without relaxed differentiation. With no
relaxation ($\kappa=0$), the direction of the contact normal switches
immediately as the closest point moves from one face to another. When the
relaxed gradients are used ($\kappa=0.01$), the contact normal smoothly
transitions between the faces.

Collision detection between convex shapes can be formulated as a convex
optimization problem, both in terms of the closest point between shapes [11],
and in terms of the minimum scale factor [21]. For the former, we introduce
two points in a world frame $p_{i}\in\mathbf{R}^{3}$, and two polytopes
described with $A_{i}p_{i}\leq b_{i}$. By constraining each point to be within
a polytope, a QP is used to solve for the closest point between these two
shapes, p_1, p_2 ∥p_1 - p_2∥ A_1 p_1≤b_1A_2 p_2≤b_2. Using a differentiable QP
solver, the gradient of the closest distance between shapes with respect to
the positions of the polytopes can be computed resulting in the contact normal
vectors.

In this example, we examine collision detection between two squares and the
behavior of these contact normals in the presence of sharp corners. As shown
in Fig. 4, the contact normals are evaluated at a strict $\kappa=0$ and a
relaxed $\kappa=0.01$. In the case of $\kappa=0$, the contact normals are
(correctly) exactly normal to the surface, and as soon as the closest point
shifts from one face to the other, the contact normals immediately rotate
$90^{\circ}$. While this is expected behavior, this discontinuity in the
gradient can prove troublesome for simulation and control algorithms that rely
on these contact normals not changing too quickly. Alternatively, with a
relaxed $\kappa=0.01$, the contact normal smoothly rotates the 90 degrees as
the face of the closest point changes. This is a result of the logarithmic
barrier smoothing out the sharp corner, and allows for continuous and smooth
gradients even in the presence of the discontinuity.

##  6 Conclusions

In this paper, we outline shortcomings with existing differentiable
optimization tools, namely the nonsmoothness of the gradients near inequality
constraints and the inability to handle infeasible problems, and propose
solutions to both of these problems. By relaxing the solution to an
optimization problem from tight tolerances to an intentionally relaxed
logarithmic barrier, unique and smooth gradients can be computed even from
sharp edges in the feasible set. This relaxation is straightforward and
leverages existing routines within existing primal-dual interior-point
solvers. We also introduce an always-feasible quadratic program where hard
constraints are converted into $\ell_{1}$-penalties, and devise a customized
algorithm for solving problems of this form with limited added computational
overhead. With both of these innovations, consistent and reliable smooth
gradients are demonstrated in common robotic tasks where smoothness is a
priority. Our fully differentiable and parallelizable solver written in JAX is
available at www.github.com/kevin-tracy/qpax.

####  6.0.1

The authors have no competing interests to declare that are relevant to the
content of this article..

## References

  * [1] Agrawal, A., Amos, B., Barratt, S., Boyd, S., Diamond, S., Kolter, Z.: Differentiable Convex Optimization Layers. Advances in Neural Information Processing Systems pp. 9558–9570 (2019)
  * [2] Agrawal, A., Barratt, S., Boyd, S., Moursi, W.M., Busseti, E.: Differentiating Through a Conic Program. Journal of Applied and Numerical Optimization 1(2), 107–115 (2019)
  * [3] Amos, B., Kolter, J.Z.: Optnet: Differentiable optimization as a layer in neural networks. In: International Conference on Machine Learning. pp. 136–145. PMLR (2017)
  * [4] Amos, B., Kolter, J.Z.: OptNet: Differentiable Optimization as a Layer in Neural Networks. arXiv:1703.00443 [cs, math, stat] (Oct 2019)
  * [5] Anitescu, M.: Optimization-based simulation of nonsmooth rigid multibody dynamics. Mathematical Programming 105(1), 113–143 (Jan 2006). https://doi.org/10.1007/s10107-005-0590-7
  * [6] Blackmore, L.: Autonomous Precision Landing of Space Rockets. The Bridge 4(46), 15–20 (2016)
  * [7] Boyd, S., Vandenberghe, L.: Convex Optimization. Cambridge University Press (2004)
  * [8] Diamond, S., Boyd, S.: CVXPY: A Python-Embedded Modeling Language for Convex Optimization
  * [9] Domahidi, A., Morari, M., Morari, M., Boyd, S.P., Boyd, S.P.: Methods and Tools for Embedded Optimization and Control. Ph.D. thesis, ETH, Zürich (2013)
  * [10] Gilbert, E., Chong Jin Ong: New distances for the separation and penetration of objects. In: Proceedings of the 1994 IEEE International Conference on Robotics and Automation. pp. 579–586. IEEE Comput. Soc. Press, San Diego, CA, USA (1994). https://doi.org/10.1109/ROBOT.1994.351237
  * [11] Gilbert, E., Johnson, D., Keerthi, S.: A fast procedure for computing the distance between complex objects in three-dimensional space. IEEE Journal on Robotics and Automation 4(2), 193–203 (Apr 1988). https://doi.org/10.1109/56.2083
  * [12] Gill, P.E., Murray, W., Saunders, M.A.: SNOPT: An SQP Algorithm for Large-Scale Constrained Optimization. SIAM Review 47(1), 99–131 (Jan 2005). https://doi.org/10.1137/S0036144504446096
  * [13] Howell, T.A., Cleac’h, S.L., Kolter, J.Z., Schwager, M., Manchester, Z.: Dojo: A Differentiable Physics Engine for Robotics. http://arxiv.org/abs/2203.00806 (Jun 2022). https://doi.org/10.48550/arXiv.2203.00806
  * [14] Kirchengast, M., Steinberger, M., Horn, M.: A Fast Quadratic Program Solver for Constrained Control Allocation. In: 2018 IEEE Conference on Decision and Control (CDC). pp. 4065–4071. IEEE, Miami Beach, FL (Dec 2018). https://doi.org/10.1109/CDC.2018.8619828
  * [15] Kuindersma, S., Permenter, F., Tedrake, R.: An Efficiently Solvable Quadratic Program for Stabilizing Dynamic Locomotion. In: Proceedings of the International Conference on Robotics and Automation (ICRA). pp. 2589–2594. IEEE, Hong Kong, China (May 2014). https://doi.org/10.1109/ICRA.2014.6907230
  * [16] Mattingley, J., Boyd, S.: CVXGEN: A code generator for embedded convex optimization. In: Optimization Engineering. pp. 1–27 (2012)
  * [17] Mehrotra, S.: On the Implementation of a Primal-Dual Interior Point Method. SIAM Journal on Optimization 2(4), 575–601 (Nov 1992). https://doi.org/10.1137/0802028
  * [18] Nocedal, J., Wright, S.J.: Numerical Optimization. Springer, 2nd edn. (2006)
  * [19] Pang, T., Tedrake, R.: A Convex Quasistatic Time-stepping Scheme for Rigid Multibody Systems with Contact and Friction. In: 2021 IEEE International Conference on Robotics and Automation (ICRA). pp. 6614–6620. IEEE, Xi’an, China (May 2021). https://doi.org/10.1109/ICRA48506.2021.9560941
  * [20] Stellato, B., Banjac, G., Goulart, P., Bemporad, A., Boyd, S.: OSQP: An Operator Splitting Solver for Quadratic Programs p. 40
  * [21] Tracy, K., Howell, T.A., Manchester, Z.: Differentiable Collision Detection for a Set of Convex Primitives. In: 2023 IEEE International Conference on Robotics and Automation (ICRA). pp. 3663–3670. IEEE, London, United Kingdom (May 2023). https://doi.org/10.1109/ICRA48891.2023.10160716
  * [22] Tracy, K., Manchester, Z., Jain, A., Go, K., Schaal, S., Erez, T., Tassa, Y.: Efficient Online Learning of Contact Force Models for Connector Insertion (Dec 2023)
  * [23] Vandenberghe, L.: The CVXOPT linear and quadratic cone program solvers p. 30
  * [24] Varin, P., Kuindersma, S.: A Constrained Kalman Filter for Rigid Body Systems with Frictional Contact. In: Morales, M., Tapia, L., Sánchez-Ante, G., Hutchinson, S. (eds.) Algorithmic Foundations of Robotics XIII, vol. 14, pp. 474–490. Springer International Publishing, Cham (2020). https://doi.org/10.1007/978-3-030-44051-0_28
  * [25] Xinjilefu, X., Feng, S., Atkeson, C.G.: Dynamic state estimation using Quadratic Programming. In: 2014 IEEE/RSJ International Conference on Intelligent Robots and Systems. pp. 989–994. IEEE, Chicago, IL, USA (Sep 2014). https://doi.org/10.1109/IROS.2014.6942679