# Optimality of Finite-Support Parameter Shift Rules for Derivatives of
Variational Quantum Circuits

Dirk Oliver Theis
Theoretical Computer Science, University of Tartu, Estonia
<EMAIL_ADDRESS>
(Tue 22 Mar 2022 11:36:19 AM CET
Compiled: , )

###### Abstract

Variational (or, parameterized) quantum circuits are quantum circuits that
contain real-number parameters, that need to be optimized/“trained” in order
to achieve the desired quantum-computational effect. For that training,
analytic derivatives (as opposed to numerical derivation) are useful.
Parameter shift rules have received attention as a way to obtain analytic
derivatives, via statistical estimators.

In this paper, using Fourier Analysis, we characterize the set of all shift
rules that realize a given fixed partial derivative of a multi-parameter VQC.
Then, using Convex Optimization, we show how the search for the shift rule
with smallest standard deviation leads to a primal-dual pair of convex
optimization problems.

We study these optimization problems theoretically, prove a strong duality
theorem, and use it to establish optimal dual solutions for several families
of VQCs. This also proves optimality for some known shift rules and answers
the odd open question.

As a byproduct, we demonstrate how optimal shift rules can be found
efficiently computationally, and how the known optimal dual solutions help
with that.

Keywords: Parameterized quantum circuits, variational circuit optimization.

###### Contents

  1. 1 Introduction
  2. 2 Function spaces from expectation values of parameterized quantum circuits
    1. 2.1 A single Hamiltonian
    2. 2.2 Several Hamiltonians
    3. 2.3 Quantum circuits for arbitrary frequency sets
  3. 3 Finite-support shift rules, feasibility
  4. 4 System of equations and existence of (finite) feasible shift rules
  5. 5 Cost of a shift rule
    1. 5.1 Cost of a shift rule as an error norm
    2. 5.2 Cost of a shift rule as worst-case standard deviation
    3. 5.3 Cost of a shift rule as a heuristic for sparsity minimization
  6. 6 “Primal”, “dual” and weak duality
    1. 6.1 Weak duality
  7. 7 The computational perspective
    1. 7.1 Solving the dual optimization problem
    2. 7.2 Finding an optimal shift rule
    3. 7.3 Finding an optimal shift rule when a dual solution is known
  8. 8 Strong duality and complementary slackness
    1. 8.1 Strong duality
    2. 8.2 Complementary slackness
  9. 9 Existence of optimal finite-support shift rules, and size of the support
  10. 10 Optimal primal and/or dual solutions
    1. 10.1 Dimension $d=1$
    2. 10.2 Dimension $d\geq 2$
  11. 11 Conclusion
  12. A Example source code
  13. B Miscellaneous math background
    1. B.1 Lattices
    2. B.2 Convexity
      1. B.2.a Convex cones
  14. C Complexity of function estimation

##  1 Introduction

Near-term quantum algorithms meant to be run on pre-fault-tolerant quantum
computing devices often employ the concept of containing parameters, which are
being optimized/“trained” to maximize the utility of the quantum algorithm for
the particular application. The concept is known, depending on the application
area, as Variational Quantum Program, Parameterized Quantum Circuit, Quantum
Neural Network, etc (see [3] and the references therein).

Typical examples include the various types of Quantum Neural Networks for
quantum machine learning applications (e.g., [15, 2]), Variational Quantum
Eigensolvers for electronic structure computations (e.g., [23]), and, to
varying degrees, Quantum Approximate Optimization Algorithms [6].

One way of training the parameters of a variational or parameterized quantum
circuit deploys some type of smooth optimization algorithms such as
(stochastic) gradient descent or second order methods or whatnot, in other
words, it requires derivatives, be they first order (e.g., [15, 18]), or
higher order (e.g., [13, 19]). While, of course, the standard formulas of
numerical analysis yield approximations for derivatives, the interest is in
so-call analytic derivatives, which don’t approximate (with some
$\varepsilon>0$ that has to be chosen carefully) the derivatives, but estimate
the exact derivative. Tickling these analytic derivatives, of any order, out
of the quantum circuit is the subject matter of this paper.

Interestingly, the interest in efficient methods for estimating derivatives
seems to have become more and more interesting, e.g., [4, 9, 1, 11, 10, 22,
12].

In this paper, we consider parameterized quantum circuits which, with
parameters set to $x\in\mathbb{R}^{d}$, apply a quantum operation
$\mathcal{E}(x)$ to an initial quantum state $\varrho_{0}$, and in the
resulting state, a fixed observable $\mu$ is measured. The expectation value
$f(x):=\operatorname{tr}(\mu\mathcal{E}(x)(\varrho_{0}))$ is the quantity of
which one wishes to find the derivatives. In other words, we are interested in
the (partial) derivatives of the expectation value function

$f\colon\mathbb{R}^{d}\to\mathbb{R}\colon
x\mapsto\operatorname{tr}(\mu\mathcal{E}(x)(\varrho_{0})).$

Moreover, we consider the situation where each of the parameters enter the
parameterized quantum circuit via one-parameter subgroups of unitary
operators, i.e., mappings $U({\textrm{\tiny$\,\square\,$}})$ with $U(0)$ equal
to the identity and for all $s,t\in\mathbb{R}$, we have $U(s+t)=U(s)U(t)$.
Each parameter of the PQC could appear in one or several places of the quantum
circuit, with the same or different one-parameter subgroups. One-parameter
subgroups are, of course, of the form
$U({\textrm{\tiny$\,\square\,$}})=e^{2{\pi i}{\textrm{\tiny$\,\square\,$}}H}$
for a Hermitian operator $H$. (Note the author’s idiosyncratic scaling factor
$2\pi$, that matches the standard Fourier transform in Harmonic Analysis.)

This is the setting that seems to be considered in the majority of derivative
estimation research, although other settings exist, e.g., the “stochastic”
parameter shift rules of [1], which applies to parameterized unitaries of the
form $t\mapsto e^{i(tA+B)}$. Our results don’t cover these types of
parameterized unitaries.

To our knowledge, analytic derivatives of Parameterized Quantum Circuits have
been introduced into near-term quantum programming in [15].

Soon after [15], the trivial observation was made that, in full generality,
the expectation value functions of parameterized quantum circuits into which
the parameters enter through one-parameter subgroups can be characterized in
terms of their Fourier properties [20, 4, 7]. This fact can be exploited
directly for the optimization/training of the PQC [20, 17], but it also has
bearing on shift rules for derivatives, most notably in [22, 12], the main
references of the present paper.

Indeed, the present paper applies Harmonic Analysis and Convex Optimization to
prove the optimality, in terms of the standard deviation of the resulting
estimator, of some parameter shift rules. Unfortunately, the shift rules for
which we prove optimality are mostly known from [22]. However, our paper adds
the certainty that no better shift rules (in terms of stddev) exist.

##### Basic approach.

We find that things get simpler if we abstract from the quantum-circuit origin
of the functions and work with their Fourier analytic properties: Expectation
value functions are real valued, and for a parameterized quantum circuit with
$d$ parameters (it’s OK to count only parameters participating in the partial
derivative), there exists a finite set $\Xi\subset\mathbb{R}^{d}$ which
contains the Fourier spectrum of the expectation value function. As the
function is real valued, the set $\Xi$ is symmetric, meaning:

$\Xi=-\Xi\text{ where $-\Xi:=\\{-\xi\mid\xi\in\Xi\\}$.}$ (1)

We refer to finite, symmetric sets $\Xi\subset\mathbb{R}^{d}$ as frequency
sets of PQCs.

While every such set can occur as the support of the expectation value
function of a Parameterized Quantum Circuit (see Proposition 4 below) we
acknowledge that, for $d>1$, only a small subclass of frequency sets are
truely of interest, mostly because a (small) set $\Xi$ containing the Fourier
spectrum has to be somehow read off from the PQC before shift rules can be
constructed and applied. See §2 below for (mostly known) examples of frequency
sets of expectation value functions of PQCs which can be found by inspecting
the parameterized quantum circuit.

Before we summarize some of the contributions of this paper, we draw the
reader’s attention to Table 1 which lists the notations frequently used in
this paper, for easy reference.

Notation |  Description
---|---
$\mathbb{N}:=\\{1,2,3,\dots\\}$ |
$d\in\mathbb{N}$ |  Number of parameters of the function / parameterized quantum circuit.
$\alpha\in\mathbb{N}^{d}$ |  The order of the multi-variable partial derivative. The reason why we don’t allow $\alpha_{j}=0$ for some $j$’s is just convenience: We can restrict to the parameters that contribute in the partial derivative.
$\partial^{\alpha}=\prod_{j=1}^{d}\partial_{j}^{\alpha_{j}}$ |  Partial derivative in the math sense, see Wikipedia.
${\lvert{\alpha}\rvert}$ |  For $\alpha\in\mathbb{N}^{d}$: ${\lvert{\alpha}\rvert}:=\sum_{j=1}^{d}\alpha_{j}$.
$\Xi\subset\mathbb{R}^{d}$ |  A finite set (of frequencies), elements are typically $\xi\in\Xi$
$\Xi_{+}$ |  “Positive” part of $\Xi$: Arbitrary but fixed set satisfying the relation $\Xi=\\{0\\}\cup\Xi_{+}\cup(-\Xi_{+})$, where the union is disjoint. In case $d=1$, we think of $\Xi\\_\subset\left]0,\infty\right[$
$A\subset\mathbb{R}^{d}$ |  A (finite) set of points, elements are typically $a\in A$. The reader will understand the reason for writing $f(a)$ for (nice) functions $f\colon\mathbb{R}^{d}\to\mathbb{R}$. We appreciate that writing both $f(x)$ and $f(a)$, for the same function $f$, must be terribly confusing for physicists.
$\,\square\,$ |  Placeholder in expressions to simplify notations for functions. The placeholder is where you plug the variable, which as no name yet. E.g., as already used above, $e^{2{\pi i}{\textrm{\tiny$\,\square\,$}}H}$ stands for $v\mapsto e^{2{\pi i}vH}$. (There’s no preference for a particular letter; to emphasize that, with $v$ we deliberately made a weird choice.)
$u\bullet v$ |  $:=\sum_{j}u_{j}v_{j}$.
$\mathscr{K}_{\Xi}$ |  Real vector space of bounded real-valued functions whose Fourier spectrum is contained in the finite set $\Xi$; i.e., functions of the form $\sum_{\xi\in\Xi}b_{\xi}e^{2{\pi i}\xi\bullet x}$ with $b_{-\xi}=b_{\xi}^{*}$ for all $\xi\in\Xi$.
Table 1: Overview of notations. Habits of notation in Physics and Math /
Theoretical Computer Science are very different, which is why we list the
notations which are used a throughout this paper, particularly the arcane
ones.

##### The contributions of this paper.

(1) Using elementary Harmonic Analysis, for any $d\in\mathbb{N}$ and any given
finite symmetric set $\Xi\subset\mathbb{R}^{d}$, and for fixed
$\alpha\in\mathbb{N}^{d}$, we characterize the set of all shift rules which
realize the partial derivative111In the math sense,
en.wikipedia.org/wiki/Partial_derivative.
$\partial^{\alpha}:=\prod_{j=1}^{d}\partial_{j}^{\alpha_{j}}$ for every real-
valued function with Fourier spectrum contained in $\Xi$. We refer to these
shift rules as “feasible”. As we explain in §4, at the heart of our
characterization is a finite, underdetermined linear system of equations, not
unlike, but more general than those that occur in [22, 12] where the systems
are used to search for shift rules.

Our more general characterization is necessary because it allows us to (2)
realize the search for the optimal shift rule (among all possible shift rules)
as a the problem of minimizing a convex function over that set (§ 6).

The word “optimal”, of course, requires the concept of a “cost” of a shift
rule. We use a cost measure closely related to (the same, really) one proposed
in [22] and called the “number of shots” there. It can be understood as (2a)
an (asymptotic) standard deviation of an estimator in the worst case over all
parameterized quantum circuits; or, (2b) for the Functional/Numerical Analyst,
as the operator norm of the shift rule when viewed as a linear operator
between the obvious normed spaces. As usual, the latter perspective allows us
to view the “cost” as the worst-case bound on how noise in the input function
(e.g., sampled estimate of the expectation value function) is blown up in the
output function (e.g., estimate of the derivative of the expectation value
function).

As a nice coincidence, our cost measure can also be understood (2c) as a
heuristic to minimize the quantity that in [22] is called the “number of
unique circuits”, i.e., the number of “shift vectors” used — we refer to it as
the size of the support of the shift rule.

See §5 for the formal definition of our cost function, and for discussions
related to (2a,2b,2c).

We emphasize that in this paper the phrase “cost of a shift rule” refers to
one fixed mathematical concept (with three different justifications) — we do
not switch between different resource requirement concepts. However, we also
emphasize that much of the technical toolkit we develop carries through to
other cost functions: For other “linearizable” cost functions (our cost
function has that property, see §7.2) most of our tools still work; for
general convex cost functions the duality concept (see below) will suffer; for
cost functions that require on/off decisions (such as for the provably-optimal
minimization of the size of the support), we exit the land of comforting
convexity, and enter the realm of NP hardness.

As we are interested only in “finite” shift rules, i.e., where the size of the
support is finite, the resulting optimization problem has a finitness
condition which, ostensibly, falls out of the scope of standard convex
optimization. Nonetheless, (3) we establish a perfectly standard “dual” convex
optimization problem (§6.1), and prove two strong duality theorems: A pair of
primal and dual solutions are both optimal, if, and only if, (4a) so-called
complementary slackness holds, (4b) the objective function values are the
same; see §8.

(5) En passant in §7, we outline how, in the cases where optimal shift rules
cannot be found on paper, standard computational convex optimization
techniques when applied to the dual convex program yield optimal solutions to
it, and how, from these dual-optimal solutions, optimal finite feasible shift
rules can be obtained computationally. The algorithms are straightforward
consequences of our work on the duality theory (2,3,4).

######  Computational Method 1.

Let $d\in\mathbb{N}$ be a fixed parameter.

There exists an efficient algorithm for the following task:

> Given as input a finite symmetric set $\Xi\subset\mathbb{Q}^{d}$ and
> $\alpha\in\mathbb{N}^{d}$, and a precision parameter $\varepsilon>0$, return
> a near-optimal near-feasible shift rule (where “near” is controlled through
> the precision parameter).

Moreover, heuristically, in cases where there are several optimal shift rules,
one can hope for a shift rule with small support.

We also have a simpler method for the cases when optimal dual solutions can be
obtained analytically (which happens often, see below): There, a single Linear
Program has to be solved (§7.3).

The computational cost of the methods depends on $d$ and $1/\varepsilon$.

The precision parameter is always necessary, firstly simply because non-
rational numbers are involved due to the trigonometric functions that occur,
and secondly because the theoretical convergence results require them. While
everything else in this paper is (meant to be) fully rigorous, we will not
prove that our algorithm(s) run in polynomial time (for222The reason why $d$
has to be fixed is that the methods involve minimizing a function with Fourier
spectrum contained in $\Xi\subset\mathbb{Q}^{d}$, which can be done
efficiently only for fixed $d$.  fixed $d$): They are typical, standard
approaches in mathematical optimization that are widely used both in practice
(cutting-plane/column generation methods) and in computational complexity
theory (proofs of poly-time running times based on the Ellipsoid method).

Instead, we invite the reader to check out the supplementary source code
written in the mathematical computation programming language
Julia333julialang.org demonstrating the methods in the form of
Pluto444plutojl.org notebooks; see Appendix A. We point out that the source
code is only for demonstration of the algorithms, and will run into numerical
instabilities when the frequency sets get a little bit nasty.

In [22, 12], searches for shift rules based on more restricted (linear) system
characterizations are performed manually. But to our knowledge, there seems to
be no method to efficiently computationally find cost-wise (near-) optimal
shift rules for non-trivial frequency sets (even in dimension $d=1$), and they
don’t seem to be known analytically, either.

(6) Finally, we use the duality theorems to prove optimal shift rules for some
situations arising from parameterized quantum circuits, see §10.

For that, we need two distinguish two types of frequency sets — ones that
generate lattices in $\mathbb{R}^{d}$, which we call lattice generating, and
ones that generate dense subgroups of sub-vector spaces of $\mathbb{R}^{d}$.
We’ll elaborate on that distinction (known to every mathematician) below
(§10). For now we just mention that (a) a function is periodic if, and only
if, its Fourier spectrum is lattice generating (Appendix B.1); (b) a set is
lattice generating if, and only if, after a change of basis (“re-scaling”), it
is integral (mathematical folklore); (c) for example, for $d=1$, if every
element of the frequency set is a rational multiple of a fixed real constant,
then the frequency set is lattice generating.

(6a) For any $d,\alpha,\Xi$, we prove that a finite feasible optimal shift
rule exists, and Wierichs et al.’s quantity “number of unique circuits” (what
we call the size of the support of the shift rule) is at most
${\lvert{\Xi}\rvert}$ (§9). (Examples are known that show that the bound is
tight.)

(6b) For $d=1$, any frequency set $\Xi\subset\mathbb{R}$, and any555In this
paper, we give the complete proof only for $\alpha\leq 4$ and we refer to a
forthcoming paper for the proof of a technical lemma for $\alpha>4$.
$\alpha\in\mathbb{N}$, we quantify the cost of optimal shift rules: It is
$(2\pi\max\Xi)^{\alpha}$. In the simplest case, if the frequency set has the
form (for a positive integer $m$ and a positive real number $\xi_{1}$)

$\Xi=\\{k\xi_{1}\mid k\in\mathbb{Z},{\lvert{k}\rvert}\leq m\\},$ (2)

we prove the optimality of the shift rule that was given in [22] for that
case. But, for all other lattice generating $\Xi$, we prove that the
corresponding shift rule is still optimal in terms of its cost! In other
words, it is impossible for a shift rule to exploit “gaps” in the set of
frequencies to obtain lower cost. This closes an issue that was left open in
[22]. However, regarding the size of the support, in view of (6a), it is clear
that for sparse frequency sets, there are optimal shift rules with different
support sizes. As outlined above, our computational method can be
heuristically hoped to find a shift rule with small support.

Particularly for $d=1$, we show that an optimal shift rule exists whose
support size is at most ${\lvert{\Xi}\rvert}-(\alpha-1)\%2$. This means that
there are optimal shift rules among those described in general forms in in
[22, 12]. We point out that our simple formula for the cost of optimal finite-
support shift rules holds in full generality, i.e., also for non-lattice-
generating $\Xi$.

(6c) For $d\geq 2$ our results are not comprehensive. We quantify the cost of
optimal shift rules only for frequency sets satisfying a restrictive technical
condition (“pointy”). The condition includes the case when $\Xi$ is a product
set, which is the only one that if considered in [22] — it is questionable,
though, whether more complicated frequency sets can be identified a priori, by
inspecting the parameterized quantum circuit or performing a small number of
expectation-value estimations that probably shouldn’t depend on any other
parameters.

Note that in this paper, we don’t consider savings that can be gained by
computing several partial derivatives, e.g., a complete Hessian, as is done in
[22].

The results (6b) and (6c), described in detail in §10.1 and §10.2 resp., are
obtained by establishing a dual-optimal solution — we consider these dual-
optimal solutions to be a major contribution of this paper. Once we have dual-
optimal solutions, using our other theorems, the (existence and) cost of
optimal shift rules can easily be derived, and one of our computational
methods (see above) demonstrates a possible simple way to computationally find
a heuristically-sparse optimal shift rule even for rugged frequency sets,
based on the dual-optimal analytic solution.

##### About non lattice generating frequency sets.

It is prudent to point out a danger of misunderstanding [22]: There, the
innocent reader might think of a dichotomy between an equi-spaced Fourier
spectrum $\\{k\xi_{1}\mid k=-K,\dots,K\\}$ on the one hand, and non-
periodicity on the other. While functions with equi-spaced Fourier spectrum
are periodic, the converse is, of course, not true. Indeed, the periodicity of
the function is equivalent to the Fourier spectrum generating a lattice.

Wierichs et al. (Appendix B.1) talk about frequencies which cannot be scaled
to be integral, and use the phrase not periodic in that context, in the case
$d=1$. It is not clear to the author of the current paper what they mean.

Note also that a frequency set that is not lattice generating cannot be
represented in a computer: A representation of the vectors/points in any basis
will always require coefficients that are not rational.

Why bother with non-lattice-generating Fourier spectra at all? Clearly, from a
practical point of view, all frequency sets are lattice generating, and the
treatment of non-lattice-generating frequency sets brings some difficulty.
(For example, it is not clear whether an optimal feasible shift rule exists,
i.e., whether the infimum of the costs is not attained.)

The reason is that lattice generating frequency sets can sometimes blow up the
complexity. Consider, say, 10 31-bit unsigned fixed-point numbers
$\xi_{1},\dots,\xi_{10}\in\left]0,2\right[$ chosen uniformly at random, and
consider a corresponding (random) real-valued function that has, as its
Fourier spectrum, the set $\Xi:=\\{0\\}\cup\\{\pm\xi_{j}\mid j=1,\dots,10\\}$.
Clearly, $\Xi$ is lattice generating, but the probability that its period is
smaller than $2^{29}$ is at roughly one in a million. Hence, treating this
frequency set as lattice generating, while correct, does not make the problem
of finding a finite optimal shift rule feasible. Considering non-lattice-
generating frequency sets, on the other hand, and trading elusive optimality
(on account of non-representability in a computer memory) for near-optimality
(off by at most an arbitrary given $\varepsilon>0$), immediately leads to
tools such as simultaneous Diophantine approximation (and so to continued
fractions, LLL algorithm, etc).

While that offer sufficient motivation, we believe, for not excluding non-
lattice-generating frequency sets, in the present paper, we do not pursue the
algorithmic approaches that they inspire.

##### Outline of the paper.

In the next section, describe the function spaces we are working with. In
Section 3, we define (feasible) shit rules rigorously. In Section 4, we give
the characterization of the set of feasible shift rules. In Section 5, we
define the cost of a shift rule, and show the equivalence of 3 possible
definitions. In Section 6, we define the “primal” and “dual” optimization
problems for shift rules and prove the Weak Duality Theorem. In Section 7, we
interject a discussion on various aspects of the optimization problems that
can be computed in practice (and, indeed, in theory). In Section 8, we state
the Strong Duality and Complementary Slackness theorems for shift rules. In
Section 9, we discuss the existence of optimal finite feasible shift rules. In
Section 10, we give concrete results for sets $\Xi$ that have been discussed
in the literature. We wrap things up in the concluding Section 11.

##  2 Function spaces from expectation values of parameterized quantum
circuits

Fix a $d\in\mathbb{N}$, and let $\Xi$ be a fixed finite, symmetric subset of
$\mathbb{R}^{d}$. We define the vector space of functions

$\mathscr{K}_{\Xi}:=\bigl{\\{}f\colon\mathbb{R}^{d}\to\mathbb{R}\bigm{|}f\text{
smooth and bounded }\operatorname{Supp}\hat{f}\subseteq\Xi\\}$

The Fourier transform of the bounded smooth function here is, technically,
that of tempered distributions, but the finitness of the spectrum along with
the boundedness renders this machinery trivial: Indeed, we have

$\mathscr{K}_{\Xi}=\bigl{\\{}f\colon x\mapsto\sum_{\xi\in\Xi}b_{\xi}e^{2{\pi
i}\xi\bullet x}\bigm{|}b\in\mathbb{C}^{\Xi}\text{ with
}b_{-\xi}=b_{\xi}^{*}\bigr{\\}};$

with $u\bullet v:=\sum_{j}u_{j}v_{j}$.

> We emphasize our conventions for the Fourier transform: $\hat{f}(\xi)=\int
> e^{-2{\pi i}\xi\bullet x}f(x)\,dx$ — the crucial part being the $2\pi$ in
> the exponent! This choice goes hand in hand with our decision to normalize
> Hamiltonian evolution to $t\mapsto e^{2{\pi i}H\cdot t}$. It leads to
> factors of powers of $2\pi$ in the derivatives and hence the standard
> deviations of the estimators of the derivatives.

As indicated above, we find it adds clarity to the subject matter to think in
terms of function spaces rather than in terms of Hamiltonians or parameterized
quantum circuits. As this also adds some generality that is not required to
address the currently known applications, we find it prudent to restrict to
function spaces which arise from practically relevant parameterized quantum
circuits whenever that simplifies the mathematics. However, we refer to the
Proposition 4 below (§2.3) as a justification for working with function spaces
with arbitrary frequency sets: Every function can all be realized as the
expectation value function of a parameterized quantum circuit.

We now go through some common examples of parameterized quantum circuits and
their function spaces. What follows is well-known and treated in detail
elsewhere, e.g., [22].

###  2.1 A single Hamiltonian

The fact that evolution under a Hamiltonian produces a Fourier spectrum
consisting of the differences between the eigenvalues has been observed long
ago, e.g., [7, 20].

For easy reference, we pack the statement into a proposition. For convenience,
let us define the difference set

$\begin{split}\operatorname{Diff}\Lambda&:=\\{\lambda^{\prime}-\lambda\mid\lambda,\lambda^{\prime}\in\Lambda\\}\text{
for $\Lambda\subset\mathbb{R}^{d}$;}\end{split}$ (3)

and the spectrum $\operatorname{Spec}(A)$, the set of eigenvalues of an
operator $A$.

######  Proposition 2.

Consider a parameterized quantum circuit on $n$ qubits of the following form.
Let $\varrho$ be a $n$-qubit density matrix, and $\mu$ an $n$-qubit
observable. Let $H$ be a Hermitian operator that acts on an arbitrary subset
of the $n$ qubits. Consider the expectation value function

$f\colon\mathbb{R}\to\mathbb{R}\colon x\mapsto\operatorname{tr}(\mu\;e^{-2\pi
iH\cdot x}\,\varrho\,e^{2\pi iH\cdot x}).$ (4)

The Fourier spectrum of $f$ is contained in the set
$\operatorname{Diff}\operatorname{Spec}(H)$.

###### Proof.

Abbreviating $\Lambda:=\operatorname{Spec}(H)$ and replacing $H$ by its
spectral decomposition $H=\sum_{\lambda\in\Lambda}\lambda P_{\lambda}$, where
the $P_{*}$ are the spectral projectors, we find, for every $x\in\mathbb{R}$,

$\operatorname{tr}(\mu\;e^{-2\pi iH\cdot x}\,\varrho\,e^{2\pi iH\cdot
x})=\sum_{\lambda^{\prime},\lambda\in\Lambda}\operatorname{tr}(\mu e^{-2{\pi
i}\lambda x}P_{\lambda}\varrho e^{2{\pi
i}\lambda^{\prime}x}P_{\lambda^{\prime}})\\\
=\sum_{\lambda^{\prime},\lambda\in\Lambda}\operatorname{tr}(\mu
P_{\lambda}\varrho P_{\lambda^{\prime}})\cdot e^{2{\pi
i}(\lambda^{\prime}-\lambda)x}.$

∎

The function space associated with the case of Proposition 2 is
$\mathscr{K}_{\operatorname{Diff}\operatorname{Spec}(H)}$, the space of all
real-valued functions whose Fourier spectum is contained in
$\operatorname{Diff}\operatorname{Spec}(H)$.

###  2.2 Several Hamiltonians

In the multi-variate case, Hamiltonians $H_{1},\dots,H_{d}$ can be paired with
variables $x_{1},\dots,x_{d}$, respectively. For example, for arbitrary
quantum operations $\mathcal{E}_{1},\dots,\mathcal{E}_{d-1}$, and denoting by
$[U]$ ($U$ a unitary) the Kraus-form quantum operation $\varrho\mapsto
U\varrho U^{\dagger}$, we might have

$f\colon\mathbb{R}^{d}\to\mathbb{R}\colon\\\
x\mapsto\operatorname{tr}(\mu\,\cdot\,\left(\bigl{[}e^{-2\pi iH_{d}\cdot
x_{d}}\bigr{]}\circ\mathcal{E}_{d-1}\circ\bigl{[}e^{-2\pi iH_{d-1}\cdot
x_{d-1}}\bigr{]}\circ\dots\circ\mathcal{E}_{1}\circ\bigl{[}e^{-2\pi
iH_{1}\cdot x_{1}}\bigr{]}\right)(\varrho)).$ (5)

Applying Proposition 2 $d$ times yields, for the Fourier spectrum, the set

$\operatorname{Diff}\operatorname{Spec}(H_{1})\times\dots\times\operatorname{Diff}\operatorname{Spec}(H_{d}).$

We can add a little more color by applying a linear transformation, i.e.,
“$x=Au$” for an arbitrary real $(d\times e)$-matrix $A$.

######  Proposition 3.

If $f$ is as in (5), and $A$ is a real $(d\times e)$-matrix, then the function
$g\colon u\mapsto f(Au)$ has Fourier spectrum contained in the set

$\\{A^{\mspace{-1.0mu}\scriptscriptstyle\top\mspace{-1.0mu}}\xi\mid\xi\in\operatorname{Diff}\operatorname{Spec}(H_{1})\times\dots\times\operatorname{Diff}\operatorname{Spec}(H_{d})\\}.$

###### Sketch of proof..

Direct verification. With
$\Xi_{f}:=\operatorname{Diff}\operatorname{Spec}(H_{1})\times\dots\times\operatorname{Diff}\operatorname{Spec}(H_{d})$,
and
$\Xi_{g}:=\\{A^{\mspace{-1.0mu}\scriptscriptstyle\top\mspace{-1.0mu}}\xi\mid\xi\in\Xi\\}$,
there are complex numbers $b_{\xi}$, $\xi\in\Xi$, such that

$g(u)=f(Au)=\sum_{\xi\in\Xi_{f}}b_{\xi}\,e^{2{\pi
i}\xi\bullet(Au)}=\sum_{\xi\in\Xi_{f}}b_{\xi}\,e^{2{\pi
i}(A^{\mspace{-1.0mu}\scriptscriptstyle\top\mspace{-1.0mu}}\xi)\bullet
u}=\sum_{\eta\in\Xi_{g}}\biggl{(}\sum_{\begin{subarray}{c}\xi\in\Xi\\\
A^{\mspace{-1.0mu}\scriptscriptstyle\top\mspace{-1.0mu}}\xi=\eta\end{subarray}}b_{\xi}\biggr{)}\;e^{2{\pi
i}\eta\bullet u}.$

∎

The function space associated with the case of Proposition 3 is
$\mathscr{K}_{\Xi_{g}}$ (re-using the notation from the proof), the space of
all real-valued functions each of whose Fourier frequencies is of the form
$A^{\mspace{-1.0mu}\scriptscriptstyle\top\mspace{-1.0mu}}\xi$ for some
$\xi\in\operatorname{Diff}\operatorname{Spec}(H_{1})\times\dots\times\operatorname{Diff}\operatorname{Spec}(H_{d})$.

###  2.3 Quantum circuits for arbitrary frequency sets

We give the following quite obvious result as an argument in support of the
generality in which the present paper is written. No attempt has been made to
minimize the resources. The proof is in Appendix C.

######  Proposition 4.

Given a $d\in\mathbb{N}$, $\Xi\subset\mathbb{R}^{d}$ finite, symmetric, and a
$f\colon\mathbb{R}^{d}\to\mathbb{R}$ with Fourier spectrum equal to $\Xi$,
there exists a parameterized quantum circuit on at most
$d\log_{2}({\lvert{\Xi}\rvert}+1)$ qubits and with $d$ parameters, whose
expectation value function equals $f$.

Hence, we see that (a) the frequency sets of parameterized quantum circuits
can be completely arbitrary, and (b) every function in $\mathscr{K}_{\Xi}$
arises from a parameterized quantum circuit.

The practical value of Proposition 4 is doubtful, as taking derivatives
through shift rules requires that the frequency sets are known.

##  3 Finite-support shift rules, feasibility

In this paper, a shift rule is defined as a (finite complex Borel) measure
$\phi$ on $\mathbb{R}^{d}$, which is applied to a function
$f\in\mathscr{K}_{\Xi}$ via convolution:

$\phi*f(x)=\int f(x-u)\,d\phi(u)$

We are interested in finite shift rules, i.e., measures with finite support.
This restricts $\phi$ to choosing a finite set $A\subset\mathbb{R}^{d}$ and
complex numbers $u\in\mathbb{C}^{A}$ (we will later restrict, wlog, to real
coefficients $u\in\mathbb{R}^{A}$), and defining

$\phi_{A,u}:=\sum_{a\in A}u_{a}\delta_{a},$ (6)

where $\delta_{a}$ is the Dirac distribution / point measure at the point $a$,
which means $\int f(x)\,d\delta_{a}(x)=f(a)$ for all reasonable functions $f$.
For the whole convolution, we get

$\phi_{A,u}*f\,(x)=\sum_{a\in A}u_{a}\,\delta_{a}*f=\sum_{a\in A}u_{a}f(x-a).$

Given $\alpha\in\mathbb{N}^{d}$, a shift rule $\phi$ is called feasible for
$\Xi,\alpha$, if convolution against it realizes the $\partial^{\alpha}$
partial derivative for all real-valued functions with Fourier spectrum
contained in $\Xi$, i.e.,

$\phi*f=\partial^{\alpha}f\text{ for all $f\in\mathscr{K}_{\Xi}$.}$ (7)

From now on, when we say “shift rule”, we mean “finite (but not necessarily
feasible) shift rule”.

##  4 System of equations and existence of (finite) feasible shift rules

Given a finite set $A\subset\mathbb{R}^{d}$ and $u\in\mathbb{C}^{A}$, applying
the Fourier transform on both sides of equation (7), we find that $\phi$ is a
feasible shift rule if and only if

$\hat{\phi}(\xi)=(2{\pi i}\xi)^{\alpha}\qquad\text{for all $\xi\in\Xi$;}$ (8)

we refer to the supplementary material for the details. Using the fact that
$\widehat{\delta_{a}}(\xi)=e^{-2{\pi i}a\bullet\xi}$ for all
$\xi\in\mathbb{R}^{d}$, we can continue to an equivalent system of equations
formulated in the notation $\phi_{A,u}$ as in (6):

$\sum_{a\in A}u_{a}e^{-2{\pi i}a\bullet\xi}=(2\pi
i\xi)^{\alpha}\qquad\text{for all $\xi\in\Xi$.}$ (9)

As one can see, the system is linear in $u$, with coefficients that are
trigonometric functions of the vectors in the set $A$.

To summarize: $\phi_{A,u}$ as in (6) is a feasible shift rule for $\Xi,\alpha$
if, and only if, it is a solution to the the system of ${\lvert{\Xi}\rvert}$
equations (9).

It is not, up to this point, obvious that for every wild combination of
$\Xi$’s and $\alpha$’s a feasible (finite!) shift rule should even exist. In
the case where $\Xi$ is equi-spaced the existence is easy (and known in the
community, referred to as “equidistant frequencies” for $d=1$ in [22]). The
general case is discussed in [22, Eqn. (90)], but the author of the current
paper wasn’t able to find in Wierichs et al. any attempt of an argument why
their parameterized family of shift rules for the “general case” in $d=1$
should contain one that’s actually feasible (it’s probably obvious and I’m
just missing it). We prove the existence of finite feasible shift rules for
all $d,\Xi$ in the supplementary material.

We also note that, provided a single feasible shift rule exists, there exists
a whole “manifold”666Read: “lots!” of them.

##  5 Cost of a shift rule

There are several possibilities to define the “cost” of a shift rule. The
first one is what [22] calls “number of unique circuits” and what in our
terminology is the size (cardinality) of the support of the measure: For
$\phi_{A,u}$ as in (6), the support is (contained in) the set $A$, and its
size is ${\lvert{A}\rvert}$ (if all $u_{*}$ are non-zero).

In this paper the quantity we use as the cost of a shift rule $\phi$ is the
usual norm777The one which underlies the total-variation: For a complex
measure $\mu$ on $\Omega$, it is defined as
$\displaystyle{\lVert{\mu}\rVert}_{1}:=\sup_{f}\int f\,d\mu$, where the
supremum ranges over measurable functions with supremum norm at most $1$.  of
a finite measure, which we denote by
${\lVert{{\textrm{\tiny$\,\square\,$}}}\rVert}_{1}$. Specializing to our
situation, for $\phi_{A,u}$ as in (6), we have

$\operatorname{cost}(A,u):={\lVert{\phi_{A,u}}\rVert}_{1}=\sum_{a\in
A}{\lvert{u_{a}}\rvert}.$ (10)

As noted above, this definition is the same as Wierichs et al.’s “number of
shots” up to factors.

Before we discuss the justification(s) for the choice of this cost function,
we need to address complex vs. real measures. The observant reader has noticed
that we allow complex coefficients $u\in\mathbb{C}^{A}$ in the shift rules
(6), while in the literature, only real numbers have been considered. It can
be proven that for every feasible shift rule $\phi$, there always exists a
feasible shift rule that is a real signed measure, i.e., has only real
coefficients when written in the form (6)), and whose cost is at most that of
$\phi$.

Hence, for the rest of the paper, we consider real coefficients
$u\in\mathbb{R}^{A}$ only.

###  5.1 Cost of a shift rule as an error norm

The first justification of the choice of the cost function comes from
Harmonic/Numerical Analysis: It is a folklore fact888See [21] for the
background. The inequality “$\geq$” follows from ${\lvert{\int
f\,d\phi}\rvert}\leq\int{\lvert{f}\rvert}\,d{\lvert{\phi}\rvert}\leq{\lVert{f}\rVert}_{\infty}\int
d{\lvert{\phi}\rvert}={\lVert{f}\rVert}_{\infty}{\lVert{\phi}\rVert}_{1}$
where ${\lvert{\phi}\rvert}$ is the total variation measure; the inequality
“$\leq$” is in the Riesz Representation Theorem. that for every finite signed
(i.e., real valued) Borel measures $\phi$ on $\mathbb{R}^{d}$, we have

${\lVert{\phi}\rVert}_{1}=\sup\biggl{\\{}{\Bigl{\lvert}{\int
f\,d\phi}\Bigr{\rvert}}\biggm{|}f\in\mathscr{C}_{b}(\mathbb{R}^{d})\text{,
}{\lVert{f}\rVert}_{\infty}\leq 1\biggr{\\}},$ (11)

where we have denoted by $\mathscr{C}_{b}$ the Banach space of continuous
bounded functions, with the supremum norm
${\lVert{f}\rVert}_{\infty}:=\sup_{x}{\lvert{f(x)}\rvert}$.

This implies that

${\lVert{\phi}\rVert}_{1}=\sup_{{\lVert{f}\rVert}_{\infty}\leq
1}{\lVert{\phi*f}\rVert}_{\infty},$

where the supremum extends over the functions in
$\mathscr{C}_{b}(\mathbb{R}^{d})$ again. In other words,
${\lVert{\phi}\rVert}_{1}$ is equal to the operator norm for the operator

$D\colon\mathscr{C}_{b}(\mathbb{R}^{d})\to\mathscr{C}_{b}(\mathbb{R}^{d})\colon
f\mapsto\phi*f.$

As a general principle, norms of operators give error estimates. Suppose for
given $\Xi\subset\mathbb{R}^{d}$, $\alpha$, we are interested in the partial
derivative $\partial^{\alpha}$ of a function $f\in\mathscr{K}_{\Xi}$, and we
have a feasible finite-support shift rule $\phi$ that accomplishes this. But
$f$ is available only with an error (or noise), so the function that we can
evaluate is $\tilde{f}:=f+h$ where $h\in\mathscr{C}_{b}(\mathbb{R}^{d})$. In
general, the Fourier frequencies of the error function, $h$, can be outside of
the set of frequencies $\Xi$. In the case of parameterized quantum circuits,
where $\tilde{f}(x)$, for given $x$, is the average of discrete samples, $h$
can be expected to have arbitrarily large frequencies.

The question how much the error, $h$, messes up the computation of the
derivative via convolution with $\phi$ is upper bounded by the operator norm.
In the worst case (i.e., supremum norm) over all $x\in\mathbb{R}^{d}$, the
noise is “blown up” by a factor of ${\lVert{\phi}\rVert}_{1}$:

${\lVert{D\tilde{f}-Df}\rVert}_{\infty}={\lVert{Dh}\rVert}_{\infty}\leq{\lVert{D}\rVert}\cdot{\lVert{h}\rVert}_{\infty}={\lVert{\phi}\rVert}_{1}\cdot{\lVert{h}\rVert}_{\infty}.$

The inequality is tight, i.e., there exist $h$’s for which equality holds.

This reasoning is, of course, typical in any kind of Harmonic/Numerical
Analysis context.

###  5.2 Cost of a shift rule as worst-case standard deviation

In this paragraph, we essentially repeat with a little more detail the
calculation and approximation of Wierichs et al. that leads to their
definition of the “number of shots” in [22, Eqn. (14)].

Let us specialize to the case that the observable ($\mu$ in (4) and (5)) has
only $\pm 1$ eigenvalues, e.g., is a (tensor product of) Pauli obserable(s).
If $f$ is such an expectation value function, then, for all
$x\in\mathbb{R}^{d}$, the variance when evaluating $f(x)$ is $1-f^{2}(x)$.

As everybody knows, when the expectation of a finite sum $\sum_{a\in A}X_{a}$
of independent random variables $X_{a}$, $a\in A$, is to be estimated by
sampling each summand $X_{a}$ independently a number of $N_{a}$ times and
taking the average, then the variance-minimizing way to distribute a budget
$N$ of shots, $N=\sum_{a}N_{a}$, is to choose $N_{a}$ proportional to the
standard deviation of $X_{a}$, for all $a$. As this is only possible in
exceptional situations,999For each $a$, the standard deviation fraction
$N\cdot\sigma_{a}/\sum_{b}\sigma_{b}$ would have to be an integer.  this
perspective is asymptotic, i.e., we take $N\to\infty$. In this limit, the
resulting variance is asymptotic to $\frac{1}{N}(\sum_{a}\sigma_{a})^{2}$,
where $\sigma_{a}$ is the standard deviation of $X_{a}$, for all $a$.

In our case, when estimating a partial derivative at $x\in\mathbb{R}^{d}$ by
convolution with $\phi_{A,u}$ as in (6) (with $u\in\mathbb{R}^{A}$), we have
$X_{a}=u_{a}\cdot F(x-a)$, $a\in A$, where for all
$(F(y))_{y\in\mathbb{R}^{d}}$ is family of independent, $\pm 1$-valued random
variables with expectations $\operatorname*{\mathbb{E}}F(y)=f(y)$ for all $y$.
The variance is ${\lvert{u_{a}}\rvert}^{2}(1-f(x-a)^{2})$, and with optimally
distributed shots, the asymptotic variance is

$\tfrac{1}{N}\left(\sum_{a\in
A}{\lvert{u_{a}}\rvert}\sqrt{1-f^{2}(x-a)\,}\right)^{2}.$ (12)

We take the square root, to obtain the asymptotic standard deviation

$\tfrac{1}{\sqrt{N}}\sum_{a\in A}{\lvert{u_{a}}\rvert}\sqrt{1-f^{2}(x-a)\,}.$

To manage this expression, we are now taking the worst case over $x$ and $f$:
We are lower bounding $f^{2}(x-a)$ by $0$, for all $a\in A$. This results in
the upper bound of the standard deviation of, asymptotically,

$\tfrac{1}{\sqrt{N}}\sum_{a\in
A}{\lvert{u_{a}}\rvert}=\operatorname{cost}(A,u)/\sqrt{N}.$

Taking the worst case is not just a convenient primitive of computational
theory, in our case, it is also a bow to the dreaded “Barren Plateaus”
phenomenon [14], where the expectation value functions are flat. (Clearly,
lower bounding $f^{2}$ by a constant(!) other than $0$ doesn’t change the
comparison relations between asymptotic standard deviations of different shift
rules, and hence has no influence on which shift rule is optimal.)

We emphasize that we are very optimistic about the ingenuity of an estimator
applying a shift rule given its support: We allow it to somehow magically know
how to distribute a budget of shots among the support of the shift rule in an
optimal way, depending on the values of the expectation value function at the
support points! Alas, in the worst case, it doesn’t help at all.

###  5.3 Cost of a shift rule as a heuristic for sparsity minimization

In the literature, the other cost factor for a shift rule is what [22] call
“number of unique circuits”, and what we refer to as the size of the support
of the Borel measure, the set $A$ in the notation of (6). Let us assume for a
moment that our goal was not to minimize the standard deviation, but to
minimize the support, while a large but finite candidate set $A$ that allows
feasible shift rules with support contained in it has been fixed. The task at
hand is to find a solution to a (finite, underdetermined) system of linear
equations whose support (number of non-zero entries of a vector) has minimal
size:

$\begin{split}{}&\min\\{a\mid u_{a}\neq 0\\}\\\ \text{subject to
}&Eu=k\end{split}$

for an appropriately chosen matrix $E$ and vector $k$; see above. This is an
NP-hard optimization problem, but, as we all know, the standard heuristic for
it is to replace the term under the minimization by the 1-norm:
$\min\sum_{a}{\lvert{u_{a}}\rvert}$:

$\begin{split}{}&\min\sum_{a}{\lvert{u_{a}}\rvert}\\\ \text{subject to
}&Eu=k\end{split}$ (13)

In many cases, provable guarantees can be made on the effectiveness of that
heuristic [5].

Hence, minimizing ${\lVert{\phi}\rVert}_{1}$ over all feasible shift rules can
be justified completely independently from error norms and worst-case standard
deviations, simply as the goal to heuristically minimize the support of the
shift rule, in other words, maximize its sparsity.

##  6 “Primal”, “dual” and weak duality

In optimization vernacular, the “primal” optimization problem is the one you
want to solve, while the “dual” optimization problem is one which helps you in
solving or proving theorems about the primal. The two are related by “duality
theorems” which make statements about pairs of primal/dual feasible solutions.

From what we have discussed, given $d,\Xi,\alpha$, the optimal finite feasible
shift rule is the optimal solution (if attained) of the following optimization
problem ($\phi$ is the variable):

Infimum of $\displaystyle{\lVert{\phi}\rVert}_{1}$ (14a) subject to
$\displaystyle\phi*f=\partial^{\alpha}f\text{ for all
$f\in\mathscr{K}_{\Xi}$,}$ (14b) $\operatorname{Supp}\phi$ finite. (14c)

We have noted above that the ostensibly infinite system of linear equations
(14b) is equivalent to each of the (finite) systems (8) and (9), but for now,
we stick with the form in (14), because it gives an easy motivation for the
Weak Duality Theorem.

In terms of optimization theory, the finiteness condition is puzzling: If
$d,\Xi,\alpha$ are such that the infimum is not attained, do we have to pass
to some kind of infinite-dimensional closure of our set of finite-support
shift rules? We want to avoid that, mainly because it isn’t necessary, as we
will see later. The puzzlement is mitigated by considering the following
“dual” optimization problem ($f$ is the variable):

Supremum of $\displaystyle(-\partial)^{\alpha}f(0)$ (15a) subject to
$\displaystyle{\lVert{f}\rVert}_{\infty}\leq 1$ (15b) $\displaystyle
f\in\mathscr{K}_{\Xi}.$ (15c)

Note that
$(-\partial)^{\alpha}=(-1)^{{\lvert{\alpha}\rvert}}\partial^{\alpha}$. (The
notation ${\lvert{\alpha}\rvert}$, for $\alpha\in\mathbb{N}^{d}$, stands for
$\sum_{j}\alpha_{j}$. )

The optimization problem (15) is of a form that is quite common in convex
optimization, and both theory and computation are manageable. Indeed, the
feasible region (meaning, the set of all dual feasible $f$, i.e., the $f$’s
satisfying the constraints (15b),(15c)) is a subset of a finite dimensional
real vector space (namely $\mathscr{K}_{\Xi}$, of dimension
${\lvert{\Xi}\rvert}$), and the subset is defined by infinitely many linear
constraints (15b):

$\text{for each $a\in\mathbb{R}^{d}$:}\quad-1\leq f(a)\leq 1.$ (16)

Moreover, the objective function $f\mapsto(-\partial)^{\alpha}f(0)$ is linear.

We have the following critical (albeit super easy) fact.

######  Theorem 5.

Let $d\in\mathbb{N}$, $\Xi\subset\mathbb{R}^{d}$ finite symmetric, and
$\alpha\in\mathbb{N}^{d}$ be given.

The dual feasible region is non-empty and compact. In particular, the
objective function attains its supremum in at least one of its points.

For the proof of the theorem we refer to the supplementary material.

###  6.1 Weak duality

To give an idea why (15) can be considered a “dual” of (15), we will prove a
“weak duality theorem”:

> The objective function values of the feasible solutions to the maximization
> problem are less than or equal to the objective function values of the
> feasible solutions to the minimization problem.

We give a slightly more general statement where a restriction to finite
support of the measure is not necessary.

######  Proposition 6 (Weak Duality Theorem).

Let $\phi$ be a finite signed Borel measure on $\mathbb{R}^{d}$ satisfying
(14b), and let $f$ be feasible for (15). Then

$(-\partial)^{\alpha}f(0)\leq{\lVert{\phi}\rVert}_{1}.$ (17)

###### Proof.

Abbreviating $x\mapsto f(-x)$ to $f(-{\textrm{\tiny$\,\square\,$}})$, we
calculate:

$\displaystyle(-\partial)^{\alpha}f(0)$
$\displaystyle=\partial^{\alpha}(f(-{\textrm{\tiny$\,\square\,$}}))(0)$ (18)
$\displaystyle=\phi*(f(-{\textrm{\tiny$\,\square\,$}}))(0)$ [by (14b) and
(15c)] $\displaystyle=\int f(u-0)\,d\phi(u)$
$\displaystyle\leq{\lVert{f}\rVert}_{\infty}\cdot{\lVert{\phi}\rVert}_{1}$ [by
(11)] $\displaystyle\leq{\lVert{\phi}\rVert}_{1}.$ [by (15b)]

In the second equation, we have also used that for all $f$, if
$f\in\mathscr{K}_{\Xi}$ then
$f(-{\textrm{\tiny$\,\square\,$}})\in\mathscr{K}_{\Xi}$, as $\Xi$ is
symmetric. ∎

##  7 The computational perspective

In this paper, we mostly use optimization as a tool for proving theorems,
ultimately proving optimality of some feasible finite shift rules. The author
finds that demonstrating how the optimization problems (14) and (15) could be
solved computationally is not only fun, but also, importantly, helps
understand the theory.

This section is not mathematically rigorous. For example, we will omit cases
in the algorithms and gloss over details in the justifications. Making it
rigorous would require to discuss too much background Convex Optimization
theory. In other words, while we will try to make as concrete and practical as
possible how the methods work, we warn the reader that understanding why they
work is impossible without solid knowledge of Convex Optimization theory.

###  7.1 Solving the dual optimization problem

We start with the dual problem (15), because its solution is canonical. From a
computational theory perspective, the thrust of the matter is that, using the
Ellipsoid algorithm, you can solve the optimization problem (to provable
optimality) efficiently, if you can (could!) solve efficiently the so-called
“separation problem”, which is an algorithmic version of the separating
hyperplane theorem: Given a point, either assert that it is feasible (i.e.,
inside the convex set described by the constraints), or find an inequality
that is valid for all feasible points and violated by the given point.

We must supply a precision, $\varepsilon>0$, to the Ellipsoid algorithm
indicating that an “$\varepsilon$-almost feasible” solution whose objective
function value is at most $\varepsilon$ worse than the mathematical optimum is
sufficient; we refer to [8], Definition (2.1.10), for the exact formulation.
The running time will depend polynomially on $\log(1/\varepsilon)$. Similarly,
the separation problem should accept as input a precision $\delta>0$
(Definition (2.1.13) in [8]).

In this paper, where the computational aspect is illustrative and not
essential for the theoretical part, we ignore the computational complexity
theory entirely, describing the practitioner’s approach. In our situation,
from a practical point of view, the separation problem is the following:

> Given as input $d,\Xi$ and a function $f\in\mathscr{K}_{\Xi}$, return either
> the assertion ${\lVert{f}\rVert}_{\infty}\leq 1$, or an $a\in\mathbb{R}^{d}$
> such that ${\lvert{f(a)}\rvert}>1$ (where both inequalities are subject to a
> precision, $\delta>0$).

The algorithm for solving (15) then consists in iteratively performing the
following steps: Solve the optimization problem with the infinitely many
inequalities reduced to a finite list defined by a finite set
$A\subset\mathbb{R}^{d}$; feed the resulting optimal (but not necessarily dual
feasible) solution $f^{\star}$ into an algorithm solving the separation
problem; update the finite set $A$ of constraints and iterate, or exit with an
optimal dual (near-)feasible solution.

Let’s start by encoding functions in $\mathscr{K}_{\Xi}$. Arbitrarily fix
$\Xi_{+}\subset\Xi$ in such a way that $\Xi$ is the disjoint union of the
three sets $\\{0\\}$, $\Xi_{+}$, $-\Xi_{+}$. Now, we represent the
“$f\in\mathscr{K}_{\Xi}$” in (15) (and also in the input to the separation
problem) by its sine-cosine expansion:

$f=\mu+\sum_{\xi\in\Xi_{+}}\cos(2\pi{\textrm{\tiny$\,\square\,$}}\bullet\xi)\cdot
x_{\xi}-\sin(2\pi{\textrm{\tiny$\,\square\,$}}\bullet\xi)\cdot y_{\xi};$ (19)

the actual optimization variables are $\mu\in\mathbb{R}$,
$x,y\in\mathbb{R}^{\Xi_{+}}$. In that way the constraint (15c) is
automatically satisfied. The number of variables is, of course, equal to the
dimension of the real vector space $\mathscr{K}_{\Xi}$, namely
$1+2{\lvert{\Xi_{+}}\rvert}={\lvert{\Xi}\rvert}$.

The objective (15a) is obtained by taking $(-\partial)^{\alpha}$ at $0$ of the
function in (19). With the convenient abbreviations

$\displaystyle c_{\xi}$ $\displaystyle:=\Re((2{\pi
i}\xi)^{\alpha}),\quad\text{ and}$ (20) $\displaystyle s_{\xi}$
$\displaystyle:=\Im((2{\pi i}\xi)^{\alpha}).$

this gives101010A small calculation is necessary to see this. The reason why
we set things up like this is because it goes nicely with duality.

$\sum_{\xi\in\Xi_{+}}c_{\xi}\cdot x_{\xi}+s_{\xi}\cdot y_{\xi}.$ (21a)

The constraint (15b) is written out in infinitely many linear inequalities as
in (16), two for every $a\in\mathbb{R}^{d}$:

$-1\leq\quad\mu+\sum_{\xi\in\Xi_{+}}\cos(2\pi a\bullet\xi)\cdot
x_{\xi}-\sin(2\pi a\bullet\xi)\cdot y_{\xi}\quad\leq+1.$ (21b)

Fig. 1 shows the complete algorithm for solving (15) to optimality — provided
we could solve the separation problem. We have added a tolerance parameter
$\delta>0$ to replace the strict inequality in the separation problem: Instead
of the computationally unpractical ${\lvert{f(a)}\rvert}>1$ we require
${\lvert{f(a)}\rvert}>1+\delta$. As a consequence, the algorithm will, of
course, return only a near-feasible solution.

Input : $d,\Xi_{+}\subset\mathbb{R}^{d},\alpha\in\mathbb{N}^{d}$; $\delta\geq
0$

Output : Optimal near-feasible solution to (15) (in $\mu,x,y$ variables)

1 Initialize $A$ to some (finite) subset of $\mathbb{R}^{d}$; e.g.,
$A:=\emptyset$.

2 Repeat

3 Solve the following Linear Program, in the variables $\mu\in\mathbb{R}$,
$x,y\in\mathbb{R}^{\Xi_{+}}$: $\displaystyle\begin{aligned}
\qquad\qquad\textsl{maximize }&\textsl{expression~{}\eqref{LP:dual:obj}}\\\
\textsl{subject to }&\textsl{inequalities~{}\eqref{LP:dual:two-ieqs}, for all
$a\in A$}\end{aligned}$

4[1ex] Obtain an optimal111111We are glossing over a detail here: The LP may
be unbounded (indeed, if you start with $A=\emptyset$, that is the expected
outcome), and the LP solver will return the data of a half-line contained in
the LP feasible region. Generally, separation algorithms can be adapted to
handle that situation.  LP solution $\mu^{\star},x^{\star},y^{\star}$, with
objective function value $v^{\star}$.

5 Separation Routine  __

6 Initialize $A_{+}:=\emptyset$

7 Find argmin $a_{\text{min}}$ and argmax $a_{\text{max}}$ of the function

$\displaystyle\qquad\qquad f^{\star}\colon
a\mapsto\mu^{\star}+\sum_{\xi\in\Xi_{+}}\cos(2\pi a\bullet\xi)\cdot
x^{\star}_{\xi}-\sin(2\pi a\bullet\xi)\cdot y^{\star}_{\xi}$

8 If $f^{\star}(a_{\text{min}})<-1-\delta$, update
$A_{+}:=\\{a_{\text{min}}\\}$

9 If $f^{\star}(a_{\text{max}})>1+\delta$, update
$A_{+}:=A_{+}\cup\\{a_{\text{max}}\\}$

10 Return $A_{+}$

11Update $A:=A\cup A_{+}$

Until  _$A_{+}=\emptyset$_

12Return Optimal near-feasible solution $\mu^{\star},x^{\star},y^{\star}$ with
objective function value $v^{\star}$.

Algorithm 1 Solve-Dual Figure 1: Algorithm “Solve-Dual”

Iterated LP algorithms with “cutting plane generation” (i.e., adding
inequalities violated by the current optimal point), often based on the
Simplex algorithm for Linear Programming, are a standard tool in computational
optimization to solve to optimality convex optimization problems with linear
objective where there is a large number of inequalities.

We invite the reader to try it out: The author has coded a toy demonstration,
for $d=2$ only, of the algorithm as a notebook; see Appendix A.

###  7.2 Finding an optimal shift rule

A similar iterative approach, referred to as “column generation” in Linear
Programming literature, would help us with the primal problem (14). Since we
have the finiteness-of-the-support condition (14c), a minimum may or may not
be attained, but we can still hope for a shift rule that is feasible and whose
cost is at most $\varepsilon$ larger than the infimum, for a given
$\varepsilon>0$. In that approach, in every iteration, for a fixed finite set
$A\subset\mathbb{R}^{d}$, we would be solving the following Linear Program.
The variables are $u^{+},u^{-}\in\mathbb{R}^{A}$:

Minimum of $\displaystyle\sum_{a\in A}u^{+}_{a}+u^{-}_{a}$ (22a) subject to
$\displaystyle\begin{aligned} t]{4}&&&\sum_{a\in A}u^{+}_{a}-u^{-}_{a}&&=0,\\\
&\forall\;\xi\in\Xi_{+}\colon&&\sum_{a\in A}(u^{+}_{a}-u^{-}_{a})\cos(2\pi
a\bullet\xi)&&=c_{\xi},\\\ &\forall\;\xi\in\Xi_{+}\colon&&\sum_{a\in
A}(u^{+}_{a}-u^{-}_{a})\sin(-2\pi a\bullet\xi)&&=s_{\xi},\\\ \end{aligned}$
(22b) $\displaystyle u^{+},u^{-}\geq 0$ (22c)

The non-linear absolute value ${\lvert{u}\rvert}$ that occurs in (14) via (10)
has been “linearized” by replacing the optimization variable $u$ from (14) by
its positive and negative parts $u^{\pm}$, i.e., $u=u^{+}-u^{-}$.

Solving the LP will return either in the optimal feasible shift rule with
support contained in $A$, or the assertion that no feasible shift rule with
support contained in $A$ exists.

We are not pursuing the column-generation approach here, mainly because it
requires a lot of LP theory while not adding much magic. Indeed: The primal
column generation method is mathematically exactly the same the dual column
generation method: Each step in one is the LP-dual of the corresponding step
in the other.

###### Finding an optimal shift rule when a dual-optimal basis is available.

Instead, we now explain how we can computationally find a feasible finite-
support shift rule whose cost equals that returned by Algorithm 1, with
$\varepsilon,\delta$ set to $0$ to simplify the exposition. Indeed, the
support of the resulting shift rule will be contained in the last set $A$ that
was used to solve the LP in line 1 there. With the notations used in Algorithm
1, take the values of the variables in the iteration where the Separation
Routine returns $A_{+}=\emptyset$, i.e., the last iteration. We also assume
that the optimal solution is a so-called “vertex solution” which comes with
what is called a “basis” in LP vernacular: A set of inequalities that are
satisfied with equality by the vertex solution, and whose left-hand-side
vectors are linearly independent spanning a space of dimension
${\lvert{\Xi}\rvert}$. (In other words: The coefficients of the selected
inequalities form a square, invertible matrix.) We index the inequalities of
the basis by the (disjoint) sets $A^{+}$ and $A^{-}$, where the elements $a\in
A^{+}$ index inequalities of the form “$f(a)\leq 1$”, and the elements of
$a\in A^{-}$ index inequalities of the form “$-1\leq f(a)$”.

To obtain a feasible shift rule with support contained in $A^{+}\cup
A^{-}\subseteq A$ whose cost equals $v^{\star}$, we need to solve for
$u^{+},u^{-}$, the system of linear equations consisting of equations (22b),
and

$\displaystyle u^{+}_{a}$ $\displaystyle=0$ for all $a\notin A^{+}$ (23)
$\displaystyle u^{-}_{a}$ $\displaystyle=0$ for all $a\notin A^{-}$.

Linear Programming theory guarantees:

  * •

The existence of a solution;

  * •

That the obtained $(u^{+},u^{-})$ is a feasible solution for the LP (22)

  * •

That the cost of $u^{+},u^{-}$ in the LP (22) is equal to $v^{\star}$, and
hence optimal by weak LP duality.

Set $u_{a}:=u_{a}^{+}-u_{a}^{-}$ for all $a\in A$, et voilà, $\phi_{A,u}$
defined as in (6) is a the desired shift rule.

From the definition of an LP-basis above, we see that the size of the support
$A^{+}\cup A^{-}$ of the shift rule is at most
$2{\lvert{\Xi_{+}}\rvert}+1={\lvert{\Xi}\rvert}$.

In our source code, see Appendix A, we don’t have to explicitly solve the
system of linear equations: That is implicitly done inside the LP solver: We
can just query the result from the “dual LP-solution”.

###  7.3 Finding an optimal shift rule when a dual solution is known

We conclude this subsection by giving a method that, in the lattice-generating
case, given an optimal dual solution “on paper” but no dual optimal basis,
recovers an optimal shift rule by solving a single Linear Program.

Let $\Xi\subset\mathbb{R}^{d}$ be (finite, symmetric, and) lattice generating,
so that, wlog, we may assume that $\Xi\subset\mathbb{Z}^{d}$, so that the
functions in $\mathscr{K}_{\Xi}$ are $\mathbb{Z}^{d}$-periodic.

Denoting the known dual optimal solution by $f^{\star}$, let

$A^{\pm}:=\\{a\in\left[-{\nicefrac{{1}}{{2}}},{\nicefrac{{1}}{{2}}}\right[^{d}\mid
f^{\star}(a)=\pm 1\\},$

and (to simplify the notation) $A:=A^{+}\cup A^{-}$.

Now, the LP is easily described: it is (22), with the additional equations
(23).

The LP is of the form (13) from §5.3 in the motivation of the cost of a shift
rule based on a sparsity minimizing Linear Program. This means that, if, for a
given $A$, if several optimal shift rules exist, we can hope that one that is
sparse will be returned. That turns out to be true in practice: We have
implemented a demonstration of the method as a notebook, see Appendix A. The
sparsest optimal shift rule is not always found.

It is important to understand that, because of “complementary slackness” (see
below), every feasible solution of the LP is an optimal shift rule, and that
the objective function of the LP here does not contribute to optimality in
terms of cost of the shift rule, only, heuristically, to sparsity.

Now finally back to mathematical rigor. Phew!

##  8 Strong duality and complementary slackness

Weak duality already has useful consequences: If you happen to stumble upon a
feasible shift rule $\phi$ along with a dual-feasible function $f$, such that
their objective function values coincide, i.e.,
$(-\partial)^{\alpha}f(0)={\lVert{\phi}\rVert}_{1}$ holds, then you can
conclude that $\phi$ is an optimal shift rule, and $f$ is optimal for the
dual. If you are not lucky enough to have that equation, you can still
conclude from weak duality that the cost of your shift rule $\phi$ is at most
${\lVert{\phi}\rVert}_{1}-(-\partial)^{\alpha}f(0)$ larger than the infimum.

###  8.1 Strong duality

Strong duality adds to weak duality the certainty that if primal and dual
optimal solutions exist, then their objective function values coincide. Our
theorem is tailored for our situation.

######  Theorem 7 (Strong duality of shift rule optimization).

Let $d\in\mathbb{N}$, $\Xi$ finite symmetric $\subset\mathbb{R}^{d}$, and
$\alpha\in\mathbb{N}^{d}$ be given.

  1. (a)

The infimum in the primal optimization problem is equal to the maximum in the
dual optimization problem.

  2. (b)

If $\phi^{\star}$ is a finite feasible shift rule that is an optimal solution
of the primal optimization problem, then there exists a dual feasible $f$ such
that

${\lVert{\phi^{\star}}\rVert}_{1}=(-\partial)^{\alpha}f(0).$ (24)

The way how Theorem 7 is phrased might appear a bit odd once the complete
picture is known: Optimal finite-support shift rules always exist. However,
Theorem 7 follows from Linear Programming theory, while for the existence of
optimal shift rules, we have to raise the stakes to graduate level math.

###  8.2 Complementary slackness

“Complementary slackness” is the Linear Programming theory equivalent of the
Karush-Kuhn-Tucker (first-order) optimality conditions in general (convex)
optimization. KKT is the inequality-constrained extension of Lagrange
multipliers in multi-variable calculus. In the complementary slackness
conditions, the dual solution takes the role of the multiplier.

######  Theorem 8 (Complementary slackness for shift rule optimization).

Let $d\in\mathbb{N}$, $\Xi$ finite symmetric $\subset\mathbb{R}^{d}$, and
$\alpha\in\mathbb{N}^{d}$ be given.

If $\phi$ is a finite feasible shift rule, and $f$ dual feasible (i.e.,
feasible for (15)), then the following are equivalent:

  * •

Both $\phi$ and $f$ are optimal;

  * •

$\operatorname{Supp}\phi\subseteq\\{a\in\mathbb{R}^{d}\mid{\lvert{f(a)}\rvert}=1\\}$
in such a way that for all $a\in\mathbb{R}^{d}$ the following logical
implications hold:

$\displaystyle\phi(a)$ $\displaystyle<0$ $\displaystyle\quad\Longrightarrow$
$\displaystyle f(a)$ $\displaystyle=-1,$ (25a) and $\displaystyle\phi(a)$
$\displaystyle>0$ $\displaystyle\Longrightarrow$ $\displaystyle f(a)$
$\displaystyle=+1.$ (25b)

The proof of the theorem is in the supplementary material.

In applications of Linear Programming theory to other fields of theoretical
computer science / mathematics, complementary slackness is a
powerful121212Even if somewhat …vintage. tool for proving theorems. Not
surprisingly, Theorem 8 is our main ingredient when proving optimality for
some shift rules.

##  9 Existence of optimal finite-support shift rules, and size of the support

The attentive reader will have noticed that the question whether optimal
(finite!) shift rules exist, is unto this point unresolved. Here is the
answer.

######  Theorem 9.

Let $d\in\mathbb{N}$, let $\alpha\in\mathbb{N}^{d}$, and
$\Xi\subset\mathbb{R}^{d}$ non-empty and symmetric.

The optimization problem (14) attains an optimal solution, i.e., there exists
a feasible (finite) shift rule $\phi$ in which the infimum is attained.

Moreover, there always exists an optimal shift rule whose support has size at
most ${\lvert{\Xi}\rvert}$.

We also prove a strengthening of the size-of-the-support statement for $d=1$
to ${\lvert{\Xi}\rvert}-(\alpha\%2)$ (where “$\%$” stands for remainder upon
division). Wierichs et al. consider only shift-rules with that support, as
they come out of their particular construction of shift rules. But they never
ask the question whether shift rules with larger support but smaller cost
(i.e., “number of shots”) might exist.

######  Theorem 10.

Let $\alpha\in\mathbb{N}$, and $\Xi\subset\mathbb{R}$ non-empty and symmetric.
There always exists an optimal shift rule whose support has size at most
${\lvert{\Xi}\rvert}-(\alpha\%2)$.

##  10 Optimal primal and/or dual solutions

Before we approach the examples, we recall131313This is standard undergrad
math. the formal definition of what we call a “lattice generating” set of
vectors, as promised in the introduction.

Consider the sub-group (wrt addition) of $\mathbb{R}^{d}$ generated by $\Xi$.
There are two cases: This sub-group is a lattice, or this sub-group is dense
in a linear subspace of $\mathbb{R}^{d}$. A lattice is a subgroup that is
generated as a group by a linearly (over the real numbers) independent set of
vectors. For $\Xi$ to generate a lattice is equivalent to the matrix whose
columns are the elements of $\Xi$ having rank over the rational numbers equal
to (instead of strictly greater than) the rank over the real numbers. If $\Xi$
generates a lattice, we say, in this paper, that it is lattice generating.

###  10.1 Dimension $d=1$

######  Theorem 11.

Let $\alpha\in\mathbb{N}$, and $\Xi\subset\mathbb{R}$ non-empty and symmetric.
The optimal solution to the “dual” optimization problem (15) is

$a\mapsto(\partial^{\alpha}\cos)(2\pi\cdot\max\Xi\cdot a).$

In particular, the cost of an optimal shift rule is $(2\pi\max\Xi)^{\alpha}$.

Note that the theorem holds without conditions on $\Xi$ (beyond symmetry): It
holds whether $\Xi$ is equi-spaced or not; whether $\Xi$ is lattice-generating
or not; whether the space $\mathscr{K}_{\Xi}$ consists of periodic functions
or not. Indeed, we prove Theorem 11 in three steps, each of which is of
interest in its own right.

  1. 1.

For $\Xi=\\{-K,\dots,K\\}\subset\mathbb{Z}$ we establish primal and dual
feasible solutions and prove complementary slackness. Complementary slackness
relies on a fact about signs of derivatives of modified Dirichlet kernels; in
this paper we give the proof of that lemma only for the practically relevant
cases $\alpha=1,2,3,4$, and refer to a separate paper for the more technical
general case.

This proves that the shift rules in [22] for equi-spaced Fourier spectra are,
indeed, optimal.

  2. 2.

The case of lattice-generating Fourier spectra can be reduced to Step 1. by a
cheap trick.

As the proof establishes, in particular, the dual optimal value, and since an
optimal finite-support shift rule exists (Theorem 9), by strong duality
(Theorem 7(b)), we also know the cost of the best possible feasible shift rule
— which answers a question left open in [22].

Moreover, knowing the dual optimal solution allows us to realize the
computational search for a sparse optimal shift rule through standard Linear
Programming techniques: Solving a single LP as described in §7.3.

  3. 3.

Finally simultaneous Diophantine approximation can be used to reduce the non-
lattice-generating case to the lattice generating one.

###  10.2 Dimension $d\geq 2$

We give optimal dual solutions only in a special case, for which we need some
notation and terminology. Our notation for the canonical projection mappings
is the following:

$\operatorname{pr}_{j}\colon\mathbb{R}^{d}\to\mathbb{R}\colon x\mapsto x_{j}.$

Let us say that a frequency set $\Xi\subset\mathbb{R}^{d}$ is pointy, if

$\text{there exists $\epsilon\in\\{\pm 1\\}^{d}$ such that
}\Bigl{(}\epsilon_{1}\cdot\max\operatorname{pr}_{1}(\Xi),\epsilon_{2}\cdot\max\operatorname{pr}_{2}(\Xi),\dots,\epsilon_{d}\cdot\max\operatorname{pr}_{d}(\Xi)\Bigr{)}^{\mspace{-1.0mu}\scriptscriptstyle\top\mspace{-1.0mu}}.$
(26)

######  Theorem 12.

Let $d\geq 2$, $\alpha\in\mathbb{N}^{d}$, and $\Xi\subset\mathbb{R}^{d}$ non-
empty and symmetric.

If $\Xi$ is pointy with $\epsilon$ as in (26), then the optimal solution to
the “dual” optimization problem (15) is

$a\mapsto\prod_{j=1}^{d}(\partial^{\alpha_{j}}\cos)(2\pi\cdot\max\operatorname{pr}_{j}(\Xi)\cdot
a_{j}).$

In particular, the cost of an optimal shift rule is
$\prod_{j=1}^{d}(2\pi\max\operatorname{pr}_{j}(\Xi))^{\alpha_{j}}$.

We prove Theorem 12 in two steps.

  1. 1.

If $\Xi$ is a product set, i.e., if there exist symmetric sets
$\Xi_{j}\subset\mathbb{R}$, $j=1,\dots,d$, with $\Xi=\prod_{j=1}^{d}\Xi_{j}$,
then we can reduce to (11) together and Theorem 9.

  2. 2.

If $\Xi$ is pointy, then it is contained in a fitting product set, so the same
cheap trick we used above gives a reduction to 1..

##  11 Conclusion

This paper contributes what the author believes to be novel and useful tools
to the study of efficient shift rules, culminating in dual optimal solutions
for a number of situations. The tools inspire efficient computational methods
for practically computing optimal, hopefully sparse, shift rules for the many
situations where (at least to date) no analytic description of the optimal
(sparse) shift rules are known.

Some questions remain, some of them “merely” mathematical, others with
practical relevance.

##### Dual optimal solutions for $d>2$.

The present papers offers dual optimal solutions for $d>2$ only for a
restricted family of frequency sets (the “pointy” ones). Is it possible to
write down and prove dual optimal solutions for all (lattice generating)
frequency sets? Is it true that (in the lattice generating case), the dual
optimal solution is always of the form
$a\mapsto\prod_{j=1}^{d}(\partial^{\alpha_{j}}\cos)(2\pi\xi_{j}\bullet a_{j})$
for a $\xi$ in the lattice generated by $\Xi$? Or even $\xi\in\Xi$?

##### Numerically stable implementations of the methods computing shift rules.

In the supplementary material (see Appendix A) we have given quick-and-dirty
implementations of the computational methods in this paper, based on cheap
Linear Programming solvers, and with no thought given to stability. Obviously,
a proper, numerically stable implementation of these methods is needed in
order to apply them practically. The most interesting method is certainly that
in §7.3, when the optimal dual solution is known on paper as it allows to
exploit all the theory of this paper.

##### Several partial derivatives at the same point

Wierichs et al. give combined shift rules for estimating the whole Hessian. Is
their construction optimal? The present paper only shows how to estimate each
individual entry of the Hessian optimally, but it does not consider estimators
which “recycle” previously made measurements to estimate several partial
derivatives at the same time. (It can be shown that, for estimating $\partial
f(x)$, measurements at $x$ cannot be used to reduce the cost, but that’s a
silly scenario.)

## Acknowledgments

This research was partly funded by the Estonian Centre of Excellence in
Computer Science (EXCITE), and by Estonian Research Agency (Eesti
Teadusagentuur, ETIS) through grant PRG946.

## APPENDIX

##  Appendix A Example source code

The author of this paper has implemented some of the techniques discussed here
in the mathematical computation programming language Julia141414julialang.org,
in the form of Pluto151515plutojl.org notebooks.

The notebooks are available on GitHub161616github.com/dojt/ShiftRulesPluto for
download. Here we go through the files in that repository with a quick
description of what the notebook does.

##### File solve-for-u.jl

The main function in this notebook, `solve_it(;`$\Xi_{+}$`, `$\alpha$`,
`$A$`)`, takes a set $\Xi_{+}\subset\mathbb{Q}$ of (positive) frequencies, and
a set $A\subset\mathbb{R}$, and finds a shift rule of the form (6). If a
feasible shift rule with that support exists, it will be found; if none
exists, one that is closest to being feasible is returned. Closest to being
feasible refers to the 2-norm, in frequency space
${\lVert{\partial^{\alpha}-\phi}\rVert}_{2,\Xi}$.

The function uses floating point arithmetic, and the precision can be chosen
between 64 bits and 16384 bits. Changing the precision can sometimes help
recognizing a set $A$ that yields no feasible shift rule, but where the error
with 64-bit precision is small, say $10^{-16}$.

##### File dual-opt.jl

This notebook implements the method that is described in the subsection §7.1.
After terminating with an optimal set $A$, the optimal shift rule is obtained:
This can be done directly from the dual variables, which, is equivalent to the
method described in §7.2, in the case when an optimal basis is known.

##### File sparse-optimal-1d.jl

In this notebook, knowledge about the optimal dual solution for lattice-
generating frequency sets with $d=1$ is used to compute an optimal feasible
shift rule using a finite, i.e., “normal”, Linear Program, as described in
§7.3.

##  Appendix B Miscellaneous math background

###  B.1 Lattices

We sketch the proof of the equivalence of the periodicity of a function and
the lattice-generating property of its Fourier spectrum:

######  Lemma 13.

If $f$ is a tempered distribution on $\mathbb{R}^{d}$ that is periodic (i.e.,
there exists an additive subgroup $\mathrm{M}$ of $\mathbb{R}^{d}$ such that
$\delta_{\mu}*f=f$ for all $\mu\in\mathrm{M}$), then the support of the
Fourier transform of $f$ is a lattice.

If a function $f\colon\mathbb{R}^{d}\to\mathbb{R}$ has a Fourier spectrum
$\Xi$ that generates a lattice $\Lambda$, then the function is periodic modulo
the dual lattice (which is, indeed, a lattice)

$\Lambda^{\bot}:=\\{\mu\in\mathbb{R}^{d}\mid\forall\lambda\in\Lambda\colon\mu\bullet\lambda=0\\},$

i.e., for all $x\in\mathbb{R}^{d}$ and every $\mu\in\Lambda^{\bot}$ we have
$f(x+\mu)=f(x)$.

For the other direction, we skip the part where you have to prove that the
Fourier spectrum is discrete, and go right to the part where we prove that it
also generates a discrete set (i.e., a lattice). If a smooth, bounded function
$f\colon\mathbb{R}^{d}\to\mathbb{R}$ is periodic, meaning, there exists an
additive subgroup $\mathrm{M}$ of $\mathbb{R}^{d}$ such that for all
$x\in\mathbb{R}^{d}$ and all $\mu\in\mathrm{M}$ we have $f(x+\mu)=f(x)$, then
either $f=0$ or $\mathrm{M}$ is a lattice. In the case when $\mathrm{M}$ is a
lattice, $\hat{f}=e^{2{\pi i}\mu\bullet{\textrm{\tiny$\,\square\,$}}}\hat{f}$,
which implies that, for all $\xi$ in the support of $\hat{f}$, we have
$\mu\bullet\xi\in\mathbb{Z}$, in other words $\Xi\subset\mathrm{M}^{\bot}$, so
that the group generated by $\Xi$ is also contained in $M^{\bot}$, and hence a
lattice.

###  B.2 Convexity

#### B.2.a. Convex cones

For a set of vectors $R$, let us denote by $\operatorname{cvxcone}(R)$ the
convex cone generated by $R$, i.e., the set of all (finite) linear
combinations with non-negative coefficients of elements of $R$. A convex cone
$C$ is called polyhedral, if there exists a finite set $R$ such that
$C=\operatorname{cvxcone}R$.

###### Lemma on infinitely-generated convex cones.

We use the following math notation only in this subsection. For a set $C$ in a
finite-dimensional real vector space, we denote: by $\bar{C}$ its topological
closure; and by $C^{\circ}$ its topological interior. Recall that, if $C$ is a
convex set (in a finite-dimensional real vector space) then $\bar{C}$ is also
convex.

The following strengthened version of Carathéodordy’s theorem is one of those
things in convexity theory about which the professor says “It’s obvious” and
then a student asks, “Why?”, and the professor spends 60min proving it.
Carathéodordy’s theorem is the case $\bar{C}=C$.

######  Lemma 14.

Let $C$ be a convex cone in an $m$-dimensional real vector space, and suppose
$C$ is generated by a set of vectors $R$. If $c\in(\bar{C})^{\circ}$ then
there exist $r_{1},\dots,r_{m}\in R$ with
$c\in\operatorname{cvxcone}\\{r_{1},\dots,r_{m}\\}$.

##  Appendix C Complexity of function estimation

In this section we briefly consider the following questions, ignoring
“subtleties” about processing real numbers in classical algorithms.

> Given a $d\in\mathbb{N}$, $\Xi\subset\mathbb{R}^{d}$ finite & symmetric, and
> a $f\colon\mathbb{R}^{d}\to\mathbb{R}$ with Fourier spectrum equal to $\Xi$:
>
>   1. (1)
>
> Does there exist a parameterized quantum circuit (with $d$ parameters) whose
> expectation value function is $f$?
>
>   2. (2)
>
> Does there exist a uniform family of these circuits, i.e., a polynomial time
> classical algorithm which, upon input of $d,\Xi,\hat{f}$, produces the PQC
> as in (1)?
>
>   3. (3)
>
> How does the minmum depth (or other complexity measure) of the quantum
> circuits in (1) relate to the complexity of estimating the function with a
> (randomized) classical circuit/algorithm?
>
>

The questions are informally inspired by the complexity theory of Boolean
functions. We give the obvious, cheap answer to (1), Proposition 4, as an
argument in support of the generality in which the present paper has been
written, making absolutely no attempt to search for an efficient parameterized
quantum circuit — and express the hope that the other questions might be
picked up by smart people.

Recall the definition of $\operatorname{Diff}$ from (3), and, for conciseness,
denote by $\operatorname{Spec}(A)$ the set of eigenvalues of a given operator
$A$.

###### Proof of Proposition 4.

For each $j=1,\dots,d$, consider the set $\Xi_{j}=\\{\xi_{j}\mid\xi\in\Xi\\}$
— the projection of $\Xi$ onto the $j$th coordinate. For each $j=1,\dots,d$,
apply Lemma 15 below to each of the $\Xi_{j}$, and define
$q_{j}:=\lceil{\log_{2}(({\lvert{\Xi_{j}}\rvert}+1)/2)}\rceil$. Letting

$q:=\sum_{j=1}^{d}q_{j}\quad<d\,\log_{2}({\lvert{\Xi}\rvert}+1),$

we take $q$ qubits, split into $d$ quantum registers, where, for
$j=1,\dots,d$, register number $j$ has $q_{j}$ qubits. Let us denote by
$H_{j}$ the operator from Lemma 15 below applied to $\Xi_{j}$, acting on the
$j$th quantum register. We then let

$H:=\sum_{j=1}^{d}H_{j},$

and consider the parameterized circuit with initial state
${\lvert\psi_{0}\rangle}$, multi-parameterized unitary evolution $U(x)$, as
follows:

$\begin{split}{\lvert\psi_{0}\rangle}&:=2^{-q/2}\sum_{\gamma\in\\{0,1\\}^{q}}{\lvert\gamma\rangle}\\\
U(x)&:=e^{2{\pi i}\sum_{j=1}^{d}x_{j}H_{j}}.\end{split}$

Now we define a Hermitian $q$-qubit operator $\mu$ as follows.

First of all, we interpret each of the $d$ quantum registers in the
computational basis as encoding an unsigned integer. For
$\ell_{j},\ell^{\prime}_{j}\in\\{0,\dots,2^{q_{j}}\\}$, $j=1,\dots,d$, we say
that $(\ell,\ell^{\prime})$ fits, if, for all $j=1,\dots,d$, we have both
$\ell^{\prime}_{j},\ell_{j}\leq({\lvert{\Xi_{j}}\rvert}-1)/2$ and
$\ell_{j}\cdot\ell^{\prime}_{j}=0$. Now, for
$\ell_{j},\ell^{\prime}_{j}\in\\{0,\dots,2^{q_{j}}\\}$, we define

${\langle\ell^{\prime}_{1},\ell^{\prime}_{2},\dots,\ell^{\prime}_{d}\rvert}\mu{\lvert\ell_{1},\ell_{2},\dots,\ell_{d}\rangle}:=\begin{cases}2^{q/2}\hat{f}(\xi),&\text{if
$(\ell,\ell^{\prime})$ fits and
}{}\Bigl{(}\xi_{\ell^{\prime}_{1}}-\xi_{\ell_{1}}\,,\;\dots\,,\;\xi_{\ell^{\prime}_{d}}-\xi_{\ell_{d}}\Bigr{)}^{\mspace{-1.0mu}\scriptscriptstyle\top\mspace{-1.0mu}}\in\Xi\\\
0,&\text{otherwise.}\end{cases}$

Note that, as $f$ is real valued, $\mu$ is indeed Hermitian.

The final quantum circuit is the following: Starting from ${\lvert 0\rangle}$,
prepare a uniform superposition over all computational basis states, then
apply the parameterized unitary (using the usual tools, as the exponent is a
linear combination of tensor-products of $Z$-operators), finally measure
$\mu$. The last step is done by diagonalizing $\mu$ using a basis change
unitary $W$ from the computational basis, i.e., $W\mu W^{\dagger}$ is diagonal
in the computational basis, and implementing $W$ as a quantum circuit, based
on the usual universality theorems (e.g., Section 4.5.2 in [16]). ∎

######  Lemma 15 (Trivial Helper Lemma).

Let $\Xi\in\mathbb{R}$ be finite and symmetric, and let
$0=:\xi^{(0)}<\xi^{(1)}<\dots<\xi^{(m)}$, where
$m:=({\lvert{\Xi}\rvert}-1)/2$, be the non-negative elements of $\Xi$. There
exists a Hermitian $\lceil{\log_{2}(m+1)}\rceil$-qubit operator $H$ such that
$\operatorname{Diff}\operatorname{Spec}(H)\supseteq\Xi$.

More precisely, the operator $H$ is diagonal in the multi-qubit Z-basis
(computational basis), and for $\ell,\ell^{\prime}\in\\{0,\dots,m\\}$ (coded
in binary):

${\langle\ell^{\prime}\rvert}H{\lvert\ell\rangle}=\delta_{\ell^{\prime},\ell}\cdot\xi^{(\ell)}.$
(27)

Let $\gamma\in\mathbb{N}^{q}$, and a 1-qubit operator $A$, we use the
shorthand

$A^{\otimes\gamma}:=A^{\gamma_{1}}\otimes A^{\gamma_{2}}\otimes\dots\otimes
A^{\gamma_{q-1}}\otimes A^{\gamma_{q}}.$

###### Proof of Lemma 15.

To define $H$, follow (27), and choose the remaining $2^{q}-m-1$ diagonal
entries arbitrarily. We have
$\operatorname{Diff}\operatorname{Spec}(H)\supseteq\operatorname{Diff}\\{\xi^{(0)},\dots,\xi^{(m)}\\}\supseteq\Xi$.
∎

## References

  * [1] Leonardo Banchi and Gavin E. Crooks.  Measuring analytic gradients of general quantum evolution with the stochastic parameter shift rule.  Quantum, 5:386, January 2021.
  * [2] Marcello Benedetti, Erika Lloyd, Stefan Sack, and Mattia Fiorentini.  Parameterized quantum circuits as machine learning models.  Quantum Science and Technology, 4(4):043001, 2019.
  * [3] Marco Cerezo, Andrew Arrasmith, Ryan Babbush, Simon C Benjamin, Suguru Endo, Keisuke Fujii, Jarrod R McClean, Kosuke Mitarai, Xiao Yuan, Lukasz Cincio, et al.  Variational quantum algorithms.  Nature Reviews Physics, 3(9):625–644, 2021.
  * [4] Gavin E. Crooks.  Gradients of parameterized quantum gates using the parameter-shift rule and gate decomposition.  Preprint, 2019.  arXiv:1905.13311.
  * [5] David L Donoho and Jared Tanner.  Sparse nonnegative solution of underdetermined linear equations by linear programming.  Proceedings of the National Academy of Sciences, 102(27):9446–9451, 2005.
  * [6] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A quantum approximate optimization algorithm.  arXiv preprint arXiv:1411.4028, 2014.
  * [7] Javier Gil Vidal and Dirk Oliver Theis.  Input redundancy for parameterized quantum circuits.  Preprint, 2019.  arXiv:1901.11434.
  * [8] M. Grötschel, L. Lovász, and A. Schrijver.  Geometric algorithms and combinatorial optimization, volume 2 of Algorithms and Combinatorics.  Springer, 2nd edition, 1993.
  * [9] Thomas Hubregtsen, Frederik Wilde, Shozab Qasim, and Jens Eisert.  Single-component gradient rules for variational quantum algorithms.  Preprint, 2021.  arXiv:2106.01388v1.
  * [10] Artur F. Izmaylov, Robert A. Lang, and Tzu-Ching Yen.  Analytic gradients in variational quantum algorithms: Algebraic extensions of the parameter-shift rule to general unitary transformations.  Preprint, 2021.  arXiv:2107.08131.
  * [11] Jakob S. Kottmann, Abhinav Anand, and Alán Aspuru-Guzik.  A feasible approach for automatically differentiable unitary coupled-cluster on quantum computers.  Chemical Science, 12(10):3497–3508, 2021.
  * [12] Oleksandr Kyriienko and Vincent E. Elfving.  Generalized quantum circuit differentiation rules.  Preprint, 2021.  arXiv:arXiv:2108.01218.
  * [13] Andrea Mari, Thomas R. Bromley, and Nathan Killoran.  Estimating the gradient and higher-order derivatives on quantum hardware.  Phys. Rev. A, 103:012405, Jan 2021.
  * [14] Jarrod R McClean, Sergio Boixo, Vadim N Smelyanskiy, Ryan Babbush, and Hartmut Neven.  Barren plateaus in quantum neural network training landscapes.  Nature communications, 9(1):4812, 2018.
  * [15] Kosuke Mitarai, Makoto Negoro, Masahiro Kitagawa, and Keisuke Fujii.  Quantum circuit learning.  Phys. Rev. A, 98:032309, 2018.  arXiv:1803.00745.
  * [16] Michael A. Nielsen and Isaac L. Chuang.  Quantum Computation and Quantum Information.  Cambridge University Press, 2000.
  * [17] Mateusz Ostaszewski, Edward Grant, and Marcello Benedetti.  Structure optimization for parameterized quantum circuits.  Quantum, 5:391, January 2021.
  * [18] Maria Schuld, Ville Bergholm, Christian Gogolin, Josh Izaac, and Nathan Killoran.  Evaluating analytic gradients on quantum hardware.  Phys. Rev. A, 99:032331, Mar 2019.
  * [19] James Stokes, Josh Izaac, Nathan Killoran, and Giuseppe Carleo.  Quantum natural gradient.  Quantum, 4:269, May 2020.
  * [20] Javier Gil Vidal and Dirk Oliver Theis.  Calculus on parameterized quantum circuits.  Preprint arXiv:1812.06323, 2018.
  * [21] Rudin Walter.  Real and Complex Analysis.  McGraw-Hill, 1987.
  * [22] David Wierichs, Josh Izaac, Cody Wang, and Cedric Yen-Yu Lin.  General parameter-shift rules for quantum gradients.  Preprint, 2021.  arXiv:2107.12390.
  * [23] Xiao Yuan, Suguru Endo, Qi Zhao, Ying Li, and Simon C Benjamin.  Theory of variational quantum simulation.  Quantum, 3:191, 2019.