# Quadratic Advantage with Quantum Randomized Smoothing Applied to Time-Series
Analysis

Nicola Franco2, Marie Kempkes3, Jakob Spiegelberg3, Jeanette Miriam Lorenz2
2Fraunhofer Institute for Cognitive Systems IKS, Munich, Germany
{nicola.franco<EMAIL_ADDRESS>3Volkswagen Group
Innovation, Volkswagen AG, Wolfsburg, Germany
{marie.kempkes<EMAIL_ADDRESS>
###### Abstract

As quantum machine learning continues to develop at a rapid pace, the
importance of ensuring the robustness and efficiency of quantum algorithms
cannot be overstated. Our research presents an analysis of quantum randomized
smoothing, how data encoding and perturbation modeling approaches can be
matched to achieve meaningful robustness certificates. By utilizing an
innovative approach integrating Grover’s algorithm, a quadratic sampling
advantage over classical randomized smoothing is achieved. This strategy
necessitates a basis state encoding, thus restricting the space of meaningful
perturbations. We show how constrained $k$-distant Hamming weight
perturbations are a suitable noise distribution here, and elucidate how they
can be constructed on a quantum computer. The efficacy of the proposed
framework is demonstrated on a time series classification task employing a
Bag-of-Words pre-processing solution. The advantage of quadratic sample
reduction is recovered especially in the regime with large number of samples.
This may allow quantum computers to efficiently scale randomized smoothing to
more complex tasks beyond the reach of classical methods.

###### Index Terms:

Quantum Machine Learning, Adversarial Robustness, Randomized Smoothing.

##  I Introduction

The integration of quantum computing (QC) principles into machine learning,
forming the field of quantum machine learning (QML), represents a
groundbreaking shift in the landscape of computational intelligence [1, 2, 3].
Theoretically, QML offers remarkable potential benefits, including the ability
to handle exponentially larger datasets and to perform certain computations
much faster than classical algorithms [4, 5]. However, QML is not without its
vulnerabilities, particularly in the realm of adversarial attacks [6, 7, 8].
These are instances where input data is deliberately manipulated to deceive
the model into making incorrect predictions or classifications, a challenge
well-documented in the classical machine learning domain [9, 10]. In response
to these adversarial threats, both quantum and classical machine learning
communities have developed probabilistic [11, 12, 13, 14, 15] and
deterministic [16, 17, 18, 19] defense strategies.

Defenses are typically developed as either formal verification of outcomes via
input space propagation or by utilizing sampling distributions linked to
differential privacy. While the latter method has shown promise for enhancing
QML applications, it is not without its limitations. For instance, the
assumptions made by Weber et al. [13] about class probability estimates
overlook the complexities of machine learning uncertainties, rendering their
robustness certificates potentially overly optimistic. Moreover, their claims
of robustness face challenges due to issues with data encoding and
scalability, raising doubts about their usefulness in practical scenarios.
Similarly, Du et al.’s [20] work is limited by its focus on depolarization
noise, thereby restricting the range of applicable data encodings. Their
approach does not accurately represent realistic perturbation distributions,
particularly for basis state encodings, leading to a mismatch between
theoretical robustness and practical applicability. Huang et al. [14] while
considering rotation noise, also suffers from similar shortcomings, with
robustness certifications heavily dependent on data encoding types and
potentially inadequate for large-scale quantum systems. A more detailed
discussion of these works can be found below in section III.

These critiques highlight a gap between the theoretical robustness measures
proposed and their practical utility across various QML applications. The
failure to integrate perturbation types with data encodings may hinder the
achievement of realistic robustness certifications. Our method aligns with
QuAdRo [15] due to its flexibility regarding the algorithm and noise type, and
its clear demonstration of how QC reduces the sample size needed for
robustness certification. However, it also has limitations: it requires data
to be encoded in basis states due to the use of the Grover oracle, limiting
data encoding options. On this background, we investigate how data encoding
and perturbation types have to be matched to get meaningful, and with QuAdRo
scalable, robustness certificates. Our contributions can be summarized as
follows:

  * •

The connection between information encoding and perturbation type is discussed
at length, we show that the interpretation of certificates obtained with
randomized smoothing changes significantly and might even be meaningless for
some combinations.

  * •

Constrained $k$-distant Hamming weight perturbations are identified as
meaningful perturbation type compatible with a basis state encoding as
necessary for QuAdRo. Additionally, an algorithm for implementing them on a
quantum computer is introduced.

  * •

We demonstrate quadratic sampling advantage using QuAdRo on a practical time
series classification task by utilizing an established Bag-of-Words pre-
processing technique.

##  II Preliminaries

###  II-A QC background

QC operates quite differently from traditional computing. At its core, QC uses
the quantum bit, or qubit, symbolized as $\ket{\psi}$. What makes a qubit
unique is its capacity to simultaneously represent two states, $\ket{0}$ and
$\ket{1}$, thanks to a property called superposition, which is described by
the equation $\ket{\psi}=\psi_{1}\ket{0}+\psi_{2}\ket{1}$. Superposition
allows a qubit to exist in a complex space, the Hilbert space, with dimensions
that grow exponentially as $2^{n}$, where $n$ is the total number of qubits.

In contrast to classical bits that represent a single state at a time, qubits’
ability to embody dual states simultaneously allows quantum computers to
process multiple computations in parallel. This processing is executed through
quantum gates, depicted by unitary matrices $U$. These matrices are crucially
unitary, meaning the complex conjugate transpose of $U$, denoted as
$U^{\dagger}$, is also its inverse, fulfilling the condition $U^{\dagger}U=I$.
This ensures the preservation of the inner product between vectors. However,
extracting a conventional outcome from quantum processes necessitates
measuring the qubits, which results in the collapse of their superposed state.
In essence, practically speaking, QC involves manipulating the quantum state
of qubits to, for example, increase the likelihood of obtaining a specific
solution.

###  II-B Embedding classical data into quantum states

In QC, a prevalent method for integrating classical data into quantum states
is called basis embedding. This technique converts classical data, formatted
as bit strings, into representations within quantum states.

######  Definition II.1 (Basis Embedding).

Given a sequence of bits $\mathbf{a}=a_{1}a_{2}...a_{n}$, basis embedding is
the process of translating this sequence into the form
$\ket{\mathbf{a}}=\ket{a_{1}}\otimes\ket{a_{2}}\otimes...\otimes\ket{a_{n}}$,
with $\ket{a_{1}a_{2}...a_{n}}$ serving as a more compact notation.

Additionally, QC utilizes amplitude embedding, which encodes data into the
amplitudes of quantum states, further expanding the range of possible
computations.

######  Definition II.2 (Amplitude Embedding).

Given a normalized vector $\mathbf{a}=(a_{1},a_{2},...,a_{N})$ where
$a_{i}\in\mathbb{R}$, $N=2^{n}$ and $\sum_{i=1}^{N}|a_{i}|^{2}=1$, amplitude
embedding is the process of mapping this vector into a quantum state
$\ket{\mathbf{a}}=a_{1}\ket{0}+a_{2}\ket{1}+...+a_{N}\ket{N-1}$, where each
amplitude $a_{i}$ corresponds to the probability amplitude of the quantum
state in the computational basis states $\ket{0},\ket{1},...,\ket{N-1}$.

Other embedding methods, such as angle and Hamiltonian embedding, are also
employed in QC to represent data in various quantum-friendly formats.

###  II-C Certifiable robustness in machine learning

Consider a binary classifier $f:\mathcal{X}\rightarrow\\{0,1\\}$ that assigns
an input vector $\mathbf{x}\in\mathcal{X}$ (Euclidean space) to a binary class
$f(\mathbf{x})\in\\{0,1\\}$. Additionally, we define as
$y:\mathcal{X}\mapsto[0,1]$ its soft version which computes the model’s logit
value and the class prediction is made through $\mathbf{1}[y>0.5]$. Further,
let $o:\mathcal{X}\times\mathcal{X}\mapsto\\{0,1\\}$ denote the oracle that
compares the semantic content of two input data points and returns $1$ if they
are semantically the same. The research area of certifiable (provable)
robustness aims to establish formal guarantees that under all admissible
attacks within the budget, none will alter the model’s prediction. Formally,
certification at a point $\mathbf{x}$ involves evaluating if,

$f(\mathbf{\tilde{x}})=f(\mathbf{x})\text{ for all
}\mathbf{\tilde{x}}\in\mathcal{X}\ \text{s.t.}\
o(\mathbf{\tilde{x}},\mathbf{x})=1.$ (1)

The assessment of semantic similarity between data points is complex and
varies by task. This complexity is mitigated by considering two inputs
semantically similar if their $\ell_{p}$-norm difference is less than a small
threshold $\epsilon$, translating the robustness certification problem at
$\mathbf{x}$ into:

$f(\mathbf{\tilde{x}})=f(\mathbf{x})\text{ for all
}\mathbf{\tilde{x}}\in\mathcal{X}\ \text{s.t.}\
\lvert\lvert\mathbf{\tilde{x}}-\mathbf{x}\rvert\rvert_{p}\leq\epsilon.$ (2)

For neural networks, robustness verification can be structured as a mixed-
integer linear program (MILP) based on the network’s weights [21]. However,
verifying neural network properties, especially with ReLU activations, is
proven to be NP-complete [22]. As a workaround, approximation methods are
used, which might not certify all non-altering perturbations [23]. Randomized
Smoothing (RS) [11, 12] has emerged as a notable method in the relaxed
certification category, enabling robustness assessment independently of the
classifier’s details.

###  II-D Randomized smoothing

Randomized Smoothing as proposed by [11, 12], constructs a smoothed classifier
$g$ by evaluating the original classifier $f$’s logits $y$ around the input,
defined as:

$\displaystyle
g(\mathbf{x})=\mathbb{E}_{\mathbf{z}\sim\mathcal{N}(0,\sigma^{2}I)}[f(\mathbf{x}+\mathbf{z})].$
(3)

Here, the focus shifts from a direct prediction at $\mathbf{x}$ to the
expected prediction over noisy samples distributed normally around
$\mathbf{x}$. In practice, the goal is to determine the minimum value of $g$
over a ball of radius $\epsilon$ centered at zero in the $\ell_{p}$ space, and
check if this minimum value exceeds $\frac{1}{2}$, formally defined as:

$\min_{\bm{\delta}\in\mathcal{B}_{p}(\epsilon)}g(\mathbf{x}+\bm{\delta})\
\text{s.t.}\
\mathcal{B}_{p}(\epsilon):=\\{\bm{\delta}\in\mathcal{X}:\lvert\lvert\bm{\delta}\rvert\rvert_{p}\leq\epsilon\\}.$
(4)

Thus, let $\Phi$ denote the CDF of the standard normal distribution
$\mathcal{N}(0,1)$, and $\underline{p_{1}}$ be a lower bound of Eq. 4, then
$g(\mathbf{x})$ as constructed in Eq. 3 is certifiably robust if
$\epsilon<\sigma\Phi^{-1}(\underline{p_{1}})$ for $1D$ classifier, i.e.,
$\mathcal{X}=\mathbb{R}$ [11].

##  III Robustness in QML

TABLE I: List of encoding types, noise types and their respective input space interpretation. An indication of suggested pairings is given. Note that we do not further specify constraints (e.g., fixed Hamming weight distance for bitflip noise). We assume initialisation in the $\ket{000\dots}$-state.  Encoding | Perturbation | Input Space Interpretation | Suggested
---|---|---|---
amplitude | depolarization noise | average feature contrast diminished | ✓
| amplitude damping | not intuitive |
| gaussian noise | noisy features, $l_{2}$ adversarial defense | ✓
| uniform noise | noisy features, $l_{1}$ adversarial defense | ✓
| rotation noise | not intuitive |
| bitflip noise | permutation of features |
rotation | depolarization noise | feature contrast diminished and occurrence of previously not accessed states |
| amplitude damping | pushes all features to $0$ |
| gaussian noise | noisy features and occurrence of previously not accessed states |
| uniform noise | noisy features and occurrence of previously not accessed states |
| rotation noise | noisy features | ✓
| bitflip noise | not intuitive |
basis state | depolarization noise | all other data values gain small likelihood |
| amplitude damping | not intuitive |
| gaussian noise | all other data values gain small likelihood stochastically |
| uniform noise | all other data values gain small likelihood stochastically |
| rotation noise | not intuitive |
| bitflip noise | randomly shift to other value, can be constrained to nearby values | ✓

Various approaches for certifiable robustness have been suggested for QML. In
this extended section, we outline known methods for issuing robustness
certificates and highlight their usage of perturbations. This leads to a
conversation about the importance of jointly designing data encoding and
perturbations.

###  III-A Existing work on robustness certificates in QML

Weber et al. [13] draw a theoretical correlation between binary quantum
hypothesis testing (QHT) and certifiable robustness: while QHT aims at
optimally distinguishing two quantum states, robust classification seeks to
provide provable guarantees that two states are not distinguishable. Based on
this connection, they derive a robustness bound for quantum classifiers
subject to worst-case input noise, which they prove to be tight for binary
classification. This bound, however, only holds in the regime where the
probability of measuring the true class probability $p_{A}$ is always larger
than the runner-up probability estimate $p_{B}$, which boils down to assuming
perfect classification accuracy in an unperturbed setting. Such an assumption
is common and necessary within the field of randomized smoothing to obtain
robustness certifications, as seen in prior works [11, 12]. However, this
assumption may not be reasonable for machine learning models that are often
overly confident and vulnerable to various uncertainties in practical
applications. To address the limitations of assuming an ideal classifier,
alternative approaches, e.g. conformal prediction [24, 25], might be required.

Another critical aspect we want to emphasize is the need of carefully
considering the type of data encoding for the robustness bounds in [13] to
provide meaningful insights on the trust that can be given to the prediction
of a perturbed data point. More precisely, as the bounds are based on the
fidelity between input states, only non-zero fidelities allow useful
conclusions on robustness. For example, when using basis state encoding, even
the smallest changes in the input data result in zero fidelity, thus
prohibiting any statement about predictions based on altered data. While this
aspect can be easily resolved by employing, e.g., the widely used rotation
encoding, there is one more issue we want to highlight, which is related to
the scalability of fidelity. As demonstrated by [26], fidelities between
random states vanish exponentially in the number of qubits and hence even
perturbations of small amplitude can cause a severe drop of certified
robustness. This includes disturbances that lead to a low level of fidelity,
yet may not directly impact the outcome of the algorithm we are analyzing. The
last point of criticism is best explained with a pathological example:
Consider a binary classification algorithm $\mathcal{A}$ defined on a
parameterized quantum circuit with two distinct subsystems, denoted $A$ and
$B$. The algorithm maps inputs $x$ to label 1 if the probability of one of the
qubits in subsystem $A$ being in the $\ket{1}$ state is larger than the
corresponding probability for subsystem $B$:

$\mathcal{A}(x)=\begin{cases}1,&\text{if
$\bra{\Psi}M_{A}^{\dagger}M_{A}\ket{\Psi}>\bra{\Psi}M_{B}^{\dagger}M_{B}\ket{\Psi}$}.\\\
0,&\text{otherwise}.\end{cases}$ (5)

For simplicity, consider three qubits with the first two constituting
subsystem $A$ and the last one subsystem $B$. The measurement operators then
take the form:

$\displaystyle M_{A}$ $\displaystyle=\bra{100}+\bra{110}+\bra{010}+\bra{000}$
$\displaystyle M_{B}$ $\displaystyle=\bra{101}+\bra{111}+\bra{011}+\bra{001}$

Assuming the state $\ket{100}$ is to be certified, a perturbation could
manipulate the state as far as $\ket{010}$, resulting in a fidelity of zero
and thus a trivial certification radius of zero according to Weber et al.
[13]. At the same time, the algorithm would still be no less certain to return
label 1 (the perturbed state is no closer to the decision boundary
$\bra{\Psi}M_{a}^{\dagger}M_{a}\ket{\Psi}>\bra{\Psi}M_{b}^{\dagger}M_{b}\ket{\Psi}$).

This example is admittedly constructed, but consider generalized versions of
it on larger systems. For instance, this might include a classification task
over many input features, where only a few are indicative of the label to be
returned, e.g., because only part of an image contains the object which we try
to classify. Then there is lots of room for perturbations to push the state
somewhere the classification outcome is unaffected, although fidelities might
drop even exactly to zero. This impedes the utility of fidelity based
robustness bounds.

In summary, the bounds derived in [13], while valid, can be expected to be
very conservative in practice, abstaining from certification in many cases
where fidelities drop, simply because of the encoding used or the size of the
system. This calls out for extensions which give tighter, more useful
robustness statements in practice.

A different line of research investigates the impact of various types of
hardware noise on robustness bounds. Huang et al. [14] show that random
rotation noise can improve robustness in quantum classifiers against
adversarial attacks by establishing a connection to differential privacy. A
certified robustness radius is derived in terms of the output probability
ratio for binary classification given that the trace distance between
perturbed and unperturbed state is smaller than some threshold. Using similar
scaling arguments as before, we remark that the trace distance between two
random quantum states vanishes exponentially in the system size and thus the
condition for the robustness bound is violated easily for a large number of
qubits.

In a similar work, it was shown for depolarization noise that robustness of
quantum classifiers can be enhanced by increasing levels of noise [20]. Here
as well we argue that the type of noise limits the data encodings which can be
used meaningfully. While it has a reasonable interpretation for amplitude
encoding, this type of noise has a questionable, highly non-trivial effect on
basis state encoded data. We further discuss the co-design of perturbation and
encoding in Section III-B.

[15] approach quantum randomized smoothing from a slightly different angle.
They introduce an algorithm based on Grover’s search providing robustness
bounds with provable quadratic sample advantage over classical methods. Our
approach closely follows [15] for several reasons. First, it is flexible in
terms of the algorithm used and the type of noise. Second, it clearly
motivates the use of quantum computers due to the reduced number of samples
necessary for a given robustness certificate. However, we do note some
shortcomings as well. First, the Grover oracle necessitates a basis state
encoding of the data, which limits the range of data encoding strategies
available. Second, not all relevant perturbation types can be
straightforwardly implemented. At last, the interplay of encoding and
perturbation type is not considered.

###  III-B Co-Design of data encoding and perturbation type

In this article, we investigate QML models applied to classical data. As such,
attacks and perturbations are applied in input space, and the perturbed $x$ is
encoded into the quantum circuit. Quantum randomized smoothing introduces
disturbances to the quantum circuit to achieve a more reliable prediction by
analyzing the prediction statistics from multiple altered versions of the same
data point.

If one is interested in a worst case robustness, as in adversarial defense,
bounded perturbations are assumed and an unstructured perturbation is applied,
e.g., Gaussian or uniform noise. If randomized smoothing is performed with
other, potentially more structured perturbations, valid certificates can still
be obtained, although certification radii can become overly conservative, as
worst case robustness often already is. To avoid overly cautious certificates,
a perturbation applied on the quantum circuit should be designed such that the
perturbed state could be obtained by directly encoding a perturbed input data
point, i.e., there should be an equivalent perturbation applied directly to
the input data domain that produces the same state as when applying the
perturbation on the quantum computer to the encoded state of the unperturbed
data.

A second case is the certification of robustness to certain types of
perturbations or noises. If these are resulting from the quantum hardware,
perturbations which do not have an equivalent input domain perturbation may be
useful, e.g., as they directly reflect the hardware noise. We put our focus on
the more common case of certifying input space robustness. In this case,
again, perturbations on the quantum circuit do need to have an input space
interpretation.

With this necessity in mind, consider Tab. I where a selection of encoding and
perturbation types is listed. For each combination, an interpretation of the
induced changes on the state vector in input space is indicated. We leave it
to the reader to work through all combinations and only give one example here:
Consider amplitude encoded data acted upon by a depolarization channel. We can
write the perturbed state as

$\displaystyle\ket{\tilde{\Psi}}$
$\displaystyle=(1-p)\sum_{i}^{N}a_{i}\ket{i}+\frac{p}{N}\mathbf{I}$
$\displaystyle=\sum_{i}^{N}\left((1-p)a_{i}+\frac{p}{N}\right)\ket{i}$
$\displaystyle=\sum_{i}^{N}\tilde{a}_{i}\ket{i}$

from which we can see that the feature contrast, the difference between any
two features $\tilde{a}_{i}$ and $\tilde{a}_{j}$ is diminished by a factor of
$1-p$ compared to the contrast of unperturbed $a_{i}$ and $a_{j}$ and
disappears completely as $p$ approaches $1$. Additionally, we see that
$\ket{\tilde{\Psi}}$ can again be written as an amplitude encoded state.
Hence, there exists an equivalent classical perturbation which could be used
to obtain the $\tilde{a}_{i}$ coefficients.

As highlighted in the last column of Tab. I, only few of the combinations of
data encoding and perturbation type suffice the requirement of
interpretability outlined above. Many others either access states which could
not be reached by an equivalent classical perturbation with subsequent
encoding, or lack a clear intuition of the induced changes. In quantum
randomized smoothing, we suggest to be aware of the interplay of encoding and
perturbation and design them simultaneously so that meaningful and not
unnecessarily loose certificates are obtained. Furthermore, note that for a
basis state encoding only bitflip noise allows a clear interpretation. We thus
dedicate the next section to its implementation with Hamming weight constrains
for later usage in QuAdRo.

##  IV $k$-Hamming Distant States

In this section, we describe how to construct a distribution around a binary
input $\mathbf{x}\in\\{0,1\\}^{n}$ of $k$-Hamming distant states. This differs
from the method of preparing Dicke states, where a Dicke state is an equal-
weight superposition of all $n$-qubit states with Hamming weight $k$ [27]. In
our approach, we prepare a state
$\psi:\\{0,1\\}^{n}\to\mathcal{D}(\mathbb{H})$ that depends on $\mathbf{x}$.

######  Definition IV.1 ($k$-Hamming distant states).

Given a binary input $\mathbf{x}\in\\{0,1\\}^{n}$, we define a constrained
$k$-Hamming distance set as:

$\mathcal{D}_{n,k}(\mathbf{x})=\\{\tilde{\mathbf{x}}\in\\{0,1\\}^{n}\,:\,||{\mathbf{x}-\tilde{\mathbf{x}}}||_{0}=k\\},$
(6)

where the zero norm111the zero ”norm” is not truly a norm. corresponds to the
non-zero elements differing from $\mathbf{x}$.

For any given positive integer Hamming weight $k$, with $k\leq n$, the number
of states that are exactly $k$-Hamming distances away from a given state is
given by the binomial coefficient $\binom{n}{k}$. We can formalize the
approach for allocating probabilities to states based on their Hamming
distances from the input state, including a maximum probability for the input
state itself ($i=0$) and decreasing probabilities for states with increasing
Hamming distances ($i=1$ to $k$). To do this, we need to ensure that:  (i) the
input state ($i=0$) has the highest probability,  (ii) probabilities for
states decrease as their Hamming distance from the input state increases and
(iii) probabilities are equally distributed among states at the same Hamming
distance.  . First, we define a base weight of $e^{-\frac{i}{\sigma}}$ for
each Hamming distance $0\leq i\leq k$, including a special consideration for
$i=0$ to ensure it receives the maximum weight of 1. Here, $\sigma>0$ works as
hyper-parameter that adjusts the distribution of weights. A larger $\sigma$
would result in a slower decrease in weights with increasing $i$, allowing for
more uniform distribution of probabilities across different Hamming distances.
Thus, the normalized weight for each Hamming distance $i$ is given by:

$w(i)=\frac{e^{-\frac{i}{\sigma}}}{\sum_{i=0}^{k}e^{-\frac{i}{\sigma}}}\
\text{for all}\ 0\leq i\leq k.$ (7)

Finally, to assign probabilities to individual states at a given Hamming
distance $i$, each state receives a fraction of the total probability
allocated to that distance:

$p(i)=\frac{w(i)}{\binom{n}{i}}\quad\text{for each state at distance }i.$ (8)

As an example, let us consider an input vector $\mathbf{x}=(0,1,1)$ which has
a length of $n=3$; there are exactly three states that are Hamming distance of
1 away from it, three of distance 2 and one state of distance 3. For
$\sigma=1$, the distribution of probabilities for states based on their
Hamming distance is: 0) the initial state with probability of $0.644$,  1. 3
states, each with a probability of $0.079$,  2. 3 states, each with a
probability of $0.029$ and  3. 1 state with a probability of $0.032$.

###  IV-A Construction of $k$-Hamming distant states

Given an input $\mathbf{x}\in\\{0,1\\}^{n}$, we construct a superposition of
$k$-Hamming distant states from $\mathbf{x}$ in the Hilbert space as:

$\ket{\psi(\mathbf{x})}=\sum_{i=0}^{k}\sqrt{p(i)}\sum_{\tilde{\mathbf{x}}\in\mathcal{D}_{n,i}(\mathbf{x})}\ket{\tilde{\mathbf{x}}}.$
(9)

These can be understood as the probabilities of staying in the original state
$\mathbf{x}$ or transitioning to the nearest state that is $k$-Hamming distant
away.

A direct approach to construct the quantum state $\ket{\psi(\mathbf{x})}$
involves first calculating the probabilities associated with the $k$-Hamming
states. Following this step, the quantum circuit is prepared using the
Möttönen state preparation method [28], a widely recognized technique for
incorporating data into the amplitudes of quantum states. Despite its
effectiveness, this approach encounters a significant limitation due to the
exponential increase in the depth of the circuit, symbolized by
$\mathcal{O}(2^{n})$, which pertains to the exponential growth in the
requirement for both CNOT and single-qubit rotational gates as the number of
qubits increases. This exponential increase presents an additional challenge,
specifically when attempting to enumerate classically all the probabilities
corresponding to states that are $k$-Hamming distances apart. As the
dimensionality of the input grows, this task becomes increasingly
unmanageable, highlighting a critical bottleneck in the scalability of this
method for larger quantum systems.

(a) Circuit preparation with controlled RY($\sigma\pi$). (b) Quasi-probability
distributions for varying $\sigma$.

Figure 1: 1-Hamming distant state preparation for $\mathbf{x}=(0,1,1)$.

#### Our state preparation

To construct the quantum state $\ket{\psi(\mathbf{x})}$ effectively, we
propose a method that includes the use of extra qubits, termed ancillary
qubits, along with controlled rotational gates. This method is enhanced by
introducing a technique that involves the inverse application of multi-
controlled basis embedding for the input. This enables the design of a more
efficient circuit that requires the addition of only one qubit for every input
dimension. In Fig. 1, we showcase the procedure for setting up a circuit for a
three-dimensional input. The process begins with the initial preparation of
the $x$ state using a basis embedding technique. Following this, $n$ ancillary
qubits are introduced, and a combination of Hadamard gates and controlled RY
gates, adjusted by $\sigma\pi$, are applied. In this scenario,
$\sigma\in[0,1]$ and an increase of its value denotes an overall increase on
the likelihood of neighbouring states. We observe that for $\sigma=0.5$ the
controlled RY rotation corresponds to a controlled Hadamard, hinting the
quantum Fourier transform state preparation. Finally, the ancillary qubits are
manipulated with a series of PauliX gates before the input undergoes a dual
embedding in the basis via a multi-controlled gate. Inverting the ancillary
qubits prior to introducing a controlled input decreases the likelihood of
states that are distant from the input.

Unlike the Möttönen state preparation method, our technique does not replicate
the distribution with exact precision. In practical terms, this situation
means that while our method can effectively adjust the likelihood of certain
quantum states being realized - specifically, those within a $k$-Hamming
distance from a given state - it also unintentionally boosts the chances of
encountering states beyond this specified Hamming distance. Essentially, when
we attempt to increase the probability of closely related states (those a few
quantum flips away), we inadvertently also make it more likely to observe
states that are further away than intended. This can introduce inaccuracies in
quantum simulations or computations that rely on precise probability
distributions.

(a) Circuit preparation for a uniform distribution. (b) Quasi-probability
distributions for varying $\sigma$.

Figure 2: Uniform distant state preparation for $\mathbf{x}=(0,1,1)$.

#### Uniform state preparation

In Fig. 2, we detail the process for setting up a uniform distribution
centered on a three-dimensional input along with its corresponding
probabilities. Continuing from the earlier discussion, the probabilities for
adjacent states are uniformly likely across all states. This method of state
preparation is more efficient since it only necessitates a single ancillary
qubit, and we can also adjust $\sigma$ to boost the probability of nearby
qubits through a RY gate on the ancillary. Nonetheless, treating the
probability of neighboring states as equally likely for all states presents a
drawback when considering semantically meaningful perturbations for the input
state.

###  IV-B Certified robustness for $k$-Hamming distant states

Here, we formally describe the computation of the certified robustness
distance in relation to a distribution that is $k$-Hamming distant from the
given input. We focus on the robustness certificate based on the
$\ell_{0}$-norm within a binary input domain $\mathcal{X}=\\{0,1\\}^{n}$.
Specifically, let $p_{i}:\mathcal{X}\to[0,1]$ denote the likelihood that
smoothing the base classifier $f$ with a $k$-Hamming distant distribution
around $\mathbf{x}$ yields class $i$, formally:

$p_{i}(\mathbf{x})=\mathbb{E}_{\tilde{\mathbf{x}}\in\mathcal{D}_{n,k}(\mathbf{x})}[f(\tilde{\mathbf{x}})=i],$
(10)

where $\mathcal{D}_{n,k}(\mathbf{x})$ represents the set of $k$-Hamming
distant vectors for $\mathbf{x}$. As established by [29], for any input
$\tilde{\mathbf{x}}$ with $||{\mathbf{x}-\tilde{\mathbf{x}}}||_{0}\leq k$
across all classes $i\in\mathcal{K}$, the following holds:

$|p_{i}(\tilde{\mathbf{x}})-p_{i}(\mathbf{x})|\leq\Delta,\;\text{where}\;\Delta=\operatorname{\rm{argmax}}_{r}p_{r}\geq
0.5,$ (11)

defines the certified robustness distance. Based on the architecture of the
base classifier, directly computing $p_{i}(\mathbf{x})$ might not be feasible.
However, we can instead generate a finite set of samples from
$\mathcal{D}_{n,k}$, in order to estimate a lower bound
$\underline{p_{i}(\mathbf{x})}$ with high confidence ($1-\alpha$). In
practice, for an input $\tilde{\mathbf{x}}$ with
$||{\mathbf{x}-\tilde{\mathbf{x}}}||_{0}\leq k$,

$\text{if}\quad\underline{p_{i}(x)}-\Delta>\frac{1}{2},\quad\text{then}\quad\operatorname{\rm{argmax}}_{j}p(x)_{j}=i,$

with probability at least $1-\alpha$ [11].

In practice, to establish a lower bound for $\underline{p_{i}(\mathbf{x})}$,
traditional methods employ the Clopper-Pearson interval [11]. Nonetheless, the
quantity of samples, which corresponds to the number of requests made to the
base classifier, increases following an O(1/$\alpha^{2}$) complexity, where
$\alpha$ represents the target error margin within the expected probability
value, delineated over the randomized smoothing area surrounding the input. In
the context of QML, [15] demonstrate that for achieving the same level of
confidence and error margin, the requisite number of inquiries to the base
classifier is reduced to $\mathcal{O}(1/\alpha)$. QuAdRo [15] employing $m$
ancilla qubits for the quantum phase estimation allows for a maximum of
$M=2^{m+1}-1$ oracle calls. In addition, the likelihood of success approaches
nearly 100% by conducting multiple trials and selecting the median outcome.

##  V Experimental Results

###  V-A Setup

Figure 3: Binary representation of discretized patterns in time series,
visualized through the Bag-of-Words representation. The method operates by
splitting in half a sample from the gun-point [30] time series dataset, then
dividing into windows of size 15 and bins of two with a word size of 2
(binary).

#### Datasets

In this experiment, we consider two datasets:  (i) Iris[31]and  (ii)
GunPointfrom the UCR time series [30] archive  . In preprocessing the Iris
dataset, the versicolor class is removed, virginica is mapped to 0 and setosa
to 1. The input feature petal width is excluded, all samples are normalized to
unit norm, and then split into a 60/40 training/testing ratio. In the
preprocessing of the GunPoint dataset from the UCR time series archive [30],
we consider using a Bag-of-Words [32, 33, 34] approach to create a new
representation of the data with reduced dimensionality. Initially, the dataset
is divided into training and testing sets, from which only the first half of
each time series is selected for further processing. This truncation serves to
focus on a specific portion of the time series that is of interest. The Bag-
of-Words model is then applied to these truncated time series as shown in Fig.
3. The model operates by dividing each time series into windows of a specified
size (15 in this case) and encoding the data within each window into words
based on a predefined number of bins (2 bins here) and a word size of 2. The
resulting data is a binary vector of size 10, i.e. $x\in\\{0,1\\}^{10}$. This
representation of time series can also be made more expressive by increasing
the number of bins and the word size. Since it remains quantized by design, a
mapping to a basis state encoding can always be undertaken easily.

With the Bag-of-Words preprocessing, input space perturbations appear as bit-
flips of the binary data representation. Specifically, if a region of the data
snippet is perturbed sufficiently strong that its amplitudes falls into a
different bin, the corresponding bit is flipped. By constraining the number of
bits that are flipped, i.e., the Hamming weight distance, the magnitude of the
amplitude deviations induced by the perturbation can be controlled. Thus, we
conclude that Hamming weight constrained bit-flip noise is particularly
suitable for this type of basis state encoded data.

#### Network Architecture

In our experiments, we employ a QNN to function as an oracle in Grover’s
algorithm. This is achieved by incorporating a strongly entangling layer [35]
as the primary architecture component. Given that our datasets are simplified
to binary classification problems, we have trained the network on a singular
output qubit to perform the classification. In Table II, we report the
accuracy, number of qubits and layers for each dataset considered.

TABLE II: QNN architecture’s parameters. Dataset | Accuracy | Embedding | N. of Qubits | N. of Layers
---|---|---|---|---
Iris | 1.00 | Basis | 3 | 2
Gun Point | 0.92 | Basis | 10 | 50

#### QuAdRo Parameters

In our application of QuAdRo [15], we employed the iterative quantum phase
estimation [36] technique to decrease the number of ancillary qubits needed
for phase estimation. Given $M$ as maximum number of calls to the oracle, we
consider $m=\lceil\log_{2}(M)\rceil$. As a result, we execute the circuit $m$
times, with each run involving $\sqrt{N}$ shots, in contrast to the $N$-times
required by the standard implementation of randomized smoothing. For each
iteration $i$, we evaluated the $2^{i}$ power of Grover’s algorithm and
estimated the phase using a classical method. Employing all available $m$
qubits for phase estimation allows for executing the circuit with only
$\sqrt{N}$ repetitions. For determining the minimum value of $p(x)$, we adopt
a 90% confidence level ($\alpha=0.1$) using the Clopper-Pearson interval.

Figure 4: Certified accuracy of quantum randomized smoothing with
perturbations limited to a Hamming distance of 1 from the original input,
across different $\sigma$ values. Each evaluation maintained a consistent
number of 109 shots.

###  V-B Evaluation of 1-Hamming distant state distributions

To conduct a practical comparison of the probability outcomes derived from
various methods of creating states with $k$-Hamming distance, we opt to
measure the certified accuracy of a QNN that has been trained on the Iris
dataset. In Fig. 4, we showcase an assessment of quantum randomized smoothing,
utilizing a fixed number of trials (one billion shots), to examine how
certified accuracy varies with different values of $\sigma$ while maintaining
a 1-Hamming distance. This analysis highlights the impact of varying the
probability of a state flipping to another state that differs by a Hamming
distance of 1. Generally, higher values of $\sigma$ indicate a greater
tolerance to state changes, reflected in the certified accuracy across the
radius spectrum. While differences between amplitude encoding and our method
of state preparation are noticeable at larger radii, both methods generally
demonstrate similar behavior, emphasizing the commonalities in how states are
prepared.

TABLE III: Circuit preparation for 1-Hamming distant states. State | N. of qubits | N. of gates | Depth
---|---|---|---
Ampl. (Möttönen [28]) | 10 | 2045 | 2036
Our | 20 | 43 | 6

As indicated in Table III, there is a considerable difference in circuit
depths between the two methods: the Möttönen amplitude state preparation
technique results in a circuit depth of 2036, whereas our method for state
preparation yields a much shallower depth of just 6 when applied to the
GunPoint dataset, which features binary input with a dimensionality of 10.
Despite this, the quantum resource requirement in terms of qubits for our
approach is twice as high. In general, our method requires only
$\mathcal{O}(2n)$ gates for $n$ qubits, while Möttönen’s state preparation
scales exponentially at $\mathcal{O}(2^{n})$.

###  V-C Comparison between uniform and 1-Hamming distant state distributions.

(a) 1-Hamming distribution. (b) Uniform distribution.

Figure 5: Certified accuracy of quantum randomized smoothing between QuAdRo
and RS in terms of shots number for 1-Hamming distant and uniform
distributions. Both tests have been conducted with a $\sigma$ of 0.5 for the
Iris dataset.

Here, we compare the 1-Hamming distance distribution with respect to the
uniform distribution in terms of certified accuracy. For both scenarios, we
execute a plain version of randomized smoothing (RS) [11], which involves a
QNN placed ahead of the input distribution, running for a number of shots that
is quadratically higher in comparison to QuAdRo. Specifically, we run a set of
shots $\\{8^{2},10^{2},12^{2},50^{2}\\}$ for RS and $\\{8,10,12,50\\}$ for
QuAdRo. In the context of QuAdRo, we implement the iterative quantum phase
estimation which makes use of one single additional qubit. Consequently, the
total quantum circuit consists of 6 qubits for the uniform distribution and 9
qubits for the 1-Hamming distribution. We conduct our experiments in the
default.qubit simulator in PennyLane222PennyLane version 0.33.1
https://github.com/PennyLaneAI/pennylane..

In Fig. 5, we show a comparison for $\sigma=0.5$ for the Iris dataset. The
analysis shows that the 1-Hamming distribution has a tendency to reach a
higher level of certified accuracy when compared with the uniform
distribution. This result is attributed to the increased probability of
locating states in proximity to the initial state, as opposed to states
dispersed uniformly at a distance. In particular, for QuAdRo, achieving high
certified accuracy requires fewer than 50 shots of the circuit, compared to
144 shots for RS. With enough shots, both distributions achieve a certified
accuracy of 1. This outcome is because the Iris classifier achieves perfect
clean accuracy and with $\sigma=0.5$, as demonstrated in Fig. 2 and Fig. 1,
maintains the probability of the original input higher than that of the
neighbors.

###  V-D Certified robustness for time-series analysis

(a) 1-Hamming distribution. (b) Uniform distribution.

Figure 6: Certified accuracy of quantum randomized smoothing between QuAdRo
and RS in terms of number of shots for 1-Hamming distant and uniform
distributions. Both tests have been conducted with a $\sigma$ of 0.5 for the
GunPoint dataset.

In Fig. 6, we present a comparison based on the GunPoint dataset with
$\sigma=0.5$. As previously described, the input was transformed using a Bag-
of-Words model, resulting in a binary format to encode the state within the
basis. In addition, we consider the iterative quantum phase estimation, which
requires 13 qubits for a uniform distribution and 22 qubits for a 1-Hamming
distribution. In the 1-Hamming distribution case (Fig. 6.a), QuAdRo
consistently outperforms RS across all radii, with iterations of 20, 50, and
120 for QuAdRo compared against 50, 100, and 500 samples for RS. Within the
uniform distribution, we observe (Fig. 6.b) that QuAdRo demonstrates a
quadratic improvement when executed for 50 and 100 iterations compared to RS.

Due to the limited number of measurements, the phase estimation of QuAdRo does
not precisely fit the exact radius. This causes QuAdRo to exhibit smoother
performance compared to RS. The likelihood of accurately measuring the phase
up to $m$ bits is $\frac{8}{\pi^{2}}$. By conducting the experiment repeatedly
and applying the median estimate [37], we can rapidly achieve a 100% accuracy
rate. In practice, QuAdRo obtains probabilistic robustness certificates,
similar to RS, but with a quadratic speed-up in the number of shots required.
The average run-time for the simulated circuit with a 1-Hamming distribution
is roughly 21 minutes per sample on a server with a 128 CPU cores and 775 GB
of RAM memory.

##  VI Discussion of Experimental Results

The challenges of implementing the quantum phase estimation algorithm are
significant for NISQ devices. To address these challenges, our approach
includes several key considerations. First, using the Bag-of-Words model to
represent data helps reduce the dimensionality of time-series data while
retaining crucial information for accurate classification. Second, our
construction of the $k$-Hamming distance for state preparation enables a
relatively shallow circuit. Although this circuit approximates the exact
distribution, the depth is three orders of magnitude smaller than that
required for amplitude state preparation. This reduction is crucial for the
overall preparation process of the quantum phase estimation. Third, the use of
iterative quantum phase estimation decreases the number of qubits needed for
precise phase estimation. This makes the QuAdRo algorithm more compatible with
current quantum hardware, which is restricted by qubit stability and coherence
times.

Experimental findings demonstrate the efficacy of the 1-Hamming distance
approach in achieving higher certified accuracy, especially in scenarios with
limited shots. This result highlights the advantages of using constrained
perturbation models like the 1-Hamming distance in situations where it is
reasonable to consider bit-flips as principal source of noise. However, it is
crucial to recognize that other sources of noise can significantly impact the
practical implementation and benefits of robustness certification. Gate noise
and decoherence, in particular, pose additional challenges to the overall
certification process, potentially affecting the algorithm’s performance on
real quantum hardware. To address these broader noise concerns and enhance the
QuAdro algorithm’s performance on NISQ devices, a comprehensive strategy is
essential. Error mitigation techniques like zero-noise extrapolation can be
applied to reduce noise effects on circuit outputs. Simultaneously, minimizing
two-qubit gates through circuit layout optimization can enhance quantum
operation fidelity. This combined approach aims to boost the algorithm’s
resilience to NISQ device noise, potentially improving its real-world
performance and accuracy.

##  VII Conclusion

In this study, we establish a connection between input transformation and data
encoding. Our analysis indicates overlooked opportunities in encoding
techniques and the related noise types. Expanding on previous methods of
preparing input distributions, we propose a new encoding strategy that
generates $k$-Hamming distant states using fewer gates than needed for
amplitude embeddings. Our method, however, does not exactly replicate the
distribution of $k$-Hamming distant states and requires $n$ additional
ancillary qubits, where $n$ is the number of qubits used for embedding the
state. Despite this drawback, the complexity of our circuit is significantly
reduced, being three orders of magnitude less deep than that required for
amplitude state preparation. In our empirical validation, preparing
$k$-Hamming states instead of uniform states resulted in higher certified
accuracy for the same number of shots. In addition, using iterative quantum
phase estimation demonstrates a quadratic improvement in speed compared to
standard sampling in verifying the robustness of QNN.

## References

  * [1] Peter Wittek  “Quantum machine learning: what quantum computing means to data mining”  Academic Press, 2014
  * [2] Maria Schuld, Ilya Sinayskiy and Francesco Petruccione  “An introduction to quantum machine learning”  In _Contemporary Physics_ 56.2 Taylor & Francis, 2015, pp. 172–185
  * [3] Jacob Biamonte et al.  “Quantum machine learning”  In _Nature_ 549.7671 Nature Publishing Group UK London, 2017, pp. 195–202
  * [4] Amira Abbas et al.  “The power of quantum neural networks”  In _Nature Computational Science_ 1.6 Nature Publishing Group US New York, 2021, pp. 403–409
  * [5] M Cerezo et al.  “Challenges and opportunities in quantum machine learning”  In _Nature Computational Science_ 2.9 Nature Publishing Group US New York, 2022, pp. 567–576
  * [6] Sirui Lu, Lu-Ming Duan and Dong-Ling Deng  “Quantum adversarial machine learning”  In _Physical Review Research_ 2.3 APS, 2020, pp. 033212
  * [7] Maxwell T West et al.  “Benchmarking adversarially robust quantum machine learning at scale”  In _Physical Review Research_ 5.2 APS, 2023, pp. 023186
  * [8] Nicola Franco et al.  “Predominant Aspects on Security for Quantum Machine Learning: Literature Review”  In _arXiv preprint arXiv:2401.07774_ , 2024
  * [9] Battista Biggio et al.  “Evasion attacks against machine learning at test time”  In _Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2013, Prague, Czech Republic, September 23-27, 2013, Proceedings, Part III 13_ , 2013, pp. 387–402  Springer
  * [10] Christian Szegedy et al.  “Intriguing properties of neural networks”  In _2nd International Conference on Learning Representations, ICLR 2014_ , 2014
  * [11] Jeremy Cohen, Elan Rosenfeld and Zico Kolter  “Certified adversarial robustness via randomized smoothing”  In _international conference on machine learning_ , 2019, pp. 1310–1320  PMLR
  * [12] Hadi Salman et al.  “Provably robust deep learning via adversarially trained smoothed classifiers”  In _Advances in Neural Information Processing Systems_ 32, 2019
  * [13] Maurice Weber et al.  “Optimal provable robustness of quantum classification via quantum hypothesis testing”  In _npj Quantum Information_ 7.1 Nature Publishing Group UK London, 2021, pp. 76
  * [14] Jhih-Cing Huang et al.  “Certified Robustness of Quantum Classifiers Against Adversarial Examples Through Quantum Noise”  In _ICASSP 2023 - 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2023, pp. 1–5  DOI: 10.1109/ICASSP49357.2023.10095030
  * [15] Aditya Sahdev and Mehul Kumar  “Adversarial Robustness based on Randomized Smoothing in Quantum Machine Learning”, 2023  URL: https://openreview.net/forum?id=o-Yxq5iicIp
  * [16] Ji Guan, Wang Fang and Mingsheng Ying  “Robustness verification of quantum classifiers”  In _Computer Aided Verification: 33rd International Conference, CAV 2021, Virtual Event, July 20–23, 2021, Proceedings, Part I 33_ , 2021, pp. 151–174  Springer
  * [17] Ji Guan, Wang Fang and Mingsheng Ying  “Verifying Fairness in Quantum Machine Learning”  In _International Conference on Computer Aided Verification_ , 2022, pp. 408–429  Springer
  * [18] Nicola Franco et al.  “Quantum robustness verification: A hybrid quantum-classical neural network certification algorithm”  In _2022 IEEE International Conference on Quantum Computing and Engineering (QCE)_ , 2022, pp. 142–153  IEEE
  * [19] Nicola Franco et al.  “Efficient MILP Decomposition in Quantum Computing for ReLU Network Robustness”  In _2023 IEEE International Conference on Quantum Computing and Engineering (QCE)_ 01, 2023, pp. 524–534  DOI: 10.1109/QCE57702.2023.00066
  * [20] Yuxuan Du et al.  “Quantum noise protects quantum classifiers against adversaries”  In _Phys. Rev. Res._ 3 American Physical Society, 2021, pp. 023153  DOI: 10.1103/PhysRevResearch.3.023153
  * [21] Vincent Tjeng, Kai Y Xiao and Russ Tedrake  “Evaluating Robustness of Neural Networks with Mixed Integer Programming”  In _International Conference on Learning Representations_ , 2018
  * [22] Guy Katz et al.  “Reluplex: An efficient SMT solver for verifying deep neural networks”  In _Computer Aided Verification: 29th International Conference, CAV 2017, Heidelberg, Germany, July 24-28, 2017, Proceedings, Part I 30_ , 2017, pp. 97–117  Springer
  * [23] Sven Gowal et al.  “On the Effectiveness of Interval Bound Propagation for Training Verifiably Robust Models”, 2019  arXiv:1810.12715 [cs.LG]
  * [24] Asaf Gendler, Tsui-Wei Weng, Luca Daniel and Yaniv Romano  “Adversarially robust conformal prediction”  In _International Conference on Learning Representations_ , 2021
  * [25] Ge Yan, Yaniv Romano and Tsui-Wei Weng  “Provably Robust Conformal Prediction with Improved Efficiency”  In _The Twelfth International Conference on Learning Representations_ , 2023
  * [26] Karol Życzkowski and Hans-Jürgen Sommers  “Average fidelity between random quantum states”  In _Phys. Rev. A_ 71 American Physical Society, 2005, pp. 032313  DOI: 10.1103/PhysRevA.71.032313
  * [27] Andreas Bärtschi and Stephan Eidenbenz  “Deterministic preparation of Dicke states”  In _International Symposium on Fundamentals of Computation Theory_ , 2019, pp. 126–139  Springer
  * [28] Mikko Möttönen, Juha J Vartiainen, Ville Bergholm and Martti M Salomaa  “Transformation of quantum states using uniformly controlled rotations”  In _Quantum Information & Computation_ 5.6 Rinton Press, Incorporated Paramus, NJ, 2005, pp. 467–473
  * [29] Guang-He Lee, Yang Yuan, Shiyu Chang and Tommi Jaakkola  “Tight certificates of adversarial robustness for randomly smoothed classifiers”  In _Advances in Neural Information Processing Systems_ 32, 2019
  * [30] Yanping Chen et al.  “The UCR Time Series Classification Archive” www.cs.ucr.edu/~eamonn/time_series_data/, 2015
  * [31] R.. Fisher  “Iris” DOI: https://doi.org/10.24432/C56C76, UCI Machine Learning Repository, 1988
  * [32] Kilian Weinberger et al.  “Feature hashing for large scale multitask learning”  In _Proceedings of the 26th annual international conference on machine learning_ , 2009, pp. 1113–1120
  * [33] Jessica Lin, Rohan Khade and Yuan Li  “Rotation-invariant similarity in time series using bag-of-patterns representation”  In _Journal of Intelligent Information Systems_ 39 Springer, 2012, pp. 287–315
  * [34] Johann Faouzi and Hicham Janati  “pyts: A Python Package for Time Series Classification”  In _Journal of Machine Learning Research_ 21.46, 2020, pp. 1–6  URL: http://jmlr.org/papers/v21/19-763.html
  * [35] Maria Schuld, Alex Bocharov, Krysta M Svore and Nathan Wiebe  “Circuit-centric quantum classifiers”  In _Physical Review A_ 101.3 APS, 2020, pp. 032308
  * [36] Miroslav Dobšıček, Göran Johansson, Vitaly Shumeiko and Göran Wendin  “Arbitrary accuracy iterative quantum phase estimation algorithm using a single ancillary qubit: A two-qubit benchmark”  In _Physical Review A_ 76.3 APS, 2007, pp. 030306
  * [37] Koichi Miyamoto  “Bermudan option pricing by quantum amplitude estimation and Chebyshev interpolation”  In _EPJ Quantum Technology_ 9.1 Springer, 2022, pp. 1–27


# Computing the Newton-step faster than Hessian accumulation

Akshay Srinivasan  Emanuel Todorov

###### Abstract

Computing the Newton-step of a generic function with $N$ decision variables
takes $O(N^{3})$ flops. In this paper, we show that given the computational
graph of the function, this bound can be reduced to $O({m\tau}^{3})$, where
$\tau,m$ are the width & size of a tree-decomposition of the graph. The
proposed algorithm generalizes nonlinear optimal-control methods based on LQR
to general optimization problems and provides non-trivial gains in iteration-
complexity even in cases where the Hessian is dense.

Machine Learning, ICML

##  1 Introduction

Newton’s method forms the basis for second-order methods in optimization.
Computing the Hessian of a generic function
${\mathpzc{f}}:\mathbb{R}^{N}\rightarrow\mathbb{R}$, requires $O(N^{2})$
flops; inverting this matrix requires a further $O(N^{3})$ flops. This super-
linear scaling in compute/memory requirements is a major obstacle in the
application of such methods despite their quadratic convergence.

The iteration-complexity of second-order methods has necessitated the
development of specialized algorithms for restricted classes of problems. In
particular, _Differential Dynamic Programming_ (DDP) methods have been
proposed for solving nonlinear optimal-control problems, and these methods,
quite miraculously, achieve quadratic convergence despite their linear
iteration-complexities (Jacobson & Mayne, 1970) (De O. Pantoja, 1988) (Wright,
1990).

DDP-methods were historically derived using _Dynamic Programming_ (DP) and
resemble LQR control-design. They were later, quite surprisingly, also found
to be related to both, the Newton-iteration on the unconstrained problem (De
O. Pantoja, 1988), and the SQP-iteration on the constrained problem (Wright,
1990). Sadly, the method is not trivially generalized to other objective
functions.

Applications of DP to other domains, esp. the solution of least-squares
problems and inference in graphical models, have seen great success, but such
techniques are not applicable to problems in optimal-control or neural-
networks, due to their inability to handle function compositions.

In this paper, we develop a method that exploits the compositional structure
in a given objective function in order to automatically derive a _fast_ Newton
update. This is done by extending the connections of DDP to constrained &
unconstrained optimization using tools from _Automatic Differentiation_ (AD),
and by using techinques from _Sparse Linear Algebra_ (SLA) to bound iteration-
complexity.

##  2 Problem setup

###  2.1 Computational graph

Let $\mathcal{G}$ be a Directed Acyclic Graph (DAG). Define,

$\displaystyle\operatorname{\operatorname{pa}}(u)$
$\displaystyle\triangleq\\{v|(v,u)\in E[\mathcal{G}]\\},\quad\mbox{(parents)}$
(1) $\displaystyle\operatorname{\operatorname{ch}}(u)$
$\displaystyle\triangleq\\{v|(u,v)\in
E[\mathcal{G}]\\}.\quad\mbox{(children)}$

Let every vertex $v\in V[\mathcal{G}]$ be associated with a _state_
$\mathtt{X}_{v}\in U_{v}\subset\mathbb{R}^{n_{v}}$ for some open set $U_{v}$.
Let $\mathtt{X}_{A}$ be the (labelled) concatenation of _state_ s associated
with vertices in set $A\subset V[\mathcal{G}]$.

Let the _input_ nodes
$\operatorname{\mathtt{Input}}\triangleq\\{u_{1},u_{2},\dots,u_{n}\\}\subset
V[\mathcal{G}]$ be the parentless vertices in $\mathcal{G}$. Let the _states_
of the non-input nodes be defined recursively by
$\mathtt{X}_{v}\triangleq\Phi_{v}(\mathtt{X}_{\operatorname{\operatorname{pa}}(v)})$
for some given function
$\Phi_{v}:\prod_{z\in\operatorname{\operatorname{pa}}(v)}U_{z}\rightarrow
U_{v}$. Since $\mathcal{G}$ is a DAG, it follows that
$\mathtt{X}_{V[\mathcal{G}]}$ is uniquely determined from the input state
$\mathtt{X}_{\operatorname{\mathtt{Input}}}$ and functions
$\\{\Phi_{v}\\}_{v\in V[\mathcal{G}]\backslash X}$.

An objective function ${\mathpzc{f}}:U_{x_{1}}\times\dots
U_{x_{n}}\rightarrow\mathbb{R}$ has the _computational structure_ given by the
tuple $(\mathcal{G},\\{\Phi_{v}\\}_{v\in
V[\mathcal{G}]\backslash\operatorname{\mathtt{Input}}},\\{{\mathpzc{l}}_{v}\\}_{v\in
V[\mathcal{G}]})$ if it can be written as the sum of local objectives
${\mathpzc{l}}_{v}:\prod_{z\in\\{v\\}\cup\operatorname{\operatorname{pa}}(v)}U_{z}\rightarrow\mathbb{R}$
on the graph $\mathcal{G}$ in the following form (2),

$\displaystyle{\mathpzc{f}}:(\mathtt{X}_{x_{1}},\dots,\mathtt{X}_{x_{n}})$
$\displaystyle\mapsto\sum_{v\in
V[\mathcal{G}]}{\mathpzc{l}}_{v}(\mathtt{X}_{v\cup\operatorname{\operatorname{pa}}(v)}),$
(2)
$\displaystyle\mathtt{X}_{v}\leftarrow\Phi_{v}(\mathtt{X}_{\operatorname{\operatorname{pa}}(v)})$
$\displaystyle,\quad\forall v\in
V[\mathcal{G}],\operatorname{\operatorname{pa}}(v)\neq\emptyset.$

The _assignment_ operator, ’$\leftarrow$’, is explicitly distinguished from
the equality operator, ’$=$’, which is taken to represent a constraint in the
program. We refer to the DAG $\mathcal{G}$ as the _computational graph_ of
$f(\cdot)$.

Notation. The symbolism
$\partial_{u}v\triangleq{\partial\mathtt{X}_{v}\over\partial\mathtt{X}_{u}}\big{|}_{\mathtt{X}_{\operatorname{\mathtt{Input}}}}$,
will be employed for succintly denoting partial derivatives. The derivative
operator _w.r.t_ the (labelled) set $A=\\{v_{1},v_{2},\dots\\}\subset
V[\mathcal{G}]$, will similarly be denoted by
$\partial_{A}\triangleq[\partial_{a_{1}},\partial_{a_{2}},\dots]$.

$u_{0}$$u_{1}$$u_{2}$$u_{3}$$u_{n-1}$$x_{1}$$x_{2}$$x_{3}$$x_{n-1}$$x_{n}$
Figure 1: (Optimal control) The nodes $\\{x_{i}\\},\\{u_{i}\\}$ represent the
states and control inputs of the dynamical system.

###  2.2 Motivation.

Consider the canonical optimal-control problem,

$\displaystyle\min_{u_{0},u_{1},\dots,u_{n-1}}\left[{\mathpzc{J}}(u_{0},\dots,u_{n})\triangleq\sum_{i=0}^{n-1}{\mathpzc{l}}_{i}(x_{i},u_{i})+{\mathpzc{l}}_{n}(x_{n})\right],$
(3) $\displaystyle\quad\forall i,x_{i+1}\leftarrow\mathtt{f}(x_{i},u_{i}),$

wherein, $\mathtt{f}(\cdot,\cdot)$ denotes the system dynamics, and
${\mathpzc{l}}_{i}(\cdot,\cdot)$’s denote the losses incurred. The order of
computation for the objective (3) can be represented by a linear-chain as
shown in Figure 1.

The computational graph Figure 1, while sparse, does not imply sparsity in the
objective. Substitution of the symbol assignments in the unconstrained
objective, has a cascading effect whereby variables get _transport_ ed down
$\mathcal{G}$,

$\displaystyle{\mathpzc{J}}(u_{0},\dots,u_{n})=$
$\displaystyle{\mathpzc{l}}_{0}(x_{0},u_{0})+$
$\displaystyle{\mathpzc{l}}_{1}(\mathtt{f}(x_{0},u_{0}),u_{1})+$
$\displaystyle{\mathpzc{l}}_{2}(\mathtt{f}(\mathtt{f}(x_{0},u_{0}),u_{1}),u_{2})+\dots.$

This results in a Hessian that is both fully dense and (very) badly
conditioned. These two properties render both iterative and direct-
factorization methods unsuitable for solving this problem.

The computational graph does, however, fully encode the sparsity of the
Lagrangian of the constrained program; the constrained form being obtained by
replacing ’$\leftarrow$’ with ’$=$’, and including $\\{x_{i}\\}^{\prime}s$ in
the domain of optimization. This results in a _fast_ linear-time SQP
(Sequential Quadratic Programming) iteration for the optimal-control problem
but increases implementation complexity by requiring _dual_ and _non-input_
updates.

This presents one with a strange set of choices: slow-simple unconstrained
optimization or fast-complex constrained optimization. DDP-like methods
resolve this dichotomy for optimal-control problems (Jacobson & Mayne, 1970)
(De O. Pantoja, 1988) (Wright, 1990). Our goal is to do the same for general
problems.

##  3 Graphical Newton

Consider the objective function in (2), defined by the tuple
$(\mathcal{G},\\{\Phi_{v}\\},\\{{\mathpzc{l}}_{v}\\})$. The optimization
problem of interest is the following,

$\displaystyle\min_{\\{\mathtt{X}_{v}|\forall
v\in\operatorname{\mathtt{Input}}\\}}$
$\displaystyle\left({\mathpzc{f}}\triangleq\sum_{v\in
V[\mathcal{G}]}{\mathpzc{l}}_{v}(\mathtt{X}_{v\cup\operatorname{\operatorname{pa}}(v)})\right),$
(4)
$\displaystyle\mathtt{X}_{v}\leftarrow\Phi_{v}(\mathtt{X}_{\operatorname{\operatorname{pa}}(v)}),$
$\displaystyle\quad\forall v\in
V[\mathcal{G}],\operatorname{\operatorname{pa}}(v)\neq\emptyset.$

The corresponding constrained problem is obtained by simply replacing the
operator ’$\leftarrow$’ with ’$=$’ in (4) and enlarging the domain of
optimization.

$\displaystyle\min_{\\{\mathtt{X}_{v}|\forall v\in V[\mathcal{G}]\\}}$
$\displaystyle\left({\mathpzc{f}}\triangleq\sum_{v\in
V[\mathcal{G}]}{\mathpzc{l}}_{v}(\mathtt{X}_{v\cup\operatorname{\operatorname{pa}}(v)})\right),$
(5)
$\displaystyle\mathtt{X}_{v}=\Phi_{v}(\mathtt{X}_{\operatorname{\operatorname{pa}}(v)}),$
$\displaystyle\quad\forall v\in
V[\mathcal{G}],\operatorname{\operatorname{pa}}(v)\neq\emptyset.$

The constrained program (5) induces the following Lagrangian function,

$\displaystyle\mathcal{L}(\mathtt{X}_{V[\mathcal{G}]},\lambda_{V[\mathcal{G}]\backslash\operatorname{\mathtt{Input}}})\triangleq\sum_{v\in
V[\mathcal{G}]}{\mathpzc{l}}_{v}(\mathtt{X}_{v\cup\operatorname{\operatorname{pa}}(v)})+\sum_{\begin{subarray}{c}v\in
V[\mathcal{G}],\\\
\operatorname{\operatorname{pa}}(v)\neq\emptyset\end{subarray}}\lambda_{v}^{\mathtt{T}}\mathbf{h}_{v}(\mathtt{X}_{v\cup\operatorname{\operatorname{pa}}(v)}),$
(6)
$\displaystyle\mbox{where,}\;\mathbf{h}_{v}(\mathtt{X}_{v\cup\operatorname{\operatorname{pa}}(v)})\triangleq\Phi_{v}(\mathtt{X}_{\operatorname{\operatorname{pa}}(v)})-\mathtt{X}_{v},\;\forall
v\in V[\mathcal{G}],\operatorname{\operatorname{pa}}(v)\neq\emptyset.$

The necessary first-order conditions for optimality are given by,

$\partial_{V}\mathcal{L}(\mathtt{X}_{V}^{*},\lambda_{V}^{*})=0,\quad
h(\mathtt{X}_{V}^{*})=0.$ (7)

Linearization of the first-order conditions for this constrained problem
(around a nominal $(\tilde{\mathtt{X}_{V}},\tilde{\lambda})$) yields a system
of KKT equations, whose solution yields the SQP search direction.

$\left[\begin{array}[]{c
c}\partial^{2}_{V}\mathcal{L}&\partial_{V}\mathbf{h}^{\mathtt{T}}\\\
\partial_{V}\mathbf{h}&0\end{array}\right]\left[\begin{array}[]{c}\delta\mathtt{X}_{V}\\\
\lambda^{+}\end{array}\right]=\left[\begin{array}[]{c}-\partial_{V}{\mathpzc{f}}\\\
-\mathbf{h}\end{array}\right],$ (8)

where $\lambda^{+}\triangleq\tilde{\lambda}+\delta\lambda$. The sequence of
iterates obtained by taking appropriate steps along
$(\delta\mathtt{X}_{V},\delta\lambda)$, converges quadratically near a
strongly-convex local minimum.

The principal result of this paper is the connection between the unconstrained
and constrained formulations,

######  Theorem 1 (Newton direction)

The Newton direction for unconstrained problem (4) is given by an iteration of
SQP for the constrained problem (5), when $\mathtt{X}_{V[\mathcal{G}]}$ is
feasible and when $\forall v,\lambda_{v}=\partial_{v}{\mathpzc{f}}$ (12).

_Proof._ See Appendix A.

1 can trivially be extended to arbitrary objective functions on the graph
$\mathcal{G}$ and is not restricted to form (5) where objective functions have
the same sparsity as the underlying computation. The dual values
$\lambda_{v}=\partial_{v}{\mathpzc{f}}$ can be computed in linear-time with
reverse-mode AD.

1 immediately leads to the sparse Newton-iteration given in Algorithm 1.

Algorithm 1 Graphical Newton

1: Input: initial $\mathtt{X}_{\operatorname{\mathtt{Input}}}^{0}$, tuple
$(\mathcal{G},\\{\Phi_{v}\\},\\{{\mathpzc{l}}_{v}\\})$

2: repeat

3: Compute non-inputs, local objectives, and their derivatives (eq. 4).

4: Set dual-values to $\lambda_{v}=\partial_{v}{\mathpzc{f}},\forall v$ by
reverse-mode AD (eq. 12).

5: Solve the KKT-system (eq. 8).

6: Compute step-length $\eta$ _via_ linesearch on inputs
$\mathtt{X}_{\operatorname{\mathtt{Input}}}$ (only).

7: Update inputs (only):
$\mathtt{X}_{\operatorname{\mathtt{Input}}}\leftarrow\mathtt{X}_{\operatorname{\mathtt{Input}}}+\eta\delta\mathtt{X}_{\operatorname{\mathtt{Input}}}$

8: until
$||{\partial_{\operatorname{\mathtt{Input}}}{\mathpzc{f}}}||\leq\epsilon$

Algorithm 1 differs markedly from both SQP and current techniques in AD. The
non-input primals are computed by running a forward pass on the computation
graph, while the dual-values are set to fixed values by running reverse-mode
AD. Once the KKT system is solved, the algorithm proceeds to update the primal
inputs only, without requiring non-input updates, dual updates, contraint
penalties, or any of the other machinery from constrained optimization.

The efficiency of Algorithm 1 stems from the fact that the KKT matrix is
typically much sparser than the Hessian. It can be further shown that the
common Hessian-accumulation/inversion method is equivalent to a particular,
generally a (very) suboptimal, pivot-ordering for solving the sparse KKT
system.

###  3.1 KKT complexity

The run-time of the Algorithm 1 depends crucially on the time taken to solve
the sparse KKT system (8) at every iteration. The complexity of solving such
sparse systems depends inturn on the support-graph of the underlying matrix.

It has been observed that many problems defined locally on graphs, can be
solved in linear-time on trees using dynamic programming. These techniques can
be extended to general graphs by grouping vertices in such a way as to _mimic_
a tree. The size of the largest cluster in this _tree-decomposition_ \- termed
_tree-width_ \- governs the dominant term in the time-complexity of the
resulting overall algorithm. Computing the tree-decomposition with minimal
tree-width is NP-hard, but heuristics for finding elimination orderings are
often known to do well in practice.

Run-time complexities in terms of the tree-width are well-established for
closely related problems such as Cholesky decomposition, but they appear to be
unknown for structured KKT systems such as (8). The following theorem
establishes the required bound for structured-KKT systems arising in Algorithm
1.

######  Theorem 2 (KKT tree-width)

The KKT system (8) associated with the constrained problem (5) can be solved
in time $O(m\operatorname{tw}(\mathcal{G})^{3})$, given the tree-
decomposition.

_Proof._ See Appendix B.

##  4 Special cases

###  4.1 Optimal control

$u_{0}$$u_{1}$$u_{2}$$u_{3}$$u_{n-1}$$x_{1}$$x_{2}$$x_{3}$$x_{n-1}$$x_{n}$

$u_{0}$$u_{1}$$u_{2}$$u_{3}$$u_{n-1}$$x_{1}$$x_{2}$$x_{3}$$x_{n-1}$$x_{n}$

Figure 2: Structure of the optimal control problem as defined in (9). Top:
Computational graph. Bottom: Constraint graph.

Consider, again, the canonical optimal control problem from (3),

$\displaystyle\min_{u_{0},u_{1},\dots,u_{n-1}}\left[{\mathpzc{J}}(u_{0},\dots,u_{n-1})\triangleq\sum_{i=0}^{n-1}{\mathpzc{l}}_{i}(x_{i},u_{i})+{\mathpzc{l}}_{n}(x_{n})\right],$
(9) $\displaystyle\quad\forall i,x_{i+1}\leftarrow\mathtt{f}(x_{i},u_{i}).$

The computational graph and its moralized relative for this problem are shown
in Figure 2. The constraint graph is chordal, and permits multiple optimal
elimination orderings.

Dynamic Programming.

Table 1: DDP methods and their key differences. Method | $\lambda_{i}$ | _back-substitution_
---|---|---
DDP (Jacobson & Mayne, 1970) | $\partial_{x_{i+1}}V_{i+1}\cdot\partial_{x_{i}}{\mathtt{f}}$ | Non-Linear (11)
Stagewise-Newton (De O. Pantoja, 1988) | $\partial_{x_{i}}{\mathpzc{J}}$ | Linear
Nonlinear Stagewise-Newton (Liao & Shoemaker, 1992) | $\partial_{x_{i}}{\mathpzc{J}}$ | Non-Linear (11)
iLQR/iLQG (Li & Todorov, ) | 0 | Non-Linear (11)

Let’s consider the LQR order, $x_{n},u_{n-1},x_{n-1}\dots,u_{0}$. The
principal minor of the KKT matrix, corresponding to the clique
$\\{x_{n},x_{n-1},u_{n-1},\lambda_{n-1}^{+}\\}$ is given by,

$\displaystyle\begin{array}[]{l c c c | c}\partial^{2}_{x}V_{n}&&&-\operatorname{\mathtt{I}}^{\operatorname{\mathtt{T}}}&-\partial_{x}V_{n}\\\ &\partial_{x}^{2}\mathcal{H}_{n-1}&\partial_{xu}^{2}\mathcal{H}_{n-1}&\partial_{x}\mathtt{f}&-\partial_{x}{\mathpzc{l}}_{n-1}\\\ &\partial_{ux}^{2}\mathcal{H}_{n-1}&\partial_{uu}^{2}\mathcal{H}_{n-1}&\partial_{u}\mathtt{f}&-\partial_{u}{\mathpzc{l}}_{n-1}\\\ -\operatorname{\mathtt{I}}^{\operatorname{\mathtt{T}}}&\partial_{x}\mathtt{f}&\partial_{u}\mathtt{f}&0&0\end{array},$ (10) $\displaystyle\mbox{where,}V_{n}=l_{n},\mathcal{H}_{n-1}(x,u)={\mathpzc{l}}_{n-1}(x,u)+\lambda_{n-1}\cdot\mathtt{f}(x,u).$

Note that we can only eliminate variables $\\{x_{n},u_{n-1}\\}$ in the above,
since $x_{n-1}$ is part of the separator _i.e_ it is connected to nodes
outside the clique.
Eliminating $\delta
x_{n}={\partial^{2}V_{n}}^{-1}(-\partial_{x}V_{n}+\lambda_{n-1}^{+})$,

$\begin{array}[]{c c c | c}\partial_{x}^{2}\mathcal{H}_{n-1}&\partial_{xu}^{2}\mathcal{H}_{n-1}&\partial_{x}\mathtt{f}^{\operatorname{\mathtt{T}}}&-\partial_{x}{\mathpzc{l}}_{n-1}\\\ \partial_{ux}^{2}\mathcal{H}_{n-1}&\partial_{uu}^{2}\mathcal{H}_{n-1}&\partial_{u}\mathtt{f}^{\operatorname{\mathtt{T}}}&-\partial_{u}{\mathpzc{l}}_{n-1}\\\ \partial_{x}\mathtt{f}&\partial_{u}\mathtt{f}&-{\partial^{2}V_{n}}^{-1}&-{\partial^{2}V_{n}}^{-1}\partial V_{n}\end{array},$

Eliminating
$\lambda^{+}_{n-1}=-\partial^{2}V_{n}(-{\partial^{2}V_{n}}^{-1}\partial
V_{n}-\partial_{x}\mathtt{f}\delta x_{n-1}-\partial_{u}\mathtt{f}_{u}\delta
u_{n-1})$,

$\begin{array}[]{c c | c}Q_{xx}&Q_{ux}^{\operatorname{\mathtt{T}}}&-q_{x}\\\ Q_{ux}&Q_{uu}&-q_{u}\end{array},$

where,

$\begin{array}[]{l}Q_{ab}=\partial_{ab}^{2}\mathcal{H}_{n-1}+\partial_{a}\mathtt{f}\;\partial^{2}V_{n}\;\partial_{b}\mathtt{f}^{\operatorname{\mathtt{T}}},\\\
q_{a}=\partial_{a}{\mathpzc{l}}_{n-1}+\partial_{a}\mathtt{f}^{\operatorname{\mathtt{T}}}\;\partial
V_{n}.\end{array}$

Eliminate $\delta u_{n-1}={Q_{uu}}^{-1}(-Q_{u}-Q_{ux}\delta x_{n-1})$,

$\begin{array}[]{c | c}Q_{xx}+Q_{ux}^{\operatorname{\mathtt{T}}}{Q_{uu}}^{-1}Q_{ux}&-(Q_{x}+Q_{ux}^{T}{Q_{uu}}^{-1}Q_{u})\end{array},$

Re-write the above as,

$\begin{array}[]{c | c}\partial_{x}^{2}V_{n-1}&-\partial_{x}^{2}V_{n-1}.\end{array},$

The eliminations of the adjoining clique can be carried out in a similar
manner.

This procedure is identical to the _backward pass_ in DDP-methods.111These
computations can be carried out more efficiently & stably in the square-root
form. In the backsubstitution phase, also known as _forward pass_ , the non-
linearities can be used directly without having to resort to the use of
linearizations,

$\displaystyle\delta x_{n}$
$\displaystyle={\partial^{2}V_{n}}^{-1}(-\partial_{x}V_{n}+\lambda_{n-1}^{+})$
(11)
$\displaystyle={\partial^{2}V_{n}}^{-1}(-\partial_{x}V_{n}-\partial^{2}V_{n}(-{\partial^{2}V_{n}}^{-1}\partial
V_{n}$ $\displaystyle\hskip 60.00009pt-\partial_{x}\mathtt{f}\delta
x_{n-1}-\partial_{u}\mathtt{f}_{u}\delta u_{n-1}))$
$\displaystyle=\partial_{x}\mathtt{f}\delta
x_{n-1}+\partial_{u}\mathtt{f}_{u}\delta u_{n-1}.$
$\displaystyle\approx\mathtt{f}(x_{n-1}+\delta x_{n-1},u_{n-1}+\delta
u_{n-1})-\mathtt{f}(x_{n-1},u_{n-1}).$

The Table 1 illustrates how variations in _forward_ & _backward_ passes,
correspond to known DP algorithms. The theory also makes it trivial to
generalize these algorithms to higher-order dynamics and constrained problems
(Srinivasan & Todorov, 2015).

##  5 Discussion

The method presented in this paper can be used to compute the Newton step in
time $O(m\operatorname{tw}^{3})$, where $\operatorname{tw},m,$ are the width
and size of a tree-decomposition of the computation graph. The method
generalizes many specialized DDP algorithms in numerical optimal-control.

However, while the work presented here provides, in a sense, optimal
iteration-complexities, many real-world problems in machine-learning also have
large tree-widths. The KKT-matrix factorization of such problems is also quite
rife with redundant computation, and indicates the necessity for partial
symbolic-condensation in sparse LDL (and QP) solvers. These are topics for
future work.

## References

  * De O. Pantoja (1988) De O. Pantoja, J.  Differential Dynamic Programming and Newton’s Method.  _International Journal of Control_ , 47(5):1539–1553, 1988.
  * Jacobson & Mayne (1970) Jacobson, D. H. and Mayne, D. Q.  _Differential Dynamic Programming_.  North-Holland, 1970.
  * (3) Li, W. and Todorov, E.  Iterative Linear Quadratic Regulator Design for Nonlinear Biological Movement Systems.
  * Liao & Shoemaker (1992) Liao, L.-z. and Shoemaker, C. A.  Advantages of Differential Dynamic Programming Over Newton’s Method for Discrete-Time Optimal Control Problems.  Technical report, Cornell University, 1992.
  * Srinivasan & Todorov (2015) Srinivasan, A. and Todorov, E.  Graphical Newton.  _arXiv preprint arXiv:1508.00952_ , 2015.
  * Wright (1990) Wright, S. J.  Solution of Discrete-Time Optimal Control Problems on Parallel Computers.  _Parallel Computing_ , 16(2):221–237, 1990.

##  Appendix A Graphical Newton

This section presents a proof of 1. We proceed by first recalling the
$1^{st},2^{nd}$-order relations from AD for the objective (4); this is then
related to the problem of computing the solutiong to the KKT system in (8)
thus completing the proof.

###  A.1 Derivative relations on $\mathcal{G}$

Reverse-mode AD. The first derivative of the objective ${\mathpzc{f}}(\cdot)$
(4) can be calculated by applying the chain rule over $\mathcal{G}$,

$\displaystyle\forall v,\quad\partial_{v}{\mathpzc{f}}=\sum_{s\in
v\cup\operatorname{\operatorname{ch}}(v)}\partial_{v}{\mathpzc{l}}_{s}+\sum_{d\in\operatorname{\operatorname{ch}}(v)}\partial_{d}{\mathpzc{f}}^{\mathtt{T}}\;\partial_{v}d;$
(12) $\displaystyle\quad\quad\quad
v\in\operatorname{\operatorname{pa}}(d)\Rightarrow\partial_{v}d\triangleq{\partial\Phi_{d}(\mathtt{X}_{\operatorname{\operatorname{pa}}(d)})\over\partial\mathtt{X}_{v}}.$

Since $\mathcal{G}$ is a DAG, there exist child-less nodes (_i.e_
$\operatorname{\operatorname{ch}}(v)=\emptyset$) from which the recursion can
be initialized. The recursion then proceeds backward on $\mathcal{G}$ in a
breadth-first order. This algorithm is known in literature as reverse-mode AD.

Hessian-vector AD: For a given infinitesimal change
$\delta{\mathtt{X}_{\operatorname{\mathtt{Input}}}}$ in the inputs, the first
derivatives exhibit a first-order change
$\delta[\partial_{v}{\mathpzc{f}}]\triangleq\partial_{\operatorname{\mathtt{Input}},v}^{2}{\mathpzc{f}}\cdot\delta\mathtt{X}_{\operatorname{\mathtt{Input}}}$,
given by the Hessian-vector product. Applying chain-rule over $\mathcal{G}$
again for all terms in (12) we obtain,

$\displaystyle\forall v,\quad\delta[\partial_{v}{\mathpzc{f}}]=\sum_{s\in
v\cup\operatorname{\operatorname{ch}}(v)}\left(\sum_{a\in
v\cup\operatorname{\operatorname{pa}}(s)}\partial_{va}^{2}{\mathpzc{l}}_{s}\cdot\delta\mathtt{X}_{a}\right)+$
(13)
$\displaystyle\sum_{d\in\operatorname{\operatorname{ch}}(v)}\left(\delta[\partial_{d}{\mathpzc{f}}]^{\mathtt{T}}\partial_{v}d+\sum_{a\in\operatorname{\operatorname{pa}}(d)}(\partial_{d}{\mathpzc{f}}^{\mathtt{T}}\;\partial_{va}^{2}d)\cdot\delta\mathtt{X}_{a})\right);$
where, $\displaystyle\forall
a,\quad\delta\mathtt{X}_{a}=\sum_{d\in\operatorname{\operatorname{pa}}(a)}\partial_{d}a\cdot\delta\mathtt{X}_{{d}}.$

These equations can be solved, for a given
$\delta\mathtt{X}_{\operatorname{\mathtt{Input}}}$ by a forward-backward
recursion similar to the one used for solving (12). Computing the Hessian-
vector product in this manner takes time
$O(|E[\mathcal{\hat{G}}]|\omega(\mathcal{\hat{G}})^{2})$, where
$\omega(\mathcal{\hat{G}})$ is the clique number of the moralization of
$\mathcal{G}$.

###  A.2 Newton direction

Computing the Newton step requires inverting the Hessian-vector AD process:
find $\delta{\mathtt{X}_{\operatorname{\mathtt{Input}}}}$ such that,
$\delta[\partial_{\operatorname{\mathtt{Input}}}{\mathpzc{f}}]=-\partial_{\operatorname{\mathtt{Input}}}{\mathpzc{f}}$.
The inversion of these relations is related to the computation of the SQP
direction via 1,

See 1

_Proof._ The second equation in (13) is equivalent to
$\partial_{V}\Phi\cdot\delta S_{V}=-\Phi$, in (8). Rearranging the first
equation from (13), and setting
$\delta[\partial_{v}{\mathpzc{f}}]=-\partial_{v}{\mathpzc{f}}$ for all inputs,
we obtain $\forall v$,

$\displaystyle 0=$ $\displaystyle\sum_{\begin{subarray}{c}s\in
v\cup\operatorname{\operatorname{ch}}(v),\\\ a\in
v\cup\operatorname{\operatorname{pa}}(s)\end{subarray}}\partial_{va}^{2}{\mathpzc{l}}_{s}\;\delta\mathtt{X}_{a}+\sum_{\begin{subarray}{c}d\in\operatorname{\operatorname{ch}}(v),\\\
a\in\operatorname{\operatorname{pa}}(d)\end{subarray}}(\partial_{d}{\mathpzc{f}}^{\mathtt{T}}\;\partial_{va}^{2}d)\;\delta\mathtt{X}_{a}+$
(14)
$\displaystyle-\left(\begin{cases}\delta[\partial_{v}{\mathpzc{f}}]&\operatorname{\operatorname{pa}}(v)\neq\emptyset\\\
-\partial_{v}{\mathpzc{f}}&\mbox{otherwise}\end{cases}\right)+\sum_{d\in\operatorname{\operatorname{ch}}(v)}(\partial_{v}d)^{\mathtt{T}}\delta[\partial_{d}{\mathpzc{f}}].$

Similarly, expanding the top block in (8) using the definitions in (4) & (8),
we obtain $\forall v$,

$\displaystyle-\partial_{v}\mathcal{L}=$
$\displaystyle\sum_{\begin{subarray}{c}s\in
v\cup\operatorname{\operatorname{ch}}(v),\\\ a\in
v\cup\operatorname{\operatorname{pa}}(s)\end{subarray}}\partial_{va}^{2}{\mathpzc{l}}_{s}\;\delta\mathtt{X}_{a}+\sum_{\begin{subarray}{c}d\in\operatorname{\operatorname{ch}}(v),\\\
a\in\operatorname{\operatorname{pa}}(d)\end{subarray}}(\lambda_{d}^{\mathtt{T}}\;\partial_{av}^{2}d)\;\delta\mathtt{X}_{a}+$
(15)
$\displaystyle-\left(\begin{cases}\delta\lambda_{v}&\operatorname{\operatorname{pa}}(v)\neq\emptyset\\\
0&\mbox{otherwise}\end{cases}\right)+\sum_{d\in\operatorname{\operatorname{ch}}(v)}(\partial_{v}d)^{\mathtt{T}}\;\delta\lambda_{d},$

where,

$\displaystyle\partial_{v}\mathcal{L}=$ $\displaystyle\begin{cases}\sum_{s\in
v\cup\operatorname{\operatorname{ch}}(v)}\partial_{v}{\mathpzc{l}}_{s}+\sum_{d\in\operatorname{\operatorname{ch}}(v)}\lambda_{d}\cdot\partial_{v}d,&\operatorname{\operatorname{pa}}(v)=\emptyset\\\
\begin{array}[]{l}\sum_{s\in
v\cup\operatorname{\operatorname{ch}}(v)}\partial_{v}{\mathpzc{l}}_{s}+\sum_{d\in\operatorname{\operatorname{ch}}(v)}\lambda_{d}\cdot\partial_{v}d\\\
-\lambda_{v}\end{array},&\mbox{otherwise}\\\ \end{cases}$ (16)

The result follows from equations (12), (14), (15) & (eq. 16).

$\hfill\square$

##  Appendix B KKT complexity

In this section, we provide a Message Passing [MP] algorithm for solving KKT
systems arising in Algorithm 1 and show that it has a theoretical run-time
bound of $O(m\operatorname{tw}^{3})$, given a tree-decomposition.

###  B.1 Hypergraph structured QPs

For a hypergraph $\mathcal{H}$, denote the adjacency and incidence matrices by
$\mathcal{A}[\mathcal{H}]$ & $\mathcal{B}[\mathcal{H}]$ respectively,

$\displaystyle\mathcal{A}[\mathcal{H}]\in\mathbb{R}^{|V[\mathcal{H}]|\times|V[\mathcal{H}]|},\quad\quad\mathcal{B}[\mathcal{H}]\in\mathbb{R}^{|E[\mathcal{H}]|\times|V[\mathcal{H}]|},$
(17) $\displaystyle\mathcal{A}[\mathcal{H}]_{uv}=\begin{cases}1&\exists e\in
E[\mathcal{H}],\;u,v\in e\\\ 0&\mbox{otherwise}\end{cases}$
$\displaystyle\mathcal{B}[\mathcal{H}]_{eu}=\begin{cases}1&u\in e\\\
0&\mbox{otherwise}\end{cases}$

Given such a hypergraph $\mathcal{H}$, the family of QPs we’re interested in
solving is the following,

$\displaystyle\min_{x}\sum_{e\in E[\mathcal{H}]}{1\over
2}\mathtt{X}_{e}^{\mathtt{T}}Q_{e}\mathtt{X}_{e}-b_{e}^{\mathtt{T}}\mathtt{X}_{e},$
(18) $\displaystyle\forall e\in E[\mathcal{H}],\quad
G_{e}\mathtt{X}_{e}=h_{e}.$

Assuming that the QP has a bounded solution and that the constraints are full
rank, the minimizer to (eq. 18) is given by the solution to the following KKT
system,

$\displaystyle\left[\begin{array}[]{c c}Q&G^{\mathtt{T}}\\\
G&0\end{array}\right]\left[\begin{array}[]{c}x\\\ \lambda\end{array}\right]$
$\displaystyle=\left[\begin{array}[]{c}b\\\ h\end{array}\right],$ (19)
$\displaystyle x,b\in\mathbb{R}^{|V|},$
$\displaystyle\lambda,h\in\mathbb{R}^{M},$

where $Q,G,\lambda,x,b$ are concatenation of terms defined in (18)
respectively. The sparsity/support of (eq. 19) is closely related to
$\mathcal{H}$ because the quadratic part of the KKT equation has the sparsity
of the adjacency matrix, and the row of the constraint $G_{i,:}$ has the same
sparsity as some edge $e\in E[\mathcal{H}]$.

Tree decomposition: Extending the notion of Dynamic Programming to non-trees
(including Hypergraphs) requires a partitioning of the graph so as to satisfy
a _lifted_ notion of being a tree. Tree decomposition captures the essence of
such graph partitions,

###### Definition 1

(Tree decomposition) A tree-decomposition of a hypergraph $\mathcal{H}$
consists of a tree $\mathcal{T}$ and a map $\chi:V[\mathcal{T}]\rightarrow
2^{V[\mathcal{H}]}$, such that,

  1. i

(_Vertex cover_) $\cup_{i\in V[\mathcal{T}]}\chi(i)=V[\mathcal{H}].$

  2. ii

(_Edge cover_) $\forall e\in E[\mathcal{H}],\;\exists i\in
V[\mathcal{T}],e\subset\chi(i).$

  3. iii

(_Induced sub-tree_) $\forall u\in
V[\mathcal{H}],\;\mathcal{T}_{u}\triangleq\mathcal{T}[\\{i\in
V[\mathcal{T}]|u\in\chi(i)\\}]\;\mbox{is a non-empty subtree}$

The tree-width of a tree-decomposition $\mathcal{T}$ is defined to be
$\operatorname{tw}(\mathcal{T})=\max_{v\in V[\mathcal{T}]}|\chi(v)|-1$. The
tree-width of a graph $\mathcal{H}$ is defined to be the minimal tree-width
attained by any tree-decomposition of $\mathcal{H}$.

We define the vertex-induced subgraph in what follows to be
$\mathcal{H}[S]\triangleq(V[\mathcal{H}],\\{e\cap S,e\in E[\mathcal{H}]\\})$.
The following lemma ensures that such a decomposition ensures _local
dependence_.

###### Lemma 1

(Edge separation) Deleting the edge $xy\in E[\mathcal{T}]$, renders
$\mathcal{H}[V\backslash(\chi(x)\cap\chi(y))]$ disconnected.

Hypertree structured QP: The tree-decomposition itself can be considered a
Hypergraph, $(V[\mathcal{H}],\\{\chi(u),\forall u\in V[\mathcal{T}]\\})$. Such
a _Hypertree_ 222There are multiple definitions of a _Hypertree_ ; we use the
term to mean a maximal Hypergraph, whose tree-decomposition can be expressed
in terms of its edges. can also be thought of as a chordal graph. We assume
henceforth that the given graph $\mathcal{H}$ is a hypertree, and that
$\mathcal{T}$ is its tree-decomposition.

A message passing [MP] algorithm for solving (19) on such Hypertrees is given
below. The gather stage of the message passing algorithm is illustrated in
algorithm 2 333Note that the addition is performed vertex label-wise in Line 6
of Algorithm algorithm 2.. The function, Factorize, computes the partial LU
decomposition of its arguments; we describe below, its operation.

Denote the vertices that are interior to $l$ by $\iota=\chi(l)\cap\chi(p)$,
and those on the boundary (_i.e_ common to $p,l$) by
$\xi=\chi(l)\backslash\chi(p)$, and let
$r=\operatorname{rank}(\tilde{Q}_{\iota,\iota})$. The function computes
Gaussian-BP messages from block pivots $\\{2,3\\}$ to $\\{1,4\\}$ in (20).
Note that, unlike Gaussian-BP, the matrices in (20) are not necessarily
positive definite, but are however invertible.

Algorithm 2 Graphical QP

1: Given:
$\mathcal{T},\mathcal{H},\\{Q_{e}\\},\\{b_{e}\\},\\{G_{e}\\},\\{h_{e}\\}.$

2:

3: function GatherMessage($l,p,\mathcal{T}$)

4:
$(\tilde{Q}_{l},\tilde{b}_{l},\tilde{G}_{l},\tilde{h}_{l})\leftarrow(Q_{l},b_{l},G_{l},h_{l})$

5: for $c\in\operatorname{\delta}_{\mathcal{T}}(l)\backslash p$ do

6: $(Q_{c\rightarrow l},G_{c\rightarrow l},b_{c\rightarrow l},h_{c\rightarrow
l})\leftarrow\mbox{GatherMessage}(c,p,\mathcal{T})$

7:
$(\tilde{Q_{l}},\tilde{b_{l}})\leftarrow(\tilde{Q_{l}},\tilde{b_{l}})+(Q_{c\rightarrow
l},b_{c\rightarrow l})$

8: $\tilde{G}_{l}\leftarrow[\tilde{G}_{l};G_{c\rightarrow
l}],\tilde{h}_{l}\leftarrow[\tilde{h}_{l};h_{c\rightarrow l}]$

9: end for

10: return
Factorize($\chi(l),\chi(p),\tilde{Q}_{l},\tilde{b}_{l},\tilde{G}_{l},\tilde{h}_{l}$)

11:

12: function
Factorize($\chi(l),\chi(p),\tilde{Q},\tilde{b},\tilde{G},\tilde{h}$)

13: $(\xi,\iota)\leftarrow(\chi(l)\backslash\chi(p),\chi(l)\cap\chi(p))$

14: $r\leftarrow\operatorname{rank}(\tilde{Q}_{\iota,\iota})$

15: return Gaussian-BP messages from (20).

16: return

${\tilde{Q}_{\xi\xi}}$${\tilde{Q}_{\iota\xi}^{\mathtt{T}}}$${\tilde{G}_{:r,\xi}^{\mathtt{T}}}$${\tilde{G}_{r:,\iota}^{\mathtt{T}}}$${\tilde{Q}_{\iota\xi}}$${\tilde{Q}_{\iota\iota}}$${\tilde{G}_{:r,\iota}^{\mathtt{T}}}$${\tilde{G}_{r:,\iota}^{\mathtt{T}}}$${\tilde{G}_{:r,\xi}}$${\tilde{G}_{:r,\iota}}$${0}$${0}$${\tilde{G}_{r:,\xi}}$${\tilde{G}_{r:,\iota}}$${0}$${0}$$\left[\vbox{\hrule
height=39.28947pt,depth=39.28947pt,width=0.0pt}\right.$$\left.\vbox{\hrule
height=39.28947pt,depth=39.28947pt,width=0.0pt}\right]$${\mathtt{X}_{\xi}}$${\mathtt{X}_{\iota}}$${\lambda_{:r}}$${\lambda_{r:}}$$\left[\vbox{\hrule
height=33.96777pt,depth=33.96777pt,width=0.0pt}\right.$$\left.\vbox{\hrule
height=33.96777pt,depth=33.96777pt,width=0.0pt}\right]$$=$${\tilde{b}_{\xi}}$${\tilde{b}_{\iota}}$${\tilde{h}_{:r}}$${\tilde{h}_{r:}}$$\left[\vbox{\hrule
height=34.63446pt,depth=34.63446pt,width=0.0pt}\right.$$\left.\vbox{\hrule
height=34.63446pt,depth=34.63446pt,width=0.0pt}\right]$ (20)

Gaussian Belief-Propagation is essentially a re-statement of LU decomposition,
and consists of messages of the form,

$\displaystyle\mu_{i\rightarrow j}$ $\displaystyle:=[J_{i\rightarrow
j},h_{i\rightarrow j}]=[J_{ii},h_{i}]-\sum_{k\in\delta(i)\backslash
j}J_{ik}J_{k\rightarrow i}^{-1}[J_{ki},h_{k\rightarrow i}],$ (21)
$\displaystyle\mu_{i}$ $\displaystyle=J_{i\rightarrow j}^{-1}(h_{i\rightarrow
j}-J_{ij}\mu_{j}),$

where $J\mu=h$ is the equation that is to be solved. These can be replaced by
appropriate square-root forms to obtain instead, an LDL decomposition.

###### Theorem 3

The linear equation (19) can be solved in time
$O(|\mathcal{H}|\operatorname{tw}(\mathcal{H})^{3})$, given the tree-
decomposition, _via_ Algorithm 2.

_Proof._ The correctness of the algorithm follows from 1. The bound holds
trivially if,
$\operatorname{rank}{\tilde{G}}\leq\operatorname{rank}{\tilde{Q}_{\iota,\iota}}$,
at every step of the algorithm. Otherwise, by realizing that
$\tilde{G}_{l\rightarrow p}$, can’t have rank more than $|\chi(p)|$, the proof
follows.

$\hfill\square$