# On the energy landscape of symmetric quantum signal processing

Jiasu Wang  Department of Mathematics, University of California, Berkeley, CA
94720, USA  Yulong Dong  Department of Mathematics, University of California,
Berkeley, CA 94720, USA  Lin Lin  Department of Mathematics, University of
California, Berkeley, CA 94720, USA  Challenge Institute for Quantum
Computation, University of California, Berkeley, CA 94720, USA  Computational
Research Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720,
USA

###### Abstract

Symmetric quantum signal processing provides a parameterized representation of
a real polynomial, which can be translated into an efficient quantum circuit
for performing a wide range of computational tasks on quantum computers. For a
given polynomial, the parameters (called phase factors) can be obtained by
solving an optimization problem. However, the cost function is non-convex, and
has a very complex energy landscape with numerous global and local minima. It
is therefore surprising that the solution can be robustly obtained in
practice, starting from a fixed initial guess $\Phi^{0}$ that contains no
information of the input polynomial. To investigate this phenomenon, we first
explicitly characterize all the global minima of the cost function. We then
prove that one particular global minimum (called the maximal solution) belongs
to a neighborhood of $\Phi^{0}$, on which the cost function is strongly convex
under suitable conditions. This explains the aforementioned success of
optimization algorithms, and solves the open problem of finding phase factors
using only standard double precision arithmetic operations.

###### Contents

  1. 1 Introduction
    1. 1.1 Main results
    2. 1.2 Background and related works
    3. 1.3 Discussion and open questions
  2. 2 Preliminaries
    1. 2.1 Chebyshev polynomials
    2. 2.2 Notation
    3. 2.3 Quantum signal processing
    4. 2.4 Optimization based method for quantum signal processing
  3. 3 Symmetric quantum signal processing
    1. 3.1 Sign of the leading coefficient of $Q(x)$
    2. 3.2 Existence of symmetric phase factors
    3. 3.3 Proof of Theorem 1
    4. 3.4 Optimization based method for finding symmetric phase factors
  4. 4 Global minimizer
    1. 4.1 Characterizing complementary polynomials
    2. 4.2 Maximal solution
    3. 4.3 Relation to other constructive methods
  5. 5 Local convergence
    1. 5.1 Sketch of the proof
    2. 5.2 Recursion relation in terms of Chebyshev coefficients
    3. 5.3 Bounding Chebyshev coefficients of the maximal solution
    4. 5.4 Proof of Theorem 5
  6. 6 Local strong convexity
    1. 6.1 Hessian matrix at the maximal solution is positive definite
    2. 6.2 Lower bounding the distance to the initial point
    3. 6.3 Local strong convexity of the optimization problem
    4. 6.4 Projected gradient descend algorithm for symmetric phase-factor evaluation
  7. 7 Additional numerical results
    1. 7.1 Hessian matrix without symmetry constraint is singular at global optima
    2. 7.2 Existence of local minima
    3. 7.3 Numerical behavior near the maximal solution
  8. A Proof of Lemma 10
  9. B Proof of Lemma 11
  10. C Proof of Lemma 12

##  1 Introduction

Given a target polynomial $f\in\mathbb{R}[x]$ satisfying (1) $\deg(f)=d$, (2)
the parity of $f$ is $d\bmod 2$, (3) $\left\lVert
f\right\rVert_{\infty}:=\max_{x\in[-1,1]}\left\lvert f(x)\right\rvert<1$, the
problem of quantum signal processing (QSP) [14] is to find a set of parameters
(called phase factors) $\Phi:=(\phi_{0},\cdots,\phi_{d})\in[-\pi,\pi)^{d+1}$
so that

$f(x)=g(x,\Phi):=\mathrm{Re}[\braket{0}{U(x,\Phi)}{0}],\quad x\in[-1,1],$ (1)

with

$U(x,\Phi):=e^{\mathrm{i}\phi_{0}Z}e^{\mathrm{i}\arccos(x)X}e^{\mathrm{i}\phi_{1}Z}e^{\mathrm{i}\arccos(x)X}\cdots
e^{\mathrm{i}\phi_{d-1}Z}e^{\mathrm{i}\arccos(x)X}e^{\mathrm{i}\phi_{d}Z}.$
(2)

Here

$X:=\left(\begin{array}[]{cc}0&1\\\ 1&0\end{array}\right),\quad
Z:=\left(\begin{array}[]{cc}1&0\\\ 0&-1\end{array}\right)$ (3)

are Pauli $X$ and $Z$ matrices, respectively. $\braket{0}{U(x,\Phi)}{0}$
stands for the upper-left entry of the matrix $U(x,\Phi)\in\mathrm{SU}(2)$,
and $g(x,\Phi)$ defines a mapping from $[-1,1]$ to $\mathbb{R}$. When the
phase factors are restricted to be symmetric, i.e.

$\Phi=(\phi_{0},\phi_{1},\phi_{2},\ldots,\phi_{2},\phi_{1},\phi_{0})\in[-\pi,\pi)^{d+1},$
(4)

this is referred to as the symmetric quantum signal processing. The simplest
example is $\Phi=(0,\ldots,0)$. This gives
$U(x,\Phi)=e^{\mathrm{i}d\arccos(x)X}$ and
$g(x,\Phi)=\cos(d\arccos(x))=T_{d}(x)$, where $T_{d}$ is the Chebyshev
polynomial of the first kind of degree $d$.

Due to the parity constraint, the number of degrees of freedom in the target
polynomial $f(x)$ is $\widetilde{d}:=\lceil\frac{d+1}{2}\rceil$. Hence $f(x)$
is entirely determined by its values on $\widetilde{d}$ distinct points.
Throughout the paper, we choose these points to be
$x_{k}=\cos\left(\frac{2k-1}{4\widetilde{d}}\pi\right)$,
$k=1,...,\widetilde{d}$, i.e.positive nodes of the Chebyshev polynomial
$T_{2\widetilde{d}}(x)$. For any target polynomial $f(x)$ defined above, the
solution to Eq. 1 exists [14, 6], and Ref. [5] suggests that the solution can
be restricted to be symmetric. The existence of the solution implies that the
problem of symmetric quantum signal processing can be equivalently solved via
the following optimization problem

$\Phi^{*}=\mathop{\mathrm{argmin}}_{\begin{subarray}{c}\Phi\in[-\pi,\pi)^{d+1},\\\
\text{symmetric.}\end{subarray}}F(\Phi),\
F(\Phi):=\frac{1}{\widetilde{d}}\sum_{k=1}^{\widetilde{d}}\left\lvert
g(x_{k},\Phi)-f(x_{k})\right\rvert^{2},$ (5)

i.e.any solution $\Phi^{*}$ to Eq. 1 achieves the global minimum of the cost
function with $F(\Phi^{*})=0$, and vice versa.

However, the energy landscape of the cost function $F(\Phi)$ is very complex,
and has numerous global as well as local minima (see Section 7). This is
already the case with the symmetry constraint. (Without the symmetry
constraint, the number of variables is larger than the number of equations and
there should be an infinite number of global minima.) Starting from a random
initial guess, an optimization algorithm can easily be trapped at a local
minima already when $d$ is small111An early attempt showed that even finding
$\Phi$ with $d>30$ can be very costly [4, Apendix H.3].. It is therefore
surprising that starting from a special symmetric initial guess

$\Phi^{0}=(\pi/4,0,0,\ldots,0,0,\pi/4),$ (6)

at least one global minimum can be robustly identified using standard
unconstrained optimization algorithms even when $d$ is as large as $10,000$
[5], and the optimization method is free from being trapped by any local
minima. Direct calculation shows that $g(x,\Phi^{0})=0$, and therefore
$\Phi^{0}$ does not contain any _a priori_ information of the target
polynomial $f(x)$!

###  1.1 Main results

In this work, we provide a theory to characterize the energy landscape of the
cost function $F(\Phi)$, and to explain this unusual behavior of numerical
optimization. To our knowledge, this leads to the first provable algorithm for
solving the QSP problem without referring to extended precision arithmetic
operations [8].

Let the domain of the symmetric phase factors be

$D_{d}=\begin{cases}[-\frac{\pi}{2},\frac{\pi}{2})^{\frac{d}{2}}\times[-\pi,\pi)\times[-\frac{\pi}{2},\frac{\pi}{2})^{\frac{d}{2}},&d\mbox{
is even,}\\\ [-\frac{\pi}{2},\frac{\pi}{2})^{d+1},&d\mbox{ is odd.}\\\
\end{cases}$ (7)

Our first main result is the existence and uniqueness of symmetric phase
factors for a class of polynomial matrices in $\mathrm{SU}(2)$.

######  Theorem 1 (Existence and uniqueness of symmetric phase factors).

Consider any $P\in\mathbb{C}[x]$ and $Q\in\mathbb{R}[x]$ satisfying the
following conditions

  1. (1)

$\deg(P)=d$ and $\deg(Q)=d-1$.

  2. (2)

$P$ has parity $(d\bmod 2)$ and $Q$ has parity $(d-1\bmod 2)$.

  3. (3)

(Normalization condition) $\forall
x\in[-1,1]:|P(x)|^{2}+(1-x^{2})|Q(x)|^{2}=1$.

  4. (4)

If $d$ is odd, then the leading coefficient of $Q$ is positive.

There exists a unique set of symmetric phase factors
$\Phi:=(\phi_{0},\phi_{1},\cdots,\phi_{1},\phi_{0})\in D_{d}$ such that

$U(x,\Phi)=\begin{pmatrix}P(x)&\mathrm{i}Q(x)\sqrt{1-x^{2}}\\\
\mathrm{i}Q(x)\sqrt{1-x^{2}}&P^{*}(x)\end{pmatrix}.$ (8)

Theorem 1 implies that there is a bijection between the global minimizers of
Eq. 5 and all possible pairs of $(P(x),Q(x))$ satisfying the assumption in
Theorem 1 and $\mathrm{Re}[P](x)=f(x)$. The conditions (1), (2) for the target
polynomial $f$ are compatible with the first two requirements in Theorem 1.
The condition (3) on the maxnorm, i.e.$\left\lVert f\right\rVert_{\infty}<1$
is compatible with the normalization condition, which is itself a natural
condition due to the unitarity of $U(x,\Phi)$.

To simplify the discussion, we introduce the following definition of
admissible pair of polynomials associated with a target polynomial. Here
$P_{\mathrm{Re}}:=\mathrm{Re}[P],P_{\mathrm{Im}}:=\mathrm{Im}[P]$.

######  Definition 2 (Admissible pair of polynomials).

Let $f\in\mathbb{R}[x]$ be a target polynomial satisfying (1) $\deg(f)=d$, (2)
the parity of $f$ is $(d\bmod 2)$, (3) $\left\lVert f\right\rVert_{\infty}<1$.
Then $(P,Q)$ is an admissible pair of polynomials associated with $f$ if the
conditions (1)-(3) in Theorem 1 are satisfied, together with

  1. (4)

The leading coefficient of $Q$ is positive,

  2. (5)

$P_{\mathrm{Re}}(x)=f(x)$.

The polynomials $P_{\mathrm{Im}},Q\in\mathbb{R}[x]$ are also called
complementary polynomials to $f$.

Comparing Theorem 1 and Definition 2, the main modification is that the
leading coefficient of the complementary polynomial $Q$ is always restricted
to be positive. This allows us to unify the discussion of odd and even values
of $d$. The bijection relation can then be concisely formulated as follows.

######  Corollary 3 (Bijection between global minima and admissible pairs).

If $d$ is odd, there is a bijection between the global minima of Eq. 5 and all
admissible pairs $(P,Q)$.

If $d$ is even, there is a bijection between the global minima of Eq. 5 and
all pairs of polynomials $(P,\pm Q)$, where $(P,Q)$ is an admissible pair.

The proof of Theorem 1 is constructive. So given an admissible pair, we have
an algorithm to evaluate the symmetric phase factor. Theorem 4 explicitly
constructs all the admissible pairs. Together with Theorem 1, we have a
complete description of all global minima of the cost function in Eq. 5.

######  Theorem 4 (Construction of admissible pairs).

Given a target polynomial $f(x)$, all admissible pairs $(P,Q)$ must take the
following form,

$\begin{split}P_{\mathrm{Im}}\left(x\right)&=\sqrt{\alpha}\frac{e\left(x+\mathrm{i}\sqrt{1-x^{2}}\right)+e\left(x-\mathrm{i}\sqrt{1-x^{2}}\right)}{2},\\\
\quad
Q\left(x\right)&=\sqrt{\alpha}\frac{e\left(x+\mathrm{i}\sqrt{1-x^{2}}\right)-e\left(x-\mathrm{i}\sqrt{1-x^{2}}\right)}{2\mathrm{i}\sqrt{1-x^{2}}},\end{split}$
(9)

where

$\mathfrak{F}(z):=1-\left[f\left(\frac{z+z^{-1}}{2}\right)\right]^{2}=\alpha
e(z)e(z^{-1}),\quad e(z):=z^{-d}\prod_{i=1}^{2d}(z-r_{i}).$ (10)

The value of $\alpha\in\mathbb{C}$ depends on $\\{r_{i}\\}$.

Here $\\{r_{i}\\}_{i=1}^{2d}$ are roots of the _Laurent polynomial_
$\mathfrak{F}(z)$, which has $4d$ roots in total. By properly choosing the
roots (see a more detailed description in Remark 16), we can construct all the
admissible pairs. In particular, Theorem 4 shows that the number of global
minima is finite, but grows combinatorially with respect to $d$. The
finiteness of the number of global minima is also a direct consequence of the
compactness of the domain $D_{d}$.

Unfortunately, the procedure described above for finding phase factors is
numerically unstable. It requires extended precision arithmetic operations and
is therefore very expensive when $d$ is large (see detailed discussions in
Section 1.2). This is noticeably different from solving the optimization
problem in Eq. 5, which is numerically stable and can be readily performed
using standard double precision arithmetic operations.

Among the myriad of global minima, Theorem 4 allows us to identify a special
global minimum, which is obtained by choosing $\\{r_{i}\\}_{i=1}^{2d}$ to be
the roots of $\mathfrak{F}(z)$ within the unit disc. The unique symmetric
phase factor associated with this admissible pair is referred to as the
_maximal solution_ (the reason for the naming is technical and is explained in
Section 4.2).

The maximal solution enjoys many desirable properties. For any target
polynomial $f$ with $\left\lVert f\right\rVert_{\infty}\leq\frac{1}{2}$, the
maximal solution lies in the neighborhood of $\Phi^{0}$:

######  Theorem 5 (Distance between the maximal solution and $\Phi^{0}$).

Let $\Phi^{*}$ be the maximal solution for the target function $f(x)$. Denote
$\widetilde{\Phi}^{*}$ and $\widetilde{\Phi}^{0}$ as the corresponding reduced
phase factors of $\Phi^{*}$ and $\Phi^{0}$ respectively. If $\left\lVert
f(x)\right\rVert_{\infty}\leq\frac{1}{2}$, then

$\left\lVert\widetilde{\Phi}^{*}-\widetilde{\Phi}^{0}\right\rVert_{2}\leq\frac{\pi}{\sqrt{3}}\left\lVert
f(x)\right\rVert_{\infty}.$ (11)

Throughout this paper, given any symmetric phase factors $\Phi$ of length
$d+1$, reduced phase factors is referred to the left half part of $\Phi$,
namely $\left(\phi_{0},\phi_{1},\cdots,\phi_{\widetilde{d}-1}\right)$, and is
denoted as $\widetilde{\Phi}$. Recall that
$\widetilde{d}:=\left\lceil\frac{d+1}{2}\right\rceil$. Due to symmetry, $\Phi$
is fully determined by $\widetilde{\Phi}$. That is also why we are able to use
Euclidean distance between $\widetilde{\Phi}^{*}$ and $\widetilde{\Phi}^{0}$
instead of that between $\Phi^{*}$ and $\Phi^{0}$. Furthermore, we choose
$\widetilde{\Phi}$ as the free variables and by abuse of notation, denote cost
function $F(\Phi)$ still by $F(\widetilde{\Phi})$.

In particular, if the target polynomial is $f=0$, then the maximal solution is
the initial guess $\Phi^{0}$. We then prove that when $\left\lVert
f\right\rVert_{\infty}$ is sufficiently small, the maximal solution belongs to
a neighborhood of $\Phi^{0}$, on which the cost function $F(\widetilde{\Phi})$
is strongly convex, i.e.the Hessian matrix denoted by
$\mathrm{Hess}(\widetilde{\Phi})$ is positive definite.

######  Theorem 6 (Local strong convexity).

If the target polynomial satisfies $\left\lVert
f\right\rVert_{\infty}\leq\frac{\sqrt{3}}{20\pi\widetilde{d}}$, for any
symmetric phase factors $\Phi$ of length $d+1$ satisfying
$\left\lVert\widetilde{\Phi}-\widetilde{\Phi}^{0}\right\rVert_{2}\leq\frac{1}{20\widetilde{d}}$,
the following estimate holds:

$\frac{1}{4}\leq\lambda_{\min}\left(\mathrm{Hess}(\widetilde{\Phi})\right)\leq\lambda_{\max}\left(\mathrm{Hess}(\widetilde{\Phi})\right)\leq\frac{25}{4}.$
(12)

The local strong convexity result in Theorem 6 immediately implies that when
$\left\lVert f\right\rVert_{\infty}$ is sufficiently small, standard
optimization algorithms, such as the projected gradient method [1, 11], can
converge in the neighborhood of $\Phi^{0}$, without being trapped by any local
minima.

######  Corollary 7 (Convergence of projected gradient method).

If the target polynomial satisfies $\left\lVert
f\right\rVert_{\infty}\leq\frac{\sqrt{3}}{20\pi\widetilde{d}}$, starting from
$\Phi^{0}$, the projected gradient method with step size $t=\frac{1}{L}$
converges exponentially to the maximal solution $\Phi^{*}$, i.e.at the
$\ell$-th iteration

$\left\lVert\widetilde{\Phi}^{\ell}-\widetilde{\Phi}^{*}\right\rVert^{2}_{2}\leq
e^{-\frac{\sigma}{L}\ell}\left\lVert\widetilde{\Phi}^{0}-\widetilde{\Phi}^{*}\right\rVert^{2}_{2}.$
(13)

Here $\sigma=\frac{1}{4}$ and $L=\frac{25}{4}$.

Therefore when $\left\lVert
f\right\rVert_{\infty}\leq\frac{\sqrt{3}}{20\pi\widetilde{d}}$, in order to
reach precision $\epsilon$, the projected gradient method can be terminated
after $\mathcal{O}(\log(1/\epsilon))$ steps independent of the details of $f$.
Since each iteration is numerically stable, the maximal solution can be
readily obtained using standard double precision arithmetic operations in
practice.

###  1.2 Background and related works

The mapping in Eqs. 1 and 2 may seem a very peculiar way for encoding a
polynomial (and it is). However, such an SU(2) representation can be directly
translated into a quantum circuit, which is so far the most concise way for
performing eigenvalue transformation $f(A)$ (when $A$ is an Hermitian matrix)
[14], and singular value transformations $f^{\mathrm{SV}}(A)$ (when $A$ is a
general matrix) on quantum computers [6]. Many tasks in quantum computation
can be formulated using such transformations. When $f$ is not a polynomial, it
can be approximated by a polynomial with bounded precision. This strategy can
be used to unify a large class of quantum algorithms including Grover’s search
and quantum phase estimation [17, 16], and to perform a wide range of
computational tasks, such as solving linear system of equations [17, 13],
eigenvalue problems [12], Hamiltonian simulation [14, 6] etc.

After the initial efforts [14, 4], significant progress has been made by using
direct methods to obtain phase factors [6, 8]. These methods are based on
finding roots of high degree polynomials to high precision and are not
numerically stable. Specifically, these algorithms require
$\mathcal{O}(d\log(d/\epsilon))$ bits of precision, where $d$ is the degree of
$f(x)$ and $\epsilon$ is the target accuracy [8]. It is worth mentioning that
the extended precision needed in these algorithms is not an artifact of the
proof technique. For instance, for $d\approx 500$, the number of bits needed
to represent each floating point number can be as large as $1000\sim 2000$
[5], which is much larger than the $64$ bits provided by standard double
precision floating point format.

Recently, alternative methods have been developed to tackle the QSP problem
[3, 5] using only standard double precision arithmetic operations (for large
$d$). The algorithm in [3] is also a heuristic algorithm and is a variation of
that in [6]. The optimization based algorithm in [5] is different from all
other proposals so far, and this work establishes the effectiveness of the
optimization based approach under the condition that $\left\lVert
f\right\rVert_{\infty}$ is sufficiently small. This solves the open problem of
finding phase factors using numerically stable algorithms, and justifies the
usage of standard double precision arithmetic operations in practice.

###  1.3 Discussion and open questions

Without the symmetry constraint, the optimization problem in Eq. 5 has $d$
degrees of freedom. This is larger than the number of degrees of freedom of
target function $f$, which is only $\widetilde{d}$ due to the parity
constraint. As a result, the Hessian matrix is always singular, preventing us
from proving the exponential convergence of optimization methods [18].

In Theorem 6, the local strong convexity is only proved when $\left\lVert
f\right\rVert_{\infty}=\mathcal{O}(d^{-1})$. Therefore for a given target
polynomial, in order to apply Corollary 7, we can only find the phase factors
for $cf(x)$, where $c=\mathcal{O}(d^{-1})$ is a scaling factor. Such a scaling
factor is undesirable and leads to suboptimal query complexities of quantum
algorithms. On the other hand, numerical results suggest that optimization
algorithms can robustly converge to a global minimum when $\left\lVert
f\right\rVert_{\infty}$ is smaller than a constant (less than $1$) that is
independent of $d$. Note that the estimate in Theorem 5 is already independent
of $d$. Hence the gap is mainly due to Theorem 6, where the eigenvalue
estimate of the Hessian matrix is obtained via the estimate of each matrix
entry, which introduces the $d$-dependence. It is a natural open question to
prove the strong convexity of the cost function near $\Phi^{0}$ where
$\left\lVert f\right\rVert_{\infty}$ can be independent of $d$.

The proof of Corollary 7 requires the projection of $\widetilde{\Phi}^{\ell}$
to the domain
$\\{\widetilde{\Phi}:\left\lVert\widetilde{\Phi}-\widetilde{\Phi}^{0}\right\rVert_{2}\leq\frac{1}{20\widetilde{d}}\\}$
at each iteration. This is because the steepest descent method may potentially
overshoot and deviate away from the locally strong convex region. Numerical
observation shows that unconstrained optimizers, such as quasi-Newton methods
can robustly achieve the global minimum, which is beyond our current
theoretical analysis.

Finally, since $\Phi^{0}$ is only one of the solution when $f=0$, we may ask
what if we run the projected gradient method starting from other solutions
when $f=0$, and to obtain optimal phase factors for more general target
polynomials. These solutions correspond to non-maximal solutions. Numerical
observations show that the convergence rate towards such non-maximal solutions
can be significantly slower than that towards the maximal solution (see Figs.
4 and 5). The analysis of such behavior requires the generalization of
estimates of in Theorems 5 and 6 to non-maximal solutions. These will be our
future works.

### Acknowledgments:

This work was partially supported by the NSF Quantum Leap Challenge Institute
(QLCI) program through grant number OMA-2016245 (J.W.,Y.D.), by Department of
Energy under Grant No. DE-SC0017867 and No. DE-AC02-05CH11231 (L.L.). L.L. is
a Simons Investigator.

##  2 Preliminaries

###  2.1 Chebyshev polynomials

We first collect some basic facts of Chebyshev polynomials, which is used to
formulate the optimization problem in Eq. 5 and is instrumental in Sections 5
and 6. Given $x\in[-1,1]$, the Chebyshev polynomial of the first kind is
$T_{n}(x)=\cos\left(n\arccos{(x)}\right)$, and that of the second kind is
given by $U_{n}(x)\sin(x)=\sin\left((n+1)\arccos{(x)}\right)$.

A useful fact is that both $T_{n}$ and $U_{n}$ form a sequence of orthogonal
polynomials. The Chebyshev polynomial of the first kind $T_{n}$ are orthogonal
with respect to the weight $\frac{1}{\sqrt{1-x^{2}}}$ over $[-1,1]$, and that
of the second kind $U_{n}$ are orthogonal with respect to the weight
$\sqrt{1-x^{2}}$ over $[-1,1]$,

$\begin{split}&\frac{1}{\pi}\int_{-1}^{1}T_{n}(x)T_{m}(x)\frac{1}{\sqrt{1-x^{2}}}\mathrm{d}x=\frac{1}{2}\delta_{nm}(1+\delta_{n0}),\\\
&\frac{1}{\pi}\int_{-1}^{1}U_{n}(x)U_{m}(x)\sqrt{1-x^{2}}\mathrm{d}x=\frac{1}{2}\delta_{nm}.\end{split}$

We refer the induced norms

$\left\lVert
f\right\rVert_{T}^{2}:=\frac{1}{\pi}\int_{-1}^{1}f^{2}(x)\frac{1}{\sqrt{1-x^{2}}}\mathrm{d}x,$

and

$\left\lVert
f\right\rVert_{U}^{2}:=\frac{1}{\pi}\int_{-1}^{1}f^{2}(x)\sqrt{1-x^{2}}\mathrm{d}x$

to as the Chebyshev norm of the first kind and the second kind respectively.
Furthermore, any function $h$ with $\left\lVert h\right\rVert_{T}<\infty$ can
be uniquely expressed as a series of Chebyshev polynomials of the first kind,

$h(x)=\sum_{n=0}^{\infty}c_{n}T_{n},\qquad
c_{n}=\frac{2}{\pi(1+\delta_{n0})}\int_{-1}^{1}h(x)T_{n}(x)\frac{1}{\sqrt{1-x^{2}}}\mathrm{d}x.$
(14)

Similarly, any function $h$ with $\left\lVert h\right\rVert_{U}<\infty$ can be
uniquely expressed as a series of Chebyshev polynomials of the second kind,

$h(x)=\sum_{n=0}^{\infty}c_{n}U_{n},\qquad
c_{n}=\frac{2}{\pi}\int_{-1}^{1}h(x)T_{n}(x)\sqrt{1-x^{2}}\mathrm{d}x.$ (15)

Apart from the weighted orthogonality, Chebyshev polynomials of the first kind
$T_{n}$ satisfies the discrete orthogonality condition:

$\sum_{j=1}^{2\widetilde{d}}T_{n}(x_{j})T_{m}(x_{j})=\begin{cases}0&n\neq m\\\
2\widetilde{d}&n=m=0\\\ \widetilde{d}&n=m\neq 0\end{cases},$ (16)

where $x_{j}=\cos\left(\pi\frac{2j-1}{4\widetilde{d}}\right)$ are the nodes of
$T_{2\widetilde{d}}$ and $\max(n,m)\leq 2\widetilde{d}$. This set of points is
also referred to as Chebyshev nodes. Furthermore, if $n$ and $m$ have the same
parity which implies $T_{n}(-x)T_{m}(-x)=T_{n}(x)T_{m}(x)$, we have

$\sum_{j=1}^{\widetilde{d}}T_{n}(x_{j})T_{m}(x_{j})=\frac{1}{2}\sum_{j=1}^{2\widetilde{d}}T_{n}(x_{j})T_{m}(x_{j})=\begin{cases}0&n\neq
m\\\ \widetilde{d}&n=m=0\\\ \frac{\widetilde{d}}{2}&n=m\neq 0\end{cases}.$
(17)

###  2.2 Notation

Given a positive integer $n$, $\mathbb{K}_{n}[x]$ denotes the set of all
polynomials of degree at most $n$ with variable $x$, and coefficients taken in
the field $\mathbb{K}$. Specifically, we are only interested in the cases
where $\mathbb{K}$ is either $\mathbb{R}$ or $\mathbb{C}$. Throughout the
paper, without otherwise noted, we assume that the length of the set of phase
factors is $d+1$. We also assume that a pair of polynomials
$(P,Q)\in\mathbb{C}_{d}[x]\times\mathbb{C}_{d-1}[x]$ satisfies the conditions
in Theorem 1. Given a polynomial $P\in\mathbb{C}[x]$, we define
$P_{\mathrm{Re}}(x):=\mathrm{Re}[P(x)]$ and
$P_{\mathrm{Im}}(x):=\mathrm{Im}\left[P(x)\right]$. The real polynomial
$P_{\mathrm{Re}}(x)$ is usually denoted by $f(x)$ when it refers to the target
of the optimization.

The QSP problem involves some delicate relations of polynomial coefficients.
Throughout the paper, we will interchangeably expand a polynomial using the
monomial basis and the Chebyshev basis, i.e.

$f(x)=\sum_{j=0}^{d}\mathbf{f}_{j}x^{j}=\sum_{j=0}^{d}f_{j}T_{j}(x),\quad\mathbf{f}_{j},f_{j}\in\mathbb{R}.$
(18)

Similarly,

$\begin{split}P(x)&=\sum_{j=0}^{d}\mathbf{p}_{j}x^{j}=\sum_{j=0}^{d}p_{j}T_{j}(x),\quad\mathbf{p}_{j},p_{j}\in\mathbb{C},\\\
P_{\mathrm{Im}}(x)&=\sum_{j=0}^{d}\mathbf{s}_{j}x^{j}=\sum_{j=0}^{d}s_{j}T_{j}(x),\quad\mathbf{s}_{j},s_{j}\in\mathbb{R},\\\
Q(x)&=\sum_{j=0}^{d-1}\mathbf{q}_{j}x^{j}=\sum_{j=0}^{d-1}q_{j}U_{j}(x)\quad\mathbf{q}_{j},q_{j}\in\mathbb{C}.\end{split}$
(19)

For convenience, we set the values of
$\mathbf{f}_{j},f_{j},\mathbf{p}_{j},p_{j},\mathbf{s}_{j},s_{j},\mathbf{q}_{j},q_{j}$
with a negative integer index to zero.

The leading coefficients in the monomial and Chebyshev basis satisfy the
relation

$\quad\mathbf{f}_{d}=2^{d-1}f_{d},\quad\mathbf{p}_{d}=2^{d-1}p_{d},\quad\mathbf{s}_{d}=2^{d-1}s_{d},\quad\mathbf{q}_{d}=2^{d}q_{d}.$
(20)

The set $[n]:=\\{0,1,\cdots,n-1\\}$ is referred to as the index set generated
by a positive integer $n$. The row and column indices of a $n$-by-$n$ matrix
run from $0$ to $n-1$, namely in the index set $[n]$. For a matrix
$A\in\mathbb{C}^{m\times n}$, the transpose, Hermitian conjugate and complex
conjugate are denoted by $A^{\top}$, $A^{{\dagger}}$, $A^{*}$, respectively.
The same notations are also used for the operations on a vector. We use the
convention in quantum computing to write the basis vectors in
$\mathbb{C}^{2}$. It is equivalent to the following identification

$\ket{0}=\left(\begin{array}[]{c}1\\\ 0\end{array}\right),\
\ket{1}=\left(\begin{array}[]{c}0\\\ 1\end{array}\right),\
\bra{0}=\ket{0}^{\dagger}=\left(1,0\right),\text{ and
}\bra{1}=\ket{1}^{\dagger}=\left(0,1\right).$

Given a $2$-by-$2$ matrix $A=(a_{ij})$, we have $\braket{i}{A}{j}=a_{ij}$ for
any $i,j\in\\{0,1\\}$.

For a matrix $A\in\mathbb{C}^{m\times n}$, we denote the operator norm
$\left\lVert\cdot\right\rVert_{2}$ induced by the vector $\ell_{2}$ norm as

$\left\lVert A\right\rVert_{2}:=\sup_{x\neq 0}\frac{\left\lVert
Ax\right\rVert_{2}}{\left\lVert x\right\rVert_{2}},$ (21)

and denote the Frobenius norm $\left\lVert\cdot\right\rVert_{F}$ as

$\left\lVert A\right\rVert_{F}:=\sqrt{\sum_{i\in[m],j\in[n]}\left\lvert
A_{ij}\right\rvert^{2}}.$ (22)

We use $\sigma_{k}(A)$ to denote its $(k+1)$-th largest singular value, i.e.,

$0\leq\sigma_{\min(m,n)-1}(A)\leq\cdots\leq\sigma_{2}(A)\leq\sigma_{0}(A).$

We define $\sigma_{\min}(A):=\sigma_{\min(m,n)-1}(A)$ and
$\sigma_{\max}(A):=\sigma_{0}(A)$. Here, the index $k$ starts from $0$ which
agrees with the convention of matrix index.

For a Hermitian matrix $A\in\mathbb{C}^{n\times n}$, we use the notation
$\lambda_{k}(A)$ to designate the $(k+1)$-th largest eigenvalue, i.e.,

$\lambda_{n-1}(A)\leq\cdots\leq\lambda_{2}(A)\leq\lambda_{0}(A).$

We define $\lambda_{\min}(A):=\lambda_{n-1}(A)$ and
$\lambda_{\max}(A):=\lambda_{0}(A)$.

###  2.3 Quantum signal processing

In the absence of symmetry constraint, the existence of phase factors for a
unitary representation of a pair of polynomials $(P,Q)$ is given by the
following theorem.

######  Theorem 8 (Quantum signal processing [6, Theorem 4]).

Let $d\in\mathbb{N}$. There exists a set of phase factors
$\Phi:=(\phi_{0},\cdots,\phi_{d})\in\mathbb{R}^{d+1}$ such that

$U(x,\Phi)=e^{\mathrm{i}\phi_{0}Z}\prod_{j=1}^{d}\left(W(x)e^{\mathrm{i}\phi_{j}Z}\right)=\left(\begin{array}[]{cc}P(x)&\mathrm{i}Q(x)\sqrt{1-x^{2}}\\\
\mathrm{i}Q^{*}(x)\sqrt{1-x^{2}}&P^{*}(x)\end{array}\right),$ (23)

where

$W(x)=e^{\mathrm{i}\arccos(x)X}=\left(\begin{array}[]{cc}{x}&{\mathrm{i}\sqrt{1-x^{2}}}\\\
{\mathrm{i}\sqrt{1-x^{2}}}&{x}\end{array}\right),$

if and only if $P,Q\in\mathbb{C}[x]$ satisfy that

  1. (1)

$\deg(P)\leq d,\deg(Q)\leq d-1$,

  2. (2)

$P(x)$ has parity $(d\mod 2)$ and $Q(x)$ has parity $(d-1\mod 2)$,

  3. (3)

$|P(x)|^{2}+(1-x^{2})|Q(x)|^{2}=1,\forall x\in[-1,1]$.

When the real component of $P$ is of interest only, the conditions for the
existence of phase factors can be further simplified.

######  Corollary 9 (Real quantum signal processing [6, Corollary 5]).

Let $f\in\mathbb{R}[x]$ be a degree-d polynomial for some $d\geq 1$ such that

  * •

$f(x)$ has parity $(d\mod 2)$,

  * •

$\left\lvert f(x)\right\rvert<1,\forall x\in[-1,1]$.

Then there exists some $P,Q\in\mathbb{C}[x]$ satisfying properties
LABEL:{itm:1}, LABEL:{itm:2} and LABEL:{itm:3} of Theorem 8 such that
$f(x)=\mathrm{Re}[P(x)]$.

###  2.4 Optimization based method for quantum signal processing

The existence result in Corollary 9 enables the following optimization based
strategy to find phase factors. Given the target polynomial
$f(x)\in\mathbb{R}[x]$, we may minimize the $\mathcal{L}_{2}$ distance between
$g(x,\Phi)$ and target function $f(x)$, i.e.

$L(\Phi)=\left(\int_{-1}^{1}\left[g(x,\Phi)-f(x)\right]^{2}\mathrm{d}x\right)^{\frac{1}{2}}.$
(24)

The desired phase factors $\Phi^{*}$ satisfy $L(\Phi^{*})=0$.

Since both $g(x,\Phi)$ and $f(\Phi)$ are polynomials of degree $d$ with
definite parity, the number of degrees of freedom is
$\widetilde{d}:=\lceil\frac{d+1}{2}\rceil$. These polynomials are determined
by their values at $x_{k}=\cos\left(\frac{2k-1}{4\widetilde{d}}\pi\right)$,
$k=1,...,\widetilde{d}$, which are positive Chebyshev nodes of
$T_{2\widetilde{d}}(x)$. This leads to the following optimization problem,

$\Phi^{*}=\mathop{\mathrm{argmin}}_{\Phi\in[-\pi,\pi)^{d+1}}F(\Phi),\
F(\Phi):=\frac{1}{\widetilde{d}}\sum_{k=1}^{\widetilde{d}}\left\lvert
g(x_{k},\Phi)-f(x_{k})\right\rvert^{2}.$ (25)

However, the energy landscape of $F(\Phi)$ is very complex. Starting from a
random initial guess, an optimization method often converges to a local minima
where $F(\Phi)>0$.

##  3 Symmetric quantum signal processing

In this section, we present some basic properties of symmetric quantum signal
processing. We first explain a constraint of the sign of the leading
coefficient of $Q(x)$ when $d$ is odd. Then we introduce a constructive method
for finding phase factors given suitable polynomials $\left(P,Q\right)$. In
the end, we prove Theorem 1, which establishes a bijection between the global
minimizers of Eq. 5 and all suitable complementary polynomials
$\left(P_{\mathrm{Im}},Q\right)$.

###  3.1 Sign of the leading coefficient of $Q(x)$

We first derive an explicit expression of the coefficients of $T_{d}$ and
$U_{d-1}$ in the Chebyshev expansion of $P$ and $Q$ respectively.

######  Lemma 10.

Let $\Phi\in\mathbb{R}^{d+1}$ be a set of phase factors, $P(x)$ and $Q(x)$ be
defined by Eq. 23. Considering the Chebyshev expansion of $P$ and $Q$ in Eq.
19, one has

$p_{d}=e^{\mathrm{i}\left(\phi_{0}+\phi_{d}\right)}\prod_{j=1}^{d-1}\cos\left(\phi_{j}\right),\quad
q_{d-1}=e^{\mathrm{i}\left(\phi_{0}-\phi_{d}\right)}\prod_{j=1}^{d-1}\cos\left(\phi_{j}\right).$
(26)

###### Proof.

See Appendix A. ∎

After imposing the symmetry constraint on phase factors, we find that there is
an additional constraint on the leading coefficient of $Q$ when $d$ is odd. It
originates from the observation that when $\Phi$ is symmetric, one has
$\phi_{0}=\phi_{d}$ and then the second equality in Eq. 26 becomes

$q_{d-1}=\prod_{j=1}^{d-1}\cos\left(\phi_{j}\right)=\prod_{j=1}^{\widetilde{d}-1}\cos^{2}\left(\phi_{j}\right)\geq
0.$ (27)

Note that $\deg(Q)$ may be less than $d-1$. The more precise statement is
given in Lemma 11.

######  Lemma 11.

Let $\Phi\in\mathbb{R}^{d+1}$ be a set of phase factors, $P(x)$ and $Q(x)$ be
defined by Eq. 23. If $\Phi$ is symmetric and $d$ is odd, then the leading
coefficient of $Q$ in the Chebyshev basis has the same sign as
$\left(-1\right)^{\frac{d-1-\deg(Q)}{2}}$.

###### Proof.

See Appendix B. ∎

As a remark, the statement in Lemma 11 also holds with the Chebyshev basis
replaced by the monomial basis, which is an consequence of directly applying
the relation in Eq. 20.

###  3.2 Existence of symmetric phase factors

We present a useful lemma, which evaluates the outmost symmetric phase factor,
and reduces the polynomial degree by $2$. It also connects the leading
coefficients in the monomial basis before and after the reduction. The result
of Lemma 12 will be used throughout the paper. To clarify, we only consider
the expansion of polynomial in the monomial basis in rest of this section.

######  Lemma 12.

Suppose that $P(x)\in\mathbb{C}[x],Q(x)\in\mathbb{R}[x]$ satisfy the
conditions in Theorem 1 and $d\geq 2$. Choose

$e^{2\mathrm{i}\phi_{0}}=\frac{\mathbf{q}_{d-1}}{\mathbf{p}_{d}^{*}}$ (28)

and define

$\begin{split}&\begin{pmatrix}P^{(1)}(x)&\mathrm{i}Q^{(1)}(x)\sqrt{1-x^{2}}\\\
\mathrm{i}Q^{(1)}(x)\sqrt{1-x^{2}}&\left(P^{(1)}\right)^{*}(x)\end{pmatrix}\\\
&:=W^{-1}e^{-\mathrm{i}\phi_{0}Z}\begin{pmatrix}P(x)&\mathrm{i}Q(x)\sqrt{1-x^{2}}\\\
\mathrm{i}Q(x)\sqrt{1-x^{2}}&P^{*}(x)\end{pmatrix}e^{-\mathrm{i}\phi_{0}Z}W^{-1}.\end{split}$
(29)

Consider the expansion of $P$, $Q$, $P^{(1)}$ and $Q^{(1)}$ in the monomial
basis,

$\begin{split}P(x)&=\sum_{i=0}^{d}\mathbf{p}_{i}x^{i},\quad\quad
Q(x)=\sum_{j=0}^{d-1}\mathbf{q}_{j}x^{j}.\\\
P^{(1)}(x)&=\sum_{i=0}^{d+2}\mathbf{p}_{i}^{(1)}x^{i},\quad\quad
Q^{(1)}(x)=\sum_{j=0}^{d+1}\mathbf{q}_{j}^{(1)}x^{j}.\end{split}$ (30)

Then $Q^{(1)}(x)\in\mathbb{R}[x]$ and $P^{(1)},Q^{(1)}$ satisfy the conditions
(1)-(3) in Theorem 1 with $d$ replaced by $d-2$.

Furthermore, if $d=2$, then

$P^{(1)}=-\mathbf{p}_{0}^{*}e^{2\mathrm{i}\phi_{0}}=-\left(\frac{\mathbf{p}_{0}}{\mathbf{p}_{2}}\right)^{*}\mathbf{q}_{1},\quad
Q^{(1)}=0.$ (31)

If $d\geq 3$, then

$\begin{split}&\mathbf{p}_{d-2}^{(1)}=\frac{\mathbf{q}_{d-1}}{4}-\frac{1}{\mathbf{q}_{d-1}}\left\lvert\mathrm{Im}[\mathbf{p}_{d-2}e^{-2\mathrm{i}\phi_{0}}]\right\rvert^{2}+\mathrm{i}\mathrm{Im}[\mathbf{p}_{d-2}e^{-2\mathrm{i}\phi_{0}}],\\\
&\mathbf{q}_{d-3}^{(1)}=\frac{\mathbf{q}_{d-1}}{4}+\frac{1}{\mathbf{q}_{d-1}}\left\lvert\mathrm{Im}[\mathbf{p}_{d-2}e^{-2\mathrm{i}\phi_{0}}]\right\rvert^{2}.\end{split}$
(32)

And $P^{(1)},Q^{(1)}$ also satisfy the condition (4) in Theorem 1.

###### Proof.

See Appendix C. ∎

For convenience we denote $\left(P^{(0)},Q^{(0)}\right)=(P,Q)$. Inspired by
Lemma 12, we choose $\phi_{\ell}$ to satisfy

$e^{2\mathrm{i}\phi_{\ell}}=\frac{\mathbf{p}^{(\ell)}_{d-2\ell}}{\mathbf{q}^{(\ell)}_{d-1-2\ell}}$
(33)

and repeat the construction

$\begin{split}&\begin{pmatrix}P^{(\ell+1)}(x)&\mathrm{i}Q^{(\ell+1)}(x)\sqrt{1-x^{2}}\\\
\mathrm{i}Q^{(\ell+1)}(x)\sqrt{1-x^{2}}&\left(P^{(\ell+1)}\right)^{*}(x)\end{pmatrix}\\\
&:=W^{-1}e^{-\mathrm{i}\phi_{\ell}Z}\begin{pmatrix}P^{(\ell)}(x)&\mathrm{i}Q^{(\ell)}(x)\sqrt{1-x^{2}}\\\
\mathrm{i}Q^{(\ell)}(x)\sqrt{1-x^{2}}&\left(P^{(\ell)}\right)^{*}(x)\end{pmatrix}e^{-\mathrm{i}\phi_{\ell}Z}W^{-1}.\end{split}$
(34)

as long as $\mathrm{deg}\left(P^{(\ell)}\right)\geq 2$. Notice that

$\deg\left(P^{(\ell)}\right)=\deg\left(P^{(\ell-1)}\right)-2=\deg\left(P^{(0)}\right)-2\ell.$

Thus, the iteration stops when $\ell=\widetilde{d}-1$. Moreover,
$\deg(P^{(\widetilde{d}-1)})=0$ if $d$ is even, and
$\deg(P^{(\widetilde{d}-1)})=1$ if $d$ is odd. For convenience, we introduce
the notation of $\hat{d}$,

$\hat{d}:=\begin{cases}\widetilde{d}-1&\text{ if }d\text{ is odd},\\\
\widetilde{d}-2&\text{ if }d\text{ is even}.\end{cases}$ (35)

######  Lemma 13.

Suppose that $P(x),Q(x)$ satisfy the conditions in Theorem 1. Then there
exists a set of symmetric phase factor
$\Phi:=(\phi_{0},\phi_{1},\cdots,\phi_{1},\phi_{0})\in D_{d}$ such that Eq. 23
holds, where $\phi_{\ell}$ satisfies

$e^{2\mathrm{i}\phi_{\ell}}=\frac{\mathbf{p}^{(\ell)}_{d-2\ell}}{\mathbf{q}^{(\ell)}_{d-1-2\ell}},\quad\forall
0\leq\ell\leq\hat{d}.$ (36)

In particular, if $d$ is even, $\phi_{\widetilde{d}-1}$ satisfies

$e^{\mathrm{i}\phi_{\widetilde{d}-1}}=\mathbf{p}^{(\widetilde{d}-1)}_{0}.$
(37)

Here $\mathbf{p}^{(\ell)}_{d-2\ell}$ and $\mathbf{q}^{(\ell)}_{d-1-2\ell}$ are
the leading coefficient of the expansion of $P^{(\ell)}$ and $Q^{(\ell)}$ in
the monomial basis respectively. We use the notation
$\left(P^{(0)},Q^{(0)}\right)=\left(P,Q\right)$ and define
$(P^{(\ell)},Q^{(\ell)})$ by Eq. 34 for $1\leq\ell\leq\widetilde{d}-1$.

###### Proof.

By repeating the construction in Eq. 34, we obtain a sequence of polynomials
$\left(P^{(\ell)},Q^{(\ell)}\right)$ where $0\leq\ell\leq\widetilde{d}-1$. Eq.
36 determines the value of $\phi_{\ell}$ over $D_{d}$ for
$0\leq\ell\leq\widetilde{d}-2$.

When $d$ is odd, $\deg\left(P^{(\widetilde{d}-1)}\right)=1$ and
$\deg\left(Q^{(\widetilde{d}-1)}\right)=0$. Lemma 12 implies that the leading
coefficient of $Q^{(\widetilde{d}-1)}$ is positive as well. By the
normalization condition, $P^{(\widetilde{d}-1)}$ must be of the form $cx$
where $\left\lvert c\right\rvert=1$ and $Q^{(\widetilde{d}-1)}$ must be 1.
Choose $\phi_{\widetilde{d}-1}\in[-\frac{\pi}{2},\frac{\pi}{2})$ according to
Eq. 36 and one can check

$\begin{pmatrix}P^{(\widetilde{d}-1)}(x)&\mathrm{i}Q^{(\widetilde{d}-1)}(x)\sqrt{1-x^{2}}\\\
\mathrm{i}Q^{(\widetilde{d}-1)}(x)\sqrt{1-x^{2}}&\left(P^{(\widetilde{d}-1)}\right)^{*}(x)\end{pmatrix}=e^{\mathrm{i}\phi_{\widetilde{d}-1}Z}We^{\mathrm{i}\phi_{\widetilde{d}-1}Z}.$

In this way,
$\Phi=(\phi_{0},\cdots,\phi_{\widetilde{d}-1},\phi_{\widetilde{d}-1},\cdots,\phi_{0})$
satisfies the requirements.

When $d$ is even, we know that $P^{(\widetilde{d}-1)}$ is a constant of
absolute value $1$ and $Q^{(\widetilde{d}-1)}=0$, i.e.

$\begin{pmatrix}P^{(\widetilde{d}-1)}(x)&\mathrm{i}Q^{(\widetilde{d}-1)}(x)\sqrt{1-x^{2}}\\\
\mathrm{i}Q^{(\widetilde{d}-1)}(x)\sqrt{1-x^{2}}&\left(P^{(\widetilde{d}-1)}\right)^{*}(x)\end{pmatrix}=e^{\mathrm{i}\phi_{\widetilde{d}-1}Z},$

where $\phi_{\widetilde{d}-1}$ satisfies Eq. 37. Furthermore, we may choose
$\phi_{\widetilde{d}-1}\in[-\pi,\pi)$ and
$\Phi=(\phi_{0},\cdots,\phi_{\widetilde{d}-2},\phi_{\widetilde{d}-1},\phi_{\widetilde{d}-2},\cdots,\phi_{0})$
satisfies the requirements. ∎

###  3.3 Proof of Theorem 1

Now we are ready to prove Theorem 1.

###### Proof of Theorem 1.

The existence is given by Lemma 13. We only need to show the uniqueness. If
there exists another set of symmetric phase factors
$\Psi:=(\psi_{0},\psi_{1},\cdots,\psi_{1},\psi_{0})\in D_{d}$ satisfying the
requirement, we may consider

$\begin{split}&\begin{pmatrix}P^{(1)}_{\Psi}(x)&\mathrm{i}Q^{(1)}_{\Psi}(x)\sqrt{1-x^{2}}\\\
\mathrm{i}Q^{(1)}_{\Psi}(x)\sqrt{1-x^{2}}&\left(P^{(1)}_{\Psi}\right)^{*}(x)\end{pmatrix}\\\
&:=W^{-1}e^{-\mathrm{i}\psi_{0}Z}\begin{pmatrix}P(x)&\mathrm{i}Q(x)\sqrt{1-x^{2}}\\\
\mathrm{i}Q(x)\sqrt{1-x^{2}}&P^{*}(x)\end{pmatrix}e^{-\mathrm{i}\psi_{0}Z}W^{-1}.\end{split}$

Since $P^{(1)}_{\Psi}(x)$ can be obtained by a set of symmetric phase factors
of length $d-2$, the coefficient of $x^{d+2}$ of $P^{(1)}_{\Psi}$ is zero,
i.e.

$\mathbf{p}^{(1)}_{d+2}=-2\mathbf{q}_{d-1}+\mathbf{p}_{d}^{*}e^{2\mathrm{i}\psi_{0}}+\mathbf{p}_{d}e^{-2\mathrm{i}\psi_{0}}=0.$

Together with $\left\lvert\mathbf{p}_{d}\right\rvert=\mathbf{q}_{d-1}$, it
leads to $e^{2\mathrm{i}\psi_{0}}=\frac{\mathbf{q}_{d-1}}{\mathbf{p}_{d}^{*}}$
and this equation has only one solution over $[-\frac{\pi}{2},\frac{\pi}{2})$,
hence $\psi_{0}=\phi_{0}$. Then $\Phi=\Psi$ follows inductively. ∎

###  3.4 Optimization based method for finding symmetric phase factors

Ref. [5] suggest imposing the symmetry constraint on the phase factors and
study the optimization problem Eq. 5 instead. Throughout this paper, we choose
the reduced phase factors of $\Phi$,
$(\phi_{0},\phi_{1},\cdots,\phi_{\widetilde{d}-1})$, as free variables. To
clarify, without otherwise noted, $\Phi$ denotes the phase factors under the
symmetry constraint and $\widetilde{\Phi}$ denotes the corresponding reduced
phase factors. By abuse of notation, we denote cost function $F(\Phi)$ still
by $F(\widetilde{\Phi})$ and denote function $g(x,\Phi)$ by
$g(x,\widetilde{\Phi})$. We use subscripts to represent the partial derivative
of any function of the set of reduced phase factors with respect to the
corresponding phase factor, for example,
$F_{i}(\widetilde{\Phi}):=\partial_{\phi_{i}}F(\widetilde{\Phi})$. These
notations and conventions are also suitable for LABEL:{sec:local_convergence}
and LABEL:{sec:local_strong_convexity}.

The Hessian matrix of the problem in Eq. 5 can be written element-wisely as

$\begin{split}F_{ij}(\widetilde{\Phi})&=\frac{2}{\widetilde{d}}\sum_{k=1}^{\widetilde{d}}g_{i}(x_{k},\widetilde{\Phi})g_{j}(x_{k},\widetilde{\Phi})-\frac{2}{\widetilde{d}}\sum_{k=1}^{\widetilde{d}}[g(x_{k},\widetilde{\Phi})-f(x_{k})]g_{ij}(x_{k},\widetilde{\Phi})\\\
&=\frac{2}{\widetilde{d}}(A^{\top}A)_{ij}-R_{ij},\end{split}$ (38)

where
$A(\widetilde{\Phi}):=[g_{j}(x_{i},\widetilde{\Phi})]\in\mathbb{R}^{\widetilde{d}\times\widetilde{d}}$
is the Jacobian matrix of
$\widetilde{\Phi}\mapsto\left(g(x_{1},\widetilde{\Phi}),\dots,g(x_{\widetilde{d}},\widetilde{\Phi})\right)^{\top}$,
and $R$ is a real $\widetilde{d}\times\widetilde{d}$ matrix closely related to
the residual with entry

$R_{ij}(\widetilde{\Phi}):=\frac{2}{\widetilde{d}}\sum_{k=1}^{\widetilde{d}}[g(x_{k},\widetilde{\Phi})-f(x_{k})]g_{ij}(x_{k},\widetilde{\Phi}).$

Due to the existence of phase factors stated in Theorem 1, at the optimal
point $\widetilde{\Phi}^{*}$, the second term vanishes, namely
$R(\widetilde{\Phi}^{*})=0$, and
$\text{Hess}(\widetilde{\Phi}^{*})=\frac{2}{\widetilde{d}}(A^{*})^{\top}A^{*}$,
where $A^{*}:=A(\widetilde{\Phi}^{*})$ and $\widetilde{\Phi}^{*}$ are the
reduced phase factors corresponding to symmetric phase factors $\Phi^{*}$.

##  4 Global minimizer

According to the last section, the problem of finding _all_ global minima of
the optimization problem is converted to that of finding all admissible pairs
of $(P,Q)$ associated with the target polynomial $f$. In this section, we
provide a characterization of _all_ possible complementary polynomials
$(P_{\mathrm{Im}},Q)$.

###  4.1 Characterizing complementary polynomials

Before constructing all possible pairs of $(P(x),Q(x))$ satisfying Definition
2, we first introduce the following lemma, which views the problem from the
perspective of Laurent polynomials.

######  Lemma 14.

Given any real polynomial $f(x)$ of degree $d$ and with $\left\lVert
f\right\rVert_{\infty}<1$,
$\mathfrak{F}(z)=1-\left(f\left(\frac{z+z^{-1}}{2}\right)\right)^{2}$ has a
decomposition

$\mathfrak{F}(z)=\alpha\prod_{i=1}^{2d}(z-r_{i})(z^{-1}-r_{i}),$ (39)

such that $\alpha,r_{i}\in\mathbb{C}\backslash\\{0\\}$.

The decomposition is unique in the following sense: if
$\widetilde{\alpha}\in\mathbb{C}$ and $\widetilde{r}_{i}\in\mathbb{C}$,
$i=1,\cdots,n$ such that

$\mathfrak{F}(z)=\widetilde{\alpha}\prod_{i=1}^{n}(z-\widetilde{r}_{i})(z^{-1}-\widetilde{r}_{i}),$
(40)

then

$n=2d,\qquad\alpha\prod_{i=1}^{2d}r_{i}=\widetilde{\alpha}\prod_{i=1}^{2d}\widetilde{r}_{i},\qquad\frac{r_{i}+r_{i}^{-1}}{2}=\frac{\widetilde{r}_{i}+\widetilde{r}_{i}^{-1}}{2},$
(41)

up to permutation of the roots.

###### Proof.

According to the fundamental theorem of algebra, $1-f(x)^{2}$ has a unique
decomposition

$1-f(x)^{2}=c\prod_{i=1}^{2d}(s_{i}-x),\quad s_{i}\in\mathbb{C}.$

Since $\left\lVert f\right\rVert_{\infty}<1$, $s_{i}\neq 0$. The equation
$s=\frac{z+z^{-1}}{2}$ always has 2 nonzero solutions, say $z_{1},z_{2}$. In
particular, $z_{1}=z_{2}^{-1}$. For any $s_{i}$, pick one solution $r_{i}$,
use equality

$r_{i}(s_{i}-x)=\frac{1}{2}(z-r_{i})(z^{-1}-r_{i}),$

and one gets

$\mathfrak{F}(z)=c\prod_{i=1}^{2d}\frac{1}{2r_{i}}(z-r_{i})(z^{-1}-r_{i})=\alpha\prod_{i=1}^{2d}(z-r_{i})(z^{-1}-r_{i}),$

where $\alpha=c\prod_{i=1}^{2d}\frac{1}{2r_{i}}$. It provides a decomposition
of $\mathfrak{F}(z)$ satisfying the requirement.

On the other hand, given any decomposition

$\mathfrak{F}(z)=\widetilde{\alpha}\prod_{i=1}^{n}(z-\widetilde{r}_{i})(z^{-1}-\widetilde{r}_{i}),$

then one has

$\mathfrak{F}(z)=\widetilde{\alpha}\prod_{i=1}^{n}(2\widetilde{r}_{i})\left(\frac{\widetilde{r}_{i}+\widetilde{r}_{i}^{-1}}{2}-\frac{z+z^{-1}}{2}\right).$

It provides a decomposition of $1-f(x)^{2}$, i.e.

$1-f(x)^{2}=\widetilde{\alpha}\left(\prod_{i=1}^{n}2\widetilde{r}_{i}\right)\prod_{i=1}^{n}\left(\frac{\widetilde{r}_{i}+\widetilde{r}_{i}^{-1}}{2}-x\right).$

By the uniqueness decomposition of $1-f(x)^{2}$ in $\mathbb{C}[x]$, we obtain
Eq. 41 up to permutation of the roots. ∎

In particular, we are interested in the following decomposition of
$\mathfrak{F}(z)$, which is referred to as an _admissible decomposition_ ,

  * •

$\mathfrak{F}(z)=\alpha\prod_{i=1}^{2d}(z-r_{i})(z^{-1}-r_{i})$,
$\alpha>0,r_{i}\in\mathbb{C}\backslash\\{0\\}$.

  * •

Multiset $\mathcal{D}:=\\{r_{i}:1\leq i\leq 2d\\}$ is closed under complex
conjugation and additive inverse.

Firstly, $\mathcal{D}$ contains half of the roots of $\mathfrak{F}(z)$ with
multiplicity. (Hence $\mathcal{D}$ is a multiset rather than set.) Secondly,
note that $\alpha>0$ naturally comes from that
$\alpha=\frac{1-f(1)^{2}}{\prod_{r\in\mathcal{D}}(1-r)^{2}}>0$ for any
$\mathcal{D}$ closed under complex conjugation. Another observation is that
the roots of $\mathfrak{F}(z)$ cannot be any number of absolute value 1.
Otherwise, $1-f\left(\frac{r_{i}+r_{i}^{-1}}{2}\right)^{2}=0$,
i.e.$\left\lVert f\right\rVert_{\infty}=1$, which contradicts the constraint
on the maxnorm of $f$.

To be more precise, we rephrase Theorem 4 as follows and present its proof.

######  Theorem 15.

The pair of polynomials $\left(P(x),Q(x)\right)$ satisfy Definition 2 if only
if there exists an admissible decomposition of $\mathfrak{F}(z)$,

$\mathfrak{F}(z)=\alpha\prod_{i=1}^{2d}(z-r_{i})(z^{-1}-r_{i}),$ (42)

such that

$\begin{split}P_{\mathrm{Im}}\left(x\right)&=\sqrt{\alpha}\frac{e\left(x+\mathrm{i}\sqrt{1-x^{2}}\right)+e\left(x-\mathrm{i}\sqrt{1-x^{2}}\right)}{2},\\\
\quad
Q\left(x\right)&=\sqrt{\alpha}\frac{e\left(x+\mathrm{i}\sqrt{1-x^{2}}\right)-e\left(x-\mathrm{i}\sqrt{1-x^{2}}\right)}{2\mathrm{i}\sqrt{1-x^{2}}},\end{split}$
(43)

where $e(z):=z^{-d}\prod_{i=1}^{2d}(z-r_{i})$.

###### Proof.

$\Longrightarrow$: Denote
$p_{\mathrm{Im}}(z):=P_{\mathrm{Im}}\left(\frac{z+z^{-1}}{2}\right)$ and
$q(z):=Q\left(\frac{z+z^{-1}}{2}\right)$. Then the normalization condition
gives

$\begin{split}\mathfrak{F}(z)&=p_{\mathrm{Im}}(z)^{2}-q(z)^{2}\left(\frac{z-z^{-1}}{2}\right)^{2}\\\
&=\left(p_{\mathrm{Im}}(z)+\frac{z-z^{-1}}{2}q(z)\right)\left(p_{\mathrm{Im}}(z)-\frac{z-z^{-1}}{2}q(z)\right).\end{split}$

Define $H(z):=p_{\mathrm{Im}}(z)+\frac{z-z^{-1}}{2}q(z)$. Since
$p_{\mathrm{Im}}(z)=p_{\mathrm{Im}}(z^{-1})$ and $q(z)=q(z^{-1})$, one gets

$\mathfrak{F}(z)=H(z)H(z^{-1}).$

Notice that $z^{d}H(z)$ is a real polynomial of degree no more than $2d$. Thus
$z^{d}H(z)$ can be represented as $c\prod_{r\in\mathcal{D}}(z-r)$, where
$c\in\mathbb{C}$ and multiset $\mathcal{D}$ contains all the roots of
$z^{d}H(z)$ with multiplicity. Since $H(z)$ is real, $\mathcal{D}$ is closed
under complex conjugation and $c\in\mathbb{R}$.

With $H(z)=cz^{-d}\prod_{r\in\mathcal{D}}(z-r)$, one gets

$\mathfrak{F}(z)=H(z)H(z^{-1})=c^{2}\prod_{r\in\mathcal{D}}(z-r)(z^{-1}-r).$

By Lemma 14, the cardinality of $\mathcal{D}$ is $2d$ and thus $z^{d}H(z)$ is
of degree $2d$. Moreover, $r\in\mathbb{C}\backslash\\{0\\}$, $\forall
r\in\mathcal{D}$.

Besides, we can check that the parity of $z^{d}H(z)$ is even by examining the
parity of $p_{\mathrm{Im}}(z)$ and $q(z)$. Thus $\mathcal{D}$ is closed under
additive inverse and Eq. 43 holds as well.

$\Longleftarrow$: Denote
$p_{\mathrm{Im}}(z):=\sqrt{\alpha}\frac{e(z)+e(z^{-1})}{2}$ and
$q(z):=\sqrt{\alpha}\frac{e(z)-e(z^{-1})}{z-z^{-1}}$. Then

$\begin{split}\mathfrak{F}(z)&=\alpha\prod_{i=1}^{2d}(z-r_{i})(z^{-1}-r_{i})=\alpha
e(z)e(z^{-1})\\\
&=\left(\sqrt{\alpha}\frac{e(z)+e(z^{-1})}{2}\right)^{2}-\left(\sqrt{\alpha}\frac{e(z)-e(z^{-1})}{2}\right)^{2}\\\
&=\left(\sqrt{\alpha}\frac{e(z)+e(z^{-1})}{2}\right)^{2}+\left(1-\left(\frac{z+z^{-1}}{2}\right)^{2}\right)\left(\sqrt{\alpha}\frac{e(z)-e(z^{-1})}{z-z^{-1}}\right)^{2}\\\
&=p_{\mathrm{Im}}(z)^{2}+\left(1-\left(\frac{z+z^{-1}}{2}\right)^{2}\right)q(z)^{2}.\end{split}$

Since $\mathcal{D}$ is closed under complex conjugation,
$\prod_{i=1}^{2d}(z-r_{i})$ is a real polynomial and
$e(z)=z^{d}+\sum_{i=1-d}^{d-1}e_{i}z^{i}+\left(\prod_{i=1}^{2d}r_{i}\right)z^{-d}$
for some $e_{i}\in\mathbb{R}$. Notice that for $k>0$

$\begin{split}&\left(x+\mathrm{i}\sqrt{1-x^{2}}\right)^{k}+\left(x+\mathrm{i}\sqrt{1-x^{2}}\right)^{-k}=2T_{k}(x),\\\
&\left(x+\mathrm{i}\sqrt{1-x^{2}}\right)^{k}-\left(x+\mathrm{i}\sqrt{1-x^{2}}\right)^{-k}=2\mathrm{i}U_{k-1}(x)\sqrt{1-x^{2}}.\\\
\end{split}$

We can verify that $p_{\mathrm{Im}}(x+\mathrm{i}\sqrt{1-x^{2}})$ and
$q(x+\mathrm{i}\sqrt{1-x^{2}})$ are real polynomial of $x$. Hence, the
definition

$P_{\mathrm{Im}}(x):=p_{\mathrm{Im}}(x+\mathrm{i}\sqrt{1-x^{2}}),\quad\quad
Q(x):=q(x+\mathrm{i}\sqrt{1-x^{2}}).$

is well defined and $\left(P(x)=f(x)+\mathrm{i}P_{\mathrm{Im}}(x),Q(x)\right)$
satisfies the normalization condition,

$1-f(x)^{2}=P_{\mathrm{Im}}(x)^{2}+(1-x^{2})Q(x)^{2}.$

Since $\mathcal{D}$ is closed under additive inverse,

$e(-z)=(-z)^{-d}\prod_{i=1}^{2d}(-z-r_{i})=(-1)^{d}z^{-d}\prod_{i=1}^{2d}(z-r_{i})=(-1)^{d}e(z).$

It leads to that $P_{\mathrm{Im}}(x)$ has parity $(d\bmod 2)$ and $Q(x)$ has
parity $(d-1\bmod 2)$.

Another important observation is that
$\prod_{r\in\widetilde{\mathcal{D}}}r<0$. It comes from

$\mathfrak{F}(z)=\alpha\prod_{i=1}^{2d}(z-r_{i})(z^{-1}-r_{i})=\alpha\left(\prod_{i=1}^{2d}r_{i}\right)z^{-2d}\prod_{i=1}^{2d}(z-r_{i})(z-r_{i}^{-1})$
(44)

and that the coefficient of $z^{2d}$ in $\mathfrak{F}(z)$ is negative.
Therefore,

$\alpha\left(\prod_{i=1}^{2d}r_{i}\right)>0\Longrightarrow\left(\prod_{i=1}^{2d}r_{i}\right)<0.$
(45)

Since

$\begin{split}\frac{e(z)-e(z^{-1})}{z-z^{-1}}&=\frac{1}{z-z^{-1}}\left[\left(1-\prod_{i=1}^{2d}r_{i}\right)(z^{d}-z^{-d})+\sum_{i=1-d}^{d-1}e_{i}(z^{i}-z^{-i})\right]\\\
&=\left(1-\prod_{i=1}^{2d}r_{i}\right)\left(z^{d-1}+z^{-d+1}\right)+l.o..\end{split}$

Here, we use $l.o.$ to denote any linear combination of $z^{j}$ with
coefficients in $\mathbb{C}$ and $\left\lvert j\right\rvert<d-1$. By comparing
with the leading coefficient of the monomial expansion for $Q(x)$,
$\mathbf{q}_{d-1}$, one gets

$\sqrt{\alpha}\left(1-\prod_{i=1}^{2d}r_{i}\right)=2^{1-d}\mathbf{q}_{d-1}.$
(46)

In this way, $\text{deg}(Q)=d-1$, $\text{deg}(P)=d$ and the leading
coefficient of $Q(x)$ is positive. ∎

######  Remark 16 (Constructing complementary polynomials).

Theorem 15 indicates that there exists a bijection between the admissible
decomposition of $\mathfrak{F}(z)$ and the admissible pair $(P,Q)$ satisfying
Definition 2. Furthermore, it provides a constructive way to find all
admissible complementary polynomials:

  1. 1.

Choose a multiset $\mathcal{D}$ such that

    1. (a)

the joint union of $\mathcal{D}$ and its reciprocal contains all roots of
$\mathfrak{F}(z)$ with multiplicity.

    2. (b)

$\mathcal{D}$ is closed under complex conjugation and additive inverse.

  2. 2.

Construct $P_{\mathrm{Im}}$ and $Q$ according to Eq. 43.

###  4.2 Maximal solution

Among all possible choices of $\mathcal{D}$, we are interested in choosing
$\mathcal{D}$ to containing all roots of $\mathfrak{F}(z)$ within the unit
disc and outside the unit disc. As will be shown below, the associated
complementary polynomials satisfy many desirable properties. Let
$\left(P_{\mathrm{Im}},Q\right)$, $\left(\bar{P}_{\mathrm{Im}}.\bar{Q}\right)$
denote the corresponding pair of complementary polynomials respectively.

######  Lemma 17 (Maximal magnitude of leading coefficients).

The leading coefficients of the expansion of $Q$ and $\bar{Q}$ in the monomial
basis have the largest magnitude among all possible complementary polynomials.

###### Proof.

Given any multiset $\mathcal{D}$ satisfying the condition mentioned in Remark
16, note that $\mathfrak{F}(z)$ can be represented as

$\mathfrak{F}(z)=\beta z^{-2d}\prod_{r\in\mathcal{D}}(z-r)(z-r^{-1})$

for some $\beta\in\mathbb{C}$. By checking the leading coefficient,
$\beta=-\frac{\mathbf{f}^{2}_{d}}{4^{d}}$ is a quantity fully determined by
$f$.

Recall the proof of Theorem 15, we know that
$\alpha\prod_{r\in\mathcal{D}}r=\beta$ by Eq. 44. Together with Eq. 46, we
obtain

$\mathbf{q}_{d-1}^{2}=4^{d-1}\left\lvert\beta\right\rvert\left(K+\frac{1}{K}+2\right),$
(47)

where $K:=-\prod_{r\in\mathcal{D}}r>0$ by Eq. 45.

When $\mathcal{D}$ contains all the roots within the unit disc or outside the
unit disc, $\frac{1}{K}+K$ is the largest among all the possibilities, that
is, the leading coefficient of $Q(x)$ and $\bar{Q}(x)$ has the largest
magnitude. ∎

The following lemma shows that $\left(P_{\mathrm{Im}},Q\right)$ and
$\left(\bar{P}_{\mathrm{Im}},\bar{Q}\right)$ are directly related to each
other.

######  Lemma 18.

$\left(P_{\mathrm{Im}},Q\right)$ and
$\left(\bar{P}_{\mathrm{Im}}.\bar{Q}\right)$ have the following relation:

$P_{\mathrm{Im}}(x)=-\bar{P}_{\mathrm{Im}}(x),\quad\quad Q(x)=\bar{Q}(x).$
(48)

###### Proof.

Denote $\mathcal{D}:=\left\\{r\text{ is a root of }\mathfrak{F}(z):\left\lvert
r\right\rvert<1\right\\}$, and

$e(z):=z^{-d}\prod_{r\in\mathcal{D}}(z-r),\quad\bar{e}(z):=z^{-d}\prod_{r\in\mathcal{D}}(z-r^{-1}).$

Then for some $\alpha,\bar{\alpha}>0$,

$\mathfrak{F}(z)=\alpha e(z)e(z^{-1})=\bar{\alpha}\bar{e}(z)\bar{e}(z^{-1}).$
(49)

Notice that

$\begin{split}\bar{e}(z)&=z^{-d}\prod_{r\in\mathcal{D}}\left(r^{-1}(rz-1)\right)=z^{d}\prod_{r\in\mathcal{D}}\left(r^{-1}(r-z^{-1})\right)\\\
&=\left(\prod_{r\in\mathcal{D}}r^{-1}\right)z^{-d}\prod_{r\in\mathcal{D}}(z^{-1}-r)=-\frac{1}{K}e(z^{-1}).\end{split}$
(50)

where $K=-\prod_{r\in\mathcal{D}}r$. Similarly,
$\bar{e}(z^{-1})=-\frac{1}{K}e(z)$. Thus, $\bar{\alpha}=K^{2}\alpha$. By
Theorem 15, we have

$\begin{split}\bar{P}_{\mathrm{Im}}\left(x\right)&=\sqrt{\bar{\alpha}}\frac{\bar{e}\left(x+\mathrm{i}\sqrt{1-x^{2}}\right)+\bar{e}\left(x-\mathrm{i}\sqrt{1-x^{2}}\right)}{2},\\\
&=-\sqrt{\alpha
K^{2}}\frac{1}{K}\frac{e\left(x-\mathrm{i}\sqrt{1-x^{2}}\right)+e\left(x+\mathrm{i}\sqrt{1-x^{2}}\right)}{2}=-P_{\mathrm{Im}}(x),\end{split}$
(51)

and

$\begin{split}\bar{Q}\left(x\right)&=\sqrt{\bar{\alpha}}\frac{\bar{e}\left(x+\mathrm{i}\sqrt{1-x^{2}}\right)-\bar{e}\left(x-\mathrm{i}\sqrt{1-x^{2}}\right)}{2\mathrm{i}\sqrt{1-x^{2}}}\\\
&=-\sqrt{\alpha
K^{2}}\frac{1}{K}\frac{e\left(x-\mathrm{i}\sqrt{1-x^{2}}\right)-e\left(x+\mathrm{i}\sqrt{1-x^{2}}\right)}{2\mathrm{i}\sqrt{1-x^{2}}}=Q(x).\end{split}$
(52)

This proves the lemma. ∎

Due to this relation, we may only focus on the admissible pair
$\left(P(x),Q(x)\right)$ obtained by choosing $\mathcal{D}$ to contains all
the roots of $\mathfrak{F}(z)$ within the unit disc. Due to the property in
Lemma 17, the corresponding phase factors in $D_{d}$ will be referred to as
the maximal solution.

######  Remark 19 (Maximal solution at $f=0$).

Lemma 18 agrees with the observation that two particular admissible pairs
associated with $f=0$ are $\left(\mathrm{i}T_{d}(x),U_{d-1}(x)\right)$ and
$\left(-\mathrm{i}T_{d}(x),U_{d-1}(x)\right)$, which are obtained by

$\Phi^{0}:=\left(\frac{\pi}{4},0,\cdots,0,\frac{\pi}{4}\right)\quad\bar{\Phi}^{0}:=\left(-\frac{\pi}{4},0,\cdots,0,-\frac{\pi}{4}\right)$
(53)

respectively. In the following discussion, we will show that the maximal
solution $\left(P,Q\right)$ will converge to $\Phi^{0}$ as $\left\lVert
f\right\rVert_{\infty}\to 0$. The same argument is suitable for
$\left(\bar{P}_{\mathrm{Im}}.\bar{Q}\right)$ with $\Phi^{0}$ replaced by
$\bar{\Phi}^{0}$.

###  4.3 Relation to other constructive methods

There are two other methods for constructing complementary polynomials, which
will be referred to as GSLW method [6] and the Haah method [8], respectively.
The two methods are equivalent. The GSLW method solves the problem without
using Laurent polynomials. The Haah method first points out that the
decomposition can be viewed more naturally in terms of Laurent polynomials.

The Haah method considers

$1-P_{\mathrm{Re}}\left(\frac{z+z^{-1}}{2}\right)^{2}+\left(\frac{z-z^{-1}}{2}\right)^{2}Q_{\mathrm{Im}}\left(\frac{z+z^{-1}}{2}\right)^{2}$
(54)

and the corresponding decomposition:

$1-P_{\mathrm{Re}}\left(\frac{z+z^{-1}}{2}\right)^{2}+\left(\frac{z-z^{-1}}{2}\right)^{2}Q_{\mathrm{Im}}\left(\frac{z+z^{-1}}{2}\right)^{2}=\alpha\prod_{r\in\mathcal{D}}(z-r)(z^{-1}-r)$
(55)

where $\mathcal{D}$ contains all the roots within the unit disc.

In our construction, we neglect term $Q_{\mathrm{Im}}$ due to the symmetry
constraint. However, one can check that our method is also suitable for the
case that $Q_{\mathrm{Im}}\neq 0$ by simply replacing $\mathfrak{F}(z)$ by Eq.
54. By the same method as that in Theorem 15, it can shown that any admissible
complementary polynomial $(P_{\mathrm{Re}},Q_{\mathrm{Re}})$ must take the
form of Eq. 43, where $Q$ is replaced by $Q_{\mathrm{Re}}$. Hence our method
can be generalized to find all admissible complementary polynomials, by going
through all possible $e(z)=z^{-d}\prod_{r\in\mathcal{D}}(z-r)$. In particular,
the Haah method considers the special case that $\mathcal{D}$ contains all the
root of $\mathfrak{F}(z)$ within the unit disc, and therefore does not
generate all admissible complementary polynomials.

We would like to clarify that Ref [8] states that we can construct
$P_{\mathrm{Im}}$ and $Q$ as long as the joint union of a conjugate-closed
multiset $\mathcal{D}$ and its reciprocal $\\{\frac{1}{r}:r\in\mathcal{D}\\}$
is the list of all the roots of Eq. 54. However, according to Theorem 15,
another necessary condition for suitable complementary polynomials is that
$\mathcal{D}$ is closed under additive inverse. Without it, the parity
condition may not hold. This can be verified by the proof of Theorem 15.

##  5 Local convergence

We prove Theorem 5 in this section. Throughout this section, polynomials are
expanded in the Chebyshev basis without otherwise noted. This is more
convenient than the expansion in the monomial basis.

According to the definition of the Chebyshev norms in Section 2.1 and the
convention of the coefficients in Section 2.2, we specifically list the
following relations which exploit the weighted orthogonality and will be
frequently used in the proof,

$\begin{split}\left\lVert
f\right\rVert_{T}^{2}=f^{2}_{0}+\frac{1}{2}\sum_{k=1}^{d}f^{2}_{k},\quad\quad\left\lVert
P_{\mathrm{Im}}\right\rVert_{T}^{2}=s_{0}^{2}+\frac{1}{2}\sum_{k=1}^{d}s_{k}^{2},\quad\quad\left\lVert
Q\right\rVert_{U}^{2}=\frac{1}{2}\sum_{k=0}^{d-1}{q_{k}}^{2}.\end{split}$ (56)

By integrating the point-wise normalization condition
$f(x)^{2}+P_{\mathrm{Im}}(x)^{2}+(1-x^{2})Q(x)^{2}=1$ with the weight function
$\frac{1}{\sqrt{1-x^{2}}}$, we obtain the following Parseval’s equality:

$f^{2}_{0}+s_{0}^{2}+\frac{1}{2}\sum_{k=1}^{d}\left(f^{2}_{k}+s_{k}^{2}\right)+\frac{1}{2}\sum_{k=0}^{d-1}{q_{k}}^{2}=1.$
(57)

###  5.1 Sketch of the proof

We first give a sketch of the proof as a road map. For simplicity, we assume
$d$ is odd in the sketch, while the proof works for any $d$.

  1. 1.

We first show that each element of $\widetilde{\Phi}^{*}-\widetilde{\Phi}^{0}$
is bounded by $\frac{\pi}{4}$ by properly restricting $\left\lVert
f\right\rVert_{\infty}$.

  2. 2.

We further upper bound the distance
$\left\lVert\widetilde{\Phi}^{*}-\widetilde{\Phi}^{0}\right\rVert_{2}$ by
$\left\lvert\sin\left(2\phi_{0}-\frac{\pi}{2}\right)\right\rvert^{2}+\sum_{\ell=1}^{\widetilde{d}-1}\left\lvert\sin(2\phi_{\ell})\right\rvert^{2}$.
It follows the inequality $\left\lvert
x\right\rvert\leq\frac{\pi}{2}\left\lvert\sin(x)\right\rvert$ which holds for
any $\left\lvert x\right\rvert\leq\frac{\pi}{2}$. The bound in the first step
ensures that the inequality can be invoked.

  3. 3.

We upper bound the phase factors by the leading Chebyshev coefficients in the
recursion sequence of polynomials. The upper bound turns out to be

$\begin{split}&\left\lvert\sin\left(2\phi_{0}-\frac{\pi}{2}\right)\right\rvert^{2}\leq\frac{8}{3}\left\lVert
f\right\rVert_{\infty}^{2},\\\
&\left\lvert\sin(2\phi_{\ell})\right\rvert^{2}\leq\frac{8}{3}\left(\left({q^{(\ell)}_{d-1-2\ell}}\right)^{2}-\left({q^{(\ell-1)}_{d+1-2\ell}}\right)^{2}\right),\quad
1\leq\ell\leq\widetilde{d}-1.\end{split}$

  4. 4.

Combining the previous results, we finally get
$\left\lVert\widetilde{\Phi}^{*}-\widetilde{\Phi}^{0}\right\rVert_{2}^{2}\leq\frac{\pi^{2}}{6}\left(\left\lVert
f\right\rVert_{\infty}^{2}+1-q_{d-1}^{2}\right)$. The upper bound can be
further simplified to be $\frac{\pi^{2}}{3}\left\lVert
f\right\rVert_{\infty}^{2}$ by exploiting the estimate to $q_{d-1}$ (see Lemma
21). Then the statement in Theorem 5 follows.

###  5.2 Recursion relation in terms of Chebyshev coefficients

In this subsection, we will convert the recursion relation in terms of
monomial coefficients to that in terms of Chebyshev coefficients. Furthermore,
we will derive some useful properties which will be used as intermediate steps
in the proof of Theorem 5.

Let $(P,Q)\in\mathbb{C}_{d}[x]\times\mathbb{R}_{d-1}[x]$ be a pair of
polynomials satisfying the conditions in Theorem 1. The symmetric reduction
procedures yield a sequence of polynomials
$(P^{(\ell)},Q^{(\ell)})\in\mathbb{C}_{d-2\ell}[x]\times\mathbb{R}_{d-2\ell-1}[x]$
satisfying the conditions and an associated sequence of phase factors
$\phi_{\ell}$ so that
$e^{2\mathrm{i}\phi_{\ell}}=\frac{\mathbf{p}^{(\ell)}_{d-2\ell}}{\mathbf{q}^{(\ell)}_{d-1-2\ell}}$.
In particular, $(P^{(0)},Q^{(0)})=(P,Q)$. An important recursion relation
generating the coefficients of new polynomials in the monomial basis is given
in Eq. 32. For convenience, we restate it as follows: if $d\geq 3$, then

$\begin{split}&\mathbf{p}_{d-2}^{(1)}=\frac{\mathbf{q}_{d-1}}{4}-\frac{1}{\mathbf{q}_{d-1}}\left\lvert\mathrm{Im}[\mathbf{p}_{d-2}e^{-2\mathrm{i}\phi_{0}}]\right\rvert^{2}+\mathrm{i}\mathrm{Im}[\mathbf{p}_{d-2}e^{-2\mathrm{i}\phi_{0}}],\\\
&\mathbf{q}_{d-3}^{(1)}=\frac{\mathbf{q}_{d-1}}{4}+\frac{1}{\mathbf{q}_{d-1}}\left\lvert\mathrm{Im}[\mathbf{p}_{d-2}e^{-2\mathrm{i}\phi_{0}}]\right\rvert^{2}.\end{split}$

Here, recall that the boldface symbols refer to the coefficients in the
monomial basis. We will first rewrite this relation using the Chebyshev basis.
As a remark, it suffices to consider the case when $\ell=0$, and the general
case can be recovered by substituting
$(P,Q,P^{(1)},Q^{(1)},\phi_{0})\leftarrow(P^{(\ell)},Q^{(\ell)},P^{(\ell+1)},Q^{(\ell+1)},\phi_{\ell})$
in the arguments as long as $\deg\left(P^{(\ell)}\right)\geq 3$.

For simplicity, let us denote the Chebyshev polynomial in the monomial basis
as $T_{d}(x)=2^{d-1}x^{d}+t_{d}x^{d-2}+l.o.$ where we use $l.o.$ to denote any
polynomial of degree no more than $d-3$ and $t_{d}$ is a real number whose
exact value is irrelevant in the proof.

Expanding the polynomials in mononial basis and Chebyshev basis, we have

$\begin{split}&P(x)e^{-2\mathrm{i}\phi_{0}}=\sum_{k=0}^{d}\mathbf{p}_{k}e^{-2\mathrm{i}\phi_{0}}x^{k}=e^{-2\mathrm{i}\phi_{0}}\mathbf{p}_{d}x^{d}+e^{-2\mathrm{i}\phi_{0}}\mathbf{p}_{d-2}x^{d-2}+l.o.\\\
&=\sum_{k=0}^{d}p_{k}e^{-2\mathrm{i}\phi_{0}}T_{k}(x)=e^{-2\mathrm{i}\phi_{0}}2^{d-1}p_{d}x^{d}+\left(e^{-2\mathrm{i}\phi_{0}}p_{d}t_{d}+e^{-2\mathrm{i}\phi_{0}}p_{d-2}2^{d-3}\right)x^{d-2}+l.o..\end{split}$

Then, by comparing the coefficients of terms $x^{d}$ and $x^{d-2}$, we have
$\mathbf{p}_{d}=2^{d-1}p_{d}$ and

$\begin{split}e^{-2\mathrm{i}\phi_{0}}\mathbf{p}_{d-2}=e^{-2\mathrm{i}\phi_{0}}p_{d}t_{d}+e^{-2\mathrm{i}\phi_{0}}p_{d-2}2^{d-3}=q_{d-1}t_{d}+e^{-2\mathrm{i}\phi_{0}}p_{d-2}2^{d-3}.\end{split}$

In the last equality, we use

$e^{2\mathrm{i}\phi_{0}}=\frac{\mathbf{p}_{d}}{\mathbf{q}_{d-1}}=\frac{2^{d-1}p_{d}}{2^{d-1}q_{d-1}}=\frac{p_{d}}{q_{d-1}}\Rightarrow
e^{-2\mathrm{i}\phi_{0}}p_{d}=q_{d-1}\in\mathbb{R}.$

An important observation is

$\mathrm{Im}[\mathbf{p}_{d-2}e^{-2\mathrm{i}\phi_{0}}]=2^{d-3}\mathrm{Im}[p_{d-2}e^{-2\mathrm{i}\phi_{0}}],$
(58)

which follows $q_{d-1}t_{d}\in\mathbb{R}$.

Plugging the previous equations in terms of Chebyshev coefficients in Eq. 32,
we have the following recursion relation

$\begin{split}&p_{d-2}^{(1)}=q_{d-1}-\frac{1}{4q_{d-1}}\left\lvert\mathrm{Im}[p_{d-2}e^{-2\mathrm{i}\phi_{0}}]\right\rvert^{2}+\mathrm{i}\mathrm{Im}[p_{d-2}e^{-2\mathrm{i}\phi_{0}}],\\\
&q_{d-3}^{(1)}=q_{d-1}+\frac{1}{4q_{d-1}}\left\lvert\mathrm{Im}[p_{d-2}e^{-2\mathrm{i}\phi_{0}}]\right\rvert^{2}.\end{split}$
(59)

The second equality implies the following sign preserving property and the
monotonicity

$\text{sgn}(q_{d-1})=\text{sgn}(q_{d-3}^{(1)})\text{, and }\left\lvert
q_{d-1}\right\rvert\leq\left\lvert q_{d-3}^{(1)}\right\rvert.$

Regarding the next phase factor $\phi_{1}$, its defined equation is converted
similarly to Chebyshev basis

$e^{2\mathrm{i}\phi_{1}}=\frac{\mathbf{p}^{(1)}_{d-2}}{\mathbf{q}^{(1)}_{d-3}}=\frac{p^{(1)}_{d-2}}{q^{(1)}_{d-3}}.$

Furthermore, we can use the derived recursion relation in Eq. 59 to obtain

$\sin\left(2\phi_{1}\right)=\mathrm{Im}\left[\frac{p^{(1)}_{d-2}}{q^{(1)}_{d-3}}\right]=\frac{\mathrm{Im}[p_{d-2}e^{-2\mathrm{i}\phi_{0}}]}{q_{d-3}^{(1)}}.$

The derived results can be generalized to $\left(P^{(\ell)},Q^{(\ell)}\right)$
where $1\leq\ell\leq\hat{d}$ ($\hat{d}$ is defined by Eq. 35). We summarize
these results as the following theorem.

######  Theorem 20.

Given $\left(P(x),Q(x)\right)$ satisfying the assumption in Theorem 1 and
$\Phi=(\phi_{0},\phi_{1},\cdots,\phi_{1},\phi_{0})$ being a set of symmetric
phase factors in $\mathcal{D}_{d}$ such that Eq. 8 is satisfied, the following
statements hold.

  1. 1.

For any $1\leq\ell\leq\hat{d}$, the leading coefficients of $P^{(\ell)}$ and
$Q^{(\ell)}$ satisfy the recursion relation

$\begin{split}&p^{(\ell)}_{d-2\ell}=q^{(\ell-1)}_{d+1-2\ell}-\frac{1}{4q^{(\ell-1)}_{d+1-2\ell}}\left\lvert\mathrm{Im}[p^{(\ell-1)}_{d-2\ell}e^{-2\mathrm{i}\phi_{\ell-1}}]\right\rvert^{2}+\mathrm{i}\mathrm{Im}[p^{(\ell)}_{d-2\ell}e^{-2\mathrm{i}\phi_{\ell-1}}],\\\
&q^{(\ell)}_{d-1-2\ell}=q^{(\ell-1)}_{d+1-2\ell}+\frac{1}{4q^{(\ell-1)}_{d+1-2\ell}}\left\lvert\mathrm{Im}[p^{(\ell-1)}_{d-2\ell}e^{-2\mathrm{i}\phi_{\ell-1}}]\right\rvert^{2}.\end{split}$
(60)

  2. 2.

For any $0\leq\ell\leq\hat{d}$,

$e^{2\mathrm{i}\phi_{\ell}}=\frac{p^{(\ell)}_{d-2\ell}}{q^{(\ell)}_{d-1-2\ell}},$
(61)

and if $1\leq\ell\leq\hat{d}$,

$\sin\left(2\phi_{\ell}\right)=\frac{\mathrm{Im}[p^{(\ell-1)}_{d-2\ell}e^{-2\mathrm{i}\phi_{\ell-1}}]}{q_{d-1-2\ell}^{(\ell)}}.$
(62)

  3. 3.

If $d$ is even,
$e^{\mathrm{i}\phi_{\widetilde{d}-1}}=p_{0}^{(\widetilde{d}-1)}$ and

$\cos(\phi_{\widetilde{d}-1})=q_{1}^{(\widetilde{d}-2)},\
\sin\left(\phi_{\widetilde{d}-1}\right)=-\mathrm{i}p^{(\widetilde{d}-2)}_{0}e^{-2\mathrm{i}\phi_{\widetilde{d}-2}}.$
(63)

  4. 4.

(The monotonicity of the leading coefficient of $Q^{(\ell)}$ in the Chebyshev
basis)

If $q^{(0)}_{d-1}>0$, then

$0<q^{(0)}_{d-1}\leq q^{(1)}_{d-3}\leq\cdots\leq
q^{(\ell)}_{d-1-2\ell}\leq\cdots\leq q^{(\hat{d})}_{d-1-2\hat{d}}\leq 1.$ (64)

###### Proof.

The statements 1, 2 and 4 follow the previous discussion. Hence, we only need
to verify Eq. 63.

Note that when $\ell=\widetilde{d}-2$, we have
$\deg\left(P^{(\widetilde{d}-2)}\right)=2$ and

$\begin{split}&\braket{0}{e^{\mathrm{i}\phi_{\widetilde{d}-2}Z}W(x)e^{\mathrm{i}\phi_{\widetilde{d}-1}}W(x)e^{\mathrm{i}\phi_{\widetilde{d}-2}Z}}{0}=e^{2\mathrm{i}\phi_{\widetilde{d}-2}}\left(2\cos\left(\phi_{\widetilde{d}-1}\right)x^{2}-e^{-\mathrm{i}\phi_{\widetilde{d}-1}}\right)\\\
&=P^{(\widetilde{d}-2)}(x)=p_{2}^{(\widetilde{d}-2)}T_{2}(x)+p_{0}^{(\widetilde{d}-2)}T_{0}(x)=2p_{2}^{(\widetilde{d}-2)}x^{2}-p_{2}^{(\widetilde{d}-2)}+p_{0}^{(\widetilde{d}-2)}.\end{split}$

The first equality is obtained by direct computation. Comparing the constant
terms, we have

$\begin{split}e^{\mathrm{i}\phi_{\widetilde{d}-1}}&=\left(p_{2}^{(\widetilde{d}-2)}e^{-2\mathrm{i}\phi_{\widetilde{d}-2}}\right)^{*}-\left(p_{0}^{(\widetilde{d}-2)}e^{-2\mathrm{i}\phi_{\widetilde{d}-2}}\right)^{*}\\\
&=q_{1}^{(\widetilde{d}-2)}-\left(p_{0}^{(\widetilde{d}-2)}e^{-2\mathrm{i}\phi_{\widetilde{d}-2}}\right)^{*}.\end{split}$

Taking the modulus on both sides, we get

$\begin{split}1&=\left(q_{1}^{(\widetilde{d}-2)}-\mathrm{Re}[p_{0}^{(\widetilde{d}-2)}e^{-2\mathrm{i}\phi_{\widetilde{d}-2}}]\right)^{2}+\left(\mathrm{Im}[p^{(\widetilde{d}-2)}_{0}e^{-2\mathrm{i}\phi_{\widetilde{d}-2}}]\right)^{2}\\\
&=\left(q_{1}^{(\widetilde{d}-2)}\right)^{2}+\left\lvert
p_{0}^{(\widetilde{d}-2)}\right\rvert^{2}-2q_{1}^{(\widetilde{d}-2)}\mathrm{Re}[p_{0}^{(\widetilde{d}-2)}e^{-2\mathrm{i}\phi_{\widetilde{d}-2}}].\end{split}$
(65)

Applying Parseval’s equality derived from $\left\lvert
P^{(\widetilde{d}-2)}(x)\right\rvert^{2}+(1-x^{2})\left\lvert
Q^{(\widetilde{d}-2)}(x)\right\rvert^{2}=1$, we get

$\left\lvert p_{0}^{(\widetilde{d}-2)}\right\rvert^{2}+\frac{1}{2}\left\lvert
p_{2}^{(\widetilde{d}-2)}\right\rvert^{2}+\frac{1}{2}\left(q_{1}^{(\widetilde{d}-2)}\right)^{2}=1.$

Note that $\left\lvert p_{2}^{(\widetilde{d}-2)}\right\rvert=\left\lvert
p_{2}^{(\widetilde{d}-2)}e^{-2\mathrm{i}\phi_{\widetilde{d}-2}}\right\rvert=\left\lvert
q_{1}^{(\widetilde{d}-2)}\right\rvert$, we have

$\left\lvert
p_{0}^{(\widetilde{d}-2)}\right\rvert^{2}+\left(q_{1}^{(\widetilde{d}-2)}\right)^{2}=1.$
(66)

Plugging the last equality into Eq. 65, we obtain
$\mathrm{Re}[p_{0}^{(\widetilde{d}-2)}e^{-2\mathrm{i}\phi_{\widetilde{d}-2}}]=0$
since $q_{1}^{(\widetilde{d}-2)}\neq 0$. Therefore,
$p^{(\widetilde{d}-2)}_{0}e^{-2\mathrm{i}\phi_{\widetilde{d}-2}}$ is pure
imaginary which yields Eq. 63. ∎

###  5.3 Bounding Chebyshev coefficients of the maximal solution

In this subsection, we assume that $(P,Q)$ is the admissible pair
corresponding to the maximal solution. We will derive bounds on the leading
coefficients of $Q$ and $P_{\mathrm{Im}}(x)$ in the Chebyshev basis, i.e.,
$q_{d-1}$ and $s_{d}$ respectively.

######  Lemma 21.

The leading coefficient of $Q$ in the Chebyshev basis satisfies the following
estimate: $1-\left\lVert f\right\rVert^{2}_{\infty}\leq q_{d-1}^{2}\leq
1-\left\lVert f\right\rVert^{2}_{T}+\frac{1}{2}f_{d}^{2}$.

###### Proof.

Expanding the polynomials in the monomial basis, the normalization condition
reduces to

$f(x)^{2}+P_{\mathrm{Im}}(x)^{2}+(1-x^{2})Q(x)^{2}=\left(\mathbf{f}_{d}^{2}+\mathbf{s}_{d}^{2}-\mathbf{q}_{d}^{2}\right)x^{2d}+l.o.=1,$

where $l.o.$ denotes any polynomial of degree less than $2d$. It implies that
$\mathbf{f}_{d}^{2}+{\mathbf{s}_{d}}^{2}=\mathbf{q}_{d-1}^{2}$ which is
further $f_{d}^{2}+{s_{d}}^{2}=q_{d-1}^{2}$. Hence, we have

$\begin{split}q_{d-1}^{2}&=\frac{1}{2}\left(f^{2}_{d}+s_{d}^{2}+q_{d-1}^{2}\right)\\\
&=1-\left(f_{0}^{2}+s_{0}^{2}+\frac{1}{2}\sum_{k=1}^{d-1}\left(f^{2}_{k}+s_{k}^{2}\right)+\frac{1}{2}\sum_{k=0}^{d-2}q_{k}^{2}\right)\\\
&\leq
1-\left(f_{0}^{2}+\frac{1}{2}\sum_{k=1}^{d-1}f^{2}_{k}\right)=1-\left\lVert
f\right\rVert^{2}_{T}+\frac{1}{2}f_{d}^{2}.\\\ \end{split}$

In the second line, we exploit Eq. 57, the Parseval’s equality derived from
the normalization condition.

To lower bound $q_{d-1}$, we first remark that Eq. 47 can be rewritten as

$q_{d-1}=\left\lvert\beta\right\rvert\left(K+\frac{1}{K}+2\right),$ (67)

where $K:=-\prod_{r\in\mathcal{D}}r>0$, and $\mathcal{D}$ contains all roots
of $\mathfrak{F}(z)$ within the unit disc for the maximal solution.

Given a polynomial fully factored into
$g(z)=\gamma\prod_{i=1}^{n}(z-a_{i})\in\mathbb{C}[z]$, its _Mahler measure_
[15] is defined as

$M(g):=\left\lvert\gamma\right\rvert\prod_{i=1}^{n}\max\\{\left\lvert
a_{i}\right\rvert,1\\}.$

This quantity is equal to the geometric mean of $\left\lvert g(z)\right\rvert$
on the unit circle by using Jensen’s formula in [10]

$M(g)=\exp\left(\frac{1}{2\pi}\int_{0}^{2\pi}\ln{\left\lvert
g(e^{\mathrm{i}\theta})\right\rvert}\mathrm{d}\theta\right).$ (68)

Remarkably, note that $\left\lvert\beta\right\rvert\frac{1}{K}$ is the Mahler
measure of a polynomial $h(z)$ defined as

$h(z):=\beta\prod_{r\in\mathcal{D}}(z-r)(z-r^{-1})=z^{2d}\mathfrak{F}(z).$
(69)

Therefore, one has

$\begin{split}\left\lvert\beta\right\rvert\frac{1}{K}&=M(h)=\exp\left(\frac{1}{2\pi}\int_{0}^{2\pi}\ln{\left\lvert
h(e^{\mathrm{i}\theta})\right\rvert}\mathrm{d}\theta\right)\\\
&=\exp\left(\frac{1}{2\pi}\int_{0}^{2\pi}\ln{\left(1-f(\cos{\theta})^{2}\right)}\mathrm{d}\theta\right)\\\
&\geq\exp\left(\frac{1}{2\pi}\int_{0}^{2\pi}\ln{\left(1-\left\lVert
f\right\rVert_{\infty}^{2}\right)}\mathrm{d}\theta\right)\\\ &=1-\left\lVert
f\right\rVert_{\infty}^{2}.\end{split}$ (70)

To conclude,

$q_{d-1}^{2}=\left\lvert\beta\right\rvert\left(K+\frac{1}{K}+2\right)\geq\left\lvert\beta\right\rvert\frac{1}{K}\geq
1-\left\lVert f\right\rVert_{\infty}^{2},$ (71)

which completes the proof. ∎

######  Remark 22 (Convergence of the maximal solution towards $\Phi^{0}$).

The derived bounds on $q_{d-1}$ consequentially yields the bounds on $s_{d}$.
Using the equality $f_{d}^{2}+s_{d}^{2}=q_{d-1}^{2}$, we get the following
estimate

$1-\left\lVert f\right\rVert^{2}_{\infty}-f_{d}^{2}\leq s_{d}^{2}\leq
1-\left\lVert f\right\rVert^{2}_{T}-\frac{1}{2}f_{d}^{2}.$ (72)

Combining the derived results, the remaining Chebyshev coefficients of
$P_{\mathrm{Im}}$ and $Q$ can be upper bounded as follows

$\begin{split}&s_{0}^{2}+\frac{1}{2}\sum_{k=1}^{d-1}s_{k}^{2}+\frac{1}{2}\sum_{k=0}^{d-2}q_{k}^{2}\\\
(\text{Parseval's equality})\
&=1-f_{0}^{2}-\frac{1}{2}\sum_{k=1}^{d}f_{k}^{2}-\frac{1}{2}s_{d}^{2}-\frac{1}{2}q_{d-1}^{2}\\\
(f_{d}^{2}+s_{d}^{2}=q_{d-1}^{2})\
&=1-q_{d-1}^{2}-f_{0}^{2}-\frac{1}{2}\sum_{k=1}^{d-1}f_{k}^{2}\\\
(\text{\lx@cref{creftype~refnum}{lma:coef_q_bound}})\ &\leq\left\lVert
f\right\rVert_{\infty}^{2}-f_{0}^{2}-\frac{1}{2}\sum_{k=1}^{d-1}f_{k}^{2}\leq\left\lVert
f\right\rVert_{\infty}^{2}.\end{split}$ (73)

As a consequence of the upper bound, we have $s_{i},q_{j}\to 0\ \forall
i\in[d],j\in[d-1]$ as $\left\lVert f\right\rVert_{\infty}\to 0$. At the same
time, $f_{d}^{2}+s_{d}^{2}=q_{d-1}^{2}$ implies that $s_{d}=q_{d-1}$ in that
limit since $f_{d}\to 0$, and the normalization condition yields
$s_{d}=q_{d-1}=1$ when the leading coefficient of $Q$ is fixed to be positive.
To conclude, as $\left\lVert f\right\rVert_{\infty}\to 0$, the maximal
solution exhibits a limit $(P,Q)\to(\mathrm{i}T_{d},U_{d-1})$ which agrees
with the polynomials generated from the set of phase factors
$\Phi^{0}=(\frac{\pi}{4},0,\cdots,0,\frac{\pi}{4})$.

###  5.4 Proof of Theorem 5

In this subsection, we will prove Theorem 5 based on the previous results. We
start from a technical lemma strengthening the results in Theorem 20 under
additional assumptions.

######  Lemma 23.

When $\left\lVert f\right\rVert_{\infty}\leq\frac{1}{2}$, it holds that

  1. 1.

$\frac{\sqrt{3}}{2}\leq\sqrt{1-\left\lVert f\right\rVert_{\infty}^{2}}\leq
q^{(0)}_{d-1}\leq q^{(1)}_{d-3}\leq\cdots\leq
q^{(\ell)}_{d-1-2\ell}\leq\cdots\leq q^{(\hat{d})}_{d-1-2\hat{d}}\leq 1$.

  2. 2.

$\left\lvert\phi_{0}-\frac{\pi}{4}\right\rvert<\frac{\pi}{4}$ and
$\left\lvert\phi_{\ell}\right\rvert<\frac{\pi}{4},\quad\forall
1\leq\ell\leq\widetilde{d}-1$.

###### Proof.

The first statement directly follows Lemma 21.

To show the second statement, we first notice that according to Theorem 20,

$\cos\left(2\phi_{0}-\frac{\pi}{2}\right)=\sin(2\phi_{0})=\frac{\mathrm{Im}[p_{d}^{(0)}]}{q_{d-1}^{(0)}}.$

Following Eq. 43, we have

$\begin{split}P_{\mathrm{Im}}(x)&=\sqrt{\alpha}\frac{e\left(x+\mathrm{i}\sqrt{1-x^{2}}\right)+e\left(x-\mathrm{i}\sqrt{1-x^{2}}\right)}{2}\\\
&=\sqrt{\alpha}\frac{1+\prod_{r\in\mathcal{D}}r}{2}\left[\left(x+\mathrm{i}\sqrt{1-x^{2}}\right)^{d}+\left(x-\mathrm{i}\sqrt{1-x^{2}}\right)^{d}\right]+l.o.\\\
&=\sqrt{\alpha}\left(1+\prod_{r\in\mathcal{D}}r\right)T_{d}(x)+l.o..\end{split}$
(74)

Here $l.o.$ represents any polynomial of degree $<d$. To achieve the second
equality, we use the expansion of $e(z)=z^{-d}\prod_{r\in\mathcal{D}}(z-r)$ in
the monomial basis and discover that only terms $z^{d}$ and
$\left(\prod_{r\in\mathcal{D}}r\right)z^{-d}$ contribute to the first term on
the second line, if we replace $z,z^{-1}$ with
$x+\mathrm{i}\sqrt{1-x^{2}},x-\mathrm{i}\sqrt{1-x^{2}}$ respectively.

In this way, one gets

$\mathrm{Im}[p_{d}^{(0)}]=\sqrt{\alpha}\left(1+\prod_{r\in\mathcal{D}}r\right)>0,$
(75)

where $\mathcal{D}$ contains the roots of $\mathfrak{F}(z)$ in the unit disc.
Therefore, $\left\lvert\phi_{0}-\frac{\pi}{4}\right\rvert<\frac{\pi}{4}$
follows $\cos\left(2\phi_{0}-\frac{\pi}{2}\right)>0$.

When $1\leq\ell\leq\hat{d}$, according to Eqs. 60 and 61, we have

$\cos\left(2\phi_{\ell}\right)=\frac{1}{q^{(\ell)}_{d-1-2\ell}}\left(q^{(\ell-1)}_{d+1-2\ell}-\frac{1}{4q^{(\ell-1)}_{d+1-2\ell}}\left\lvert\mathrm{Im}[p^{(\ell-1)}_{d-2\ell}e^{-2\mathrm{i}\phi_{\ell-1}}]\right\rvert^{2}\right).$
(76)

Using Parseval’s equality, the normalization condition $\left\lvert
P^{(\ell-1)}\right\rvert^{2}+(1-x^{2})\left\lvert
Q^{(\ell-1)}\right\rvert^{2}=1$ yields

$\begin{split}1&=\left\lvert
p^{(\ell-1)}_{0}\right\rvert^{2}+\frac{1}{2}\sum_{j=1}^{d-2\ell+2}\left\lvert
p_{j}^{(\ell-1)}\right\rvert^{2}+\frac{1}{2}\sum_{j=0}^{d-2\ell+1}\left(q_{j}^{(\ell-1)}\right)^{2}\\\
&\geq\frac{1}{2}\left\lvert
p^{(\ell-1)}_{d-2\ell}\right\rvert^{2}+\frac{1}{2}\left\lvert
p^{(\ell-1)}_{d-2\ell+2}\right\rvert^{2}+\frac{1}{2}\left(q_{d-2\ell+1}^{(\ell-1)}\right)^{2}.\end{split}$

Furthermore, the cancellation of the leading order in the normalization
condition implies $\left\lvert
p^{(\ell-1)}_{d-2\ell+2}\right\rvert^{2}=\left(q_{d-2\ell+1}^{(\ell-1)}\right)^{2}$.
Thus

$\frac{1}{2}\left\lvert
p^{(\ell-1)}_{d-2\ell}\right\rvert^{2}+\left(q_{d-2\ell+1}^{(\ell-1)}\right)^{2}\leq
1.$

Using the first statement frequently, we have

$\begin{split}&\left\lvert\mathrm{Im}[p^{(\ell-1)}_{d-2\ell}e^{-2\mathrm{i}\phi_{\ell-1}}]\right\rvert^{2}\leq\left\lvert
p^{(\ell-1)}_{d-2\ell}\right\rvert^{2}\leq
2\left(1-\left(q_{d-2\ell+1}^{(\ell-1)}\right)^{2}\right)\\\ &\leq
2\left(1-\left(q_{d-1}^{(0)}\right)^{2}\right)\leq 2\left\lVert
f\right\rVert_{\infty}^{2}\leq 4\left\lVert
f\right\rVert^{2}_{\infty}<4\left(1-\left\lVert
f\right\rVert_{\infty}^{2}\right)\leq
4\left(q^{(\ell-1)}_{d+1-2\ell}\right)^{2}.\end{split}$ (77)

Plugging Eq. 77 into Eq. 76, we have $\cos\left(2\phi_{\ell}\right)>0$. Hence,
$\left\lvert\phi_{\ell}\right\rvert<\frac{\pi}{4}$.

If $d$ is even, following Theorem 20 statement 3, we have

$\cos(\phi_{\widetilde{d}-1})=q_{1}^{(\widetilde{d}-2)}=q_{1}^{(\hat{d})}\geq\frac{\sqrt{3}}{2}.$

Here, $\hat{d}=\widetilde{d}-2$ when $d$ is even and the first statement is
invoked. It implies that
$\left\lvert\phi_{\widetilde{d}-1}\right\rvert<\frac{\pi}{4}$. ∎

The second statement of the previous lemma ensures that the inequality
$\frac{\pi}{2}\left\lvert\sin(x)\right\rvert\geq\left\lvert x\right\rvert$ for
any $\left\lvert x\right\rvert\leq\frac{\pi}{2}$ can be invoked. Now we are
ready to prove Theorem 5.

###### Proof of Theorem 5.

We first estimate $\left\lvert\sin(2\phi_{0}-\frac{\pi}{2})\right\rvert^{2}$
as follows

$\left\lvert\sin(2\phi_{0}-\frac{\pi}{2})\right\rvert^{2}=\left\lvert\frac{\mathrm{Re}[p_{d}]}{q_{d-1}}\right\rvert^{2}=\frac{f_{d}^{2}}{q_{d-1}^{2}}\leq\frac{4}{3}f_{d}^{2}\leq\frac{8}{3}\left\lVert
f\right\rVert_{\infty}^{2},$ (78)

where we use the inequality $f_{d}^{2}\leq 2\left\lVert
f\right\rVert_{T}^{2}\leq 2\left\lVert f\right\rVert_{\infty}^{2}$.

For any $1\leq\ell\leq\hat{d}$, by Eq. 62, we have

$\begin{split}\left\lvert\sin(2\phi_{\ell})\right\rvert^{2}&=\left\lvert\frac{\mathrm{Im}[p^{(\ell-1)}_{d-2\ell}e^{-2\mathrm{i}\phi_{\ell-1}}]}{q_{d-1-2\ell}^{(\ell)}}\right\rvert^{2}\\\
\text{(\lx@cref{creftype~refnum}{lma:supplement}
1.})&\leq\frac{4}{3}\left\lvert\mathrm{Im}[p^{(\ell-1)}_{d-2\ell}e^{-2\mathrm{i}\phi_{\ell-1}}]\right\rvert^{2}\\\
\text{(\lx@cref{creftype~refnum}{thm:cheby_version}
1.})&=\frac{16}{3}\left(q^{(\ell)}_{d-1-2\ell}-q^{(\ell-1)}_{d+1-2\ell}\right)q^{(\ell-1)}_{d+1-2\ell}\\\
\text{(\lx@cref{creftype~refnum}{thm:cheby_version}
4.})&\leq\frac{8}{3}\left(q^{(\ell)}_{d-1-2\ell}-q^{(\ell-1)}_{d+1-2\ell}\right)\left(q^{(\ell)}_{d-1-2\ell}+q^{(\ell-1)}_{d+1-2\ell}\right)\\\
&=\frac{8}{3}\left(\left({q^{(\ell)}_{d-1-2\ell}}\right)^{2}-\left({q^{(\ell-1)}_{d+1-2\ell}}\right)^{2}\right).\\\
\end{split}$

Here, the monotonicity in Theorem 20 implies that
$q^{(\ell)}_{d-1-2\ell}-q^{(\ell-1)}_{d+1-2\ell}\geq 0$.

If $d$ is odd, one can derive an upper bound by telescoping and Lemma 23

$\begin{split}\sum_{\ell=1}^{\hat{d}}\left\lvert\sin(2\phi_{\ell})\right\rvert^{2}&\leq\frac{8}{3}\sum_{\ell=1}^{\hat{d}}\left(\left(q^{(\ell)}_{d-1-2\ell}\right)^{2}-\left(q^{(\ell-1)}_{d+1-2\ell}\right)^{2}\right)\\\
&\leq\frac{8}{3}\left(1-\left({q_{d-1}^{(0)}}\right)^{2}\right)\leq\frac{8}{3}\left\lVert
f\right\rVert_{\infty}^{2}.\end{split}$

Therefore, we have

$\begin{split}\left\lVert\widetilde{\Phi}^{*}-\widetilde{\Phi}^{0}\right\rVert_{2}^{2}&=\left(\phi_{0}-\frac{\pi}{4}\right)^{2}+\sum_{\ell=1}^{\hat{d}}\phi_{\ell}^{2}\\\
&\leq\frac{1}{4}\frac{\pi^{2}}{4}\left\lvert\sin(2\phi_{0}-\frac{\pi}{2})\right\rvert^{2}+\frac{1}{4}\frac{\pi^{2}}{4}\sum_{\ell=1}^{\hat{d}}\left\lvert\sin(2\phi_{\ell})\right\rvert^{2}\\\
&\leq\frac{1}{4}\frac{\pi^{2}}{4}\left(\frac{8}{3}\left\lVert
f\right\rVert_{\infty}^{2}+\frac{8}{3}\left\lVert
f\right\rVert_{\infty}^{2}\right)=\frac{1}{3}\pi^{2}\left\lVert
f\right\rVert_{\infty}^{2}.\end{split}$

If $d$ is even, we have to treat $\phi_{\widetilde{d}-1}$ separately. Using
Theorem 20 statement 3, we have

$\sin^{2}\left(\phi_{\widetilde{d}-1}\right)=1-\cos^{2}\left(\phi_{\widetilde{d}-1}\right)=1-\left(q^{(\widetilde{d}-2)}_{1}\right)^{2}.$

Recall that $\left\lvert\phi_{\widetilde{d}-1}\right\rvert<\frac{\pi}{4}$, we
apply a finer estimate $\left\lvert
x\right\rvert\leq\frac{\pi}{2\sqrt{2}}\left\lvert\sin(x)\right\rvert$ for any
$\left\lvert x\right\rvert\leq\frac{\pi}{4}$ and get

$\left\lvert\phi_{\widetilde{d}-1}\right\rvert^{2}\leq\frac{\pi^{2}}{8}\sin^{2}\left(\phi_{\widetilde{d}-1}\right).$
(79)

Similarly,

$\sum_{\ell=1}^{\hat{d}}\left\lvert\sin(2\phi_{\ell})\right\rvert^{2}\leq\frac{8}{3}\sum_{\ell=1}^{\hat{d}}\left(\left(q^{(\ell)}_{d-1-2\ell}\right)^{2}-\left(q^{(\ell-1)}_{d+1-2\ell}\right)^{2}\right)=\frac{8}{3}\left(q^{(\widetilde{d}-2)}_{1}\right)^{2}-\frac{8}{3}\left(q_{d-1}^{(0)}\right)^{2}.$

Furthermore, combining the derived results, one has

$\begin{split}\left\lVert\widetilde{\Phi}^{*}-\widetilde{\Phi}^{0}\right\rVert_{2}^{2}&=\left(\phi_{0}-\frac{\pi}{4}\right)^{2}+\sum_{\ell=1}^{\hat{d}}\phi_{\ell}^{2}+\phi_{\widetilde{d}-1}^{2}\\\
&\leq\frac{1}{4}\frac{\pi^{2}}{4}\left\lvert\sin(2\phi_{0})\right\rvert^{2}+\frac{1}{4}\frac{\pi^{2}}{4}\sum_{\ell=1}^{\hat{d}}\left\lvert\sin(2\phi_{\ell})\right\rvert^{2}+\frac{\pi^{2}}{8}\left\lvert\sin\left(\phi_{\widetilde{d}-1}\right)\right\rvert^{2}\\\
&\leq\frac{\pi^{2}}{4}\left(\frac{2}{3}\left\lVert
f\right\rVert_{\infty}^{2}+\frac{2}{3}\left(q^{(\widetilde{d}-2)}_{1}\right)^{2}-\frac{2}{3}\left(q_{d-1}^{(0)}\right)^{2}+\frac{1}{2}-\frac{1}{2}\left(q^{(\widetilde{d}-2)}_{1}\right)^{2}\right)\\\
(\text{\lx@cref{creftype~refnum}{thm:cheby_version}
4})&\leq\frac{\pi^{2}}{4}\left(\frac{2}{3}\left\lVert
f\right\rVert_{\infty}^{2}+\frac{2}{3}-\frac{2}{3}\left(q_{d-1}^{(0)}\right)^{2}\right)\\\
(\text{\lx@cref{creftype~refnum}{lma:supplement}
1})&\leq\frac{\pi^{2}}{4}\left(\frac{2}{3}\left\lVert
f\right\rVert_{\infty}^{2}+\frac{2}{3}\left\lVert
f\right\rVert_{\infty}^{2}\right)=\frac{1}{3}\pi^{2}\left\lVert
f\right\rVert_{\infty}^{2}.\end{split}$

In conclusion, we have

$\left\lVert\widetilde{\Phi}^{*}-\widetilde{\Phi}^{0}\right\rVert_{2}\leq\frac{\pi}{\sqrt{3}}\left\lVert
f\right\rVert_{\infty}.$

∎

##  6 Local strong convexity

###  6.1 Hessian matrix at the maximal solution is positive definite

We first show that the Hessian matrix at the initial point
$\widetilde{\Phi}^{0}$ is positive definite.

######  Theorem 24.

At the initial point
$\widetilde{\Phi}^{0}:=\left(\frac{\pi}{4},0,\cdots,0\right)\in\mathbb{R}^{\widetilde{d}}$,
the Hessian matrix of the optimization problem Eq. 5 is

$\mathrm{Hess}(\widetilde{\Phi}^{0})=\frac{2}{\widetilde{d}}A(\widetilde{\Phi}^{0})^{\top}A(\widetilde{\Phi}^{0})=\left\\{\begin{array}[]{ll}4I&,\text{if
}d\text{ is odd},\\\ \mathrm{diag}\\{4,\cdots,4,2\\}&,\text{if }d\text{ is
even}.\end{array}\right.$ (80)

###### Proof.

According to [5, Lemma 3],

$\mathrm{Re}\left[\langle
0|U(x,\widetilde{\Phi})|0\rangle\right]=-\mathrm{Im}\left[\langle
0|e^{-\mathrm{i}\frac{\pi}{4}Z}U(x,\widetilde{\Phi})e^{-\mathrm{i}\frac{\pi}{4}Z}|0\rangle\right],$

which suggests us to extract $\frac{\pi}{4}$ from the phase factors at the
edges and consider $e^{-\mathrm{i}\frac{\pi}{4}Z}U(x,\widetilde{\Phi}^{0})\\\
e^{-\mathrm{i}\frac{\pi}{4}Z}=W(x)^{d}$. The integer power of the matrix
$W(x):=e^{\mathrm{i}\arccos(x)X}$ exactly generates the Chebyshev polynomial
of the first kind as its upper-left matrix element. Furthermore, its upper-
left is real. The claim follows the property of the exponential map of Pauli
matrices
$W(x)^{d}=e^{\mathrm{i}d\arccos(x)X}=\cos(d\arccos(x))+\mathrm{i}\sin(d\arccos(x))X$.
Then, noting that $\braket{0}{X}{0}=0$, it yields
$\braket{0}{W(x)^{d}}{0}=\cos(d\arccos(x))=T_{d}(x)$.

To proceed, let us consider the monomial
$M^{d_{1},d_{2},d_{3}}(x):=W(x)^{d_{1}}(\mathrm{i}Z)W(x)^{d_{2}}(\mathrm{i}Z)W(x)^{d_{3}}$
for any $d_{1},d_{2},d_{3}\geq 0$. We will show that the upper-left element
$\braket{0}{M^{d_{1},d_{2},d_{3}}(x)}{0}$ is real. Following the
anticommutation relation $ZXZ=-X$, a useful equality can be derived
$ZW(x)Z=e^{\mathrm{i}\arccos(x)ZXZ}=e^{-\mathrm{i}\arccos(x)X}=W(x)^{-1}$.
Consequentially, by moving $Z$ gates to the same site for cancellation, it can
be shown that $\braket{0}{M^{d_{1},d_{2},d_{3}}(x)}{0}=-T_{\left\lvert
d_{1}+d_{3}-d_{2}\right\rvert}(x)$ is real.

Note that taking second-order derivative on the unitary defined in Equation 5
at $\Phi^{0}$ results in a integer-coefficient superposition of monomials of
the defined form. To illustrate, given an odd integer $d=2\widetilde{d}-1$ and
$0\leq i<j<\widetilde{d}$, one can verify that

$\frac{\partial^{2}}{\partial\phi_{i}\partial\phi_{j}}U(x,\widetilde{\Phi})\bigg{|}_{\widetilde{\Phi}_{0}}=M^{i,j-i,d-j}(x)+M^{i,d-j-i,j}(x)+M^{d-j,j-i,i}(x)+M^{j,d-j-i,i}(x).$

Then, by taking the imaginary component, the second-order derivative $g_{ij}$
is vanishing at $\widetilde{\Phi}^{0}$. Thus, the second term in the right-
hand side of Equation 38 vanishes since $g_{ij}(x,\widetilde{\Phi}^{0})=0$.

When $d=2(\widetilde{d}-1)$ is even,
$g_{i}(x,\widetilde{\Phi}^{0})=-2\mathrm{Im}\left[\mathrm{i}W(x)^{i}ZW(x)^{d-i}\right]=-2T_{d-2i}(x)$
for $i=0,\cdots,\widetilde{d}-2$ and
$g_{\widetilde{d}-1}(x,\widetilde{\Phi}^{0})=-1$. According to the discrete
orthogonality of Chebyshev nodes, we have

$\sum_{k=1}^{\widetilde{d}}g_{i}(x_{k},\widetilde{\Phi}^{0})g_{j}(x_{k},\widetilde{\Phi}^{0})=\widetilde{d}\left(2\delta_{ij}-\delta_{i,\widetilde{d}-1}\delta_{j,\widetilde{d}-1}\right).$

When $d=2\widetilde{d}-1$ is odd, we have similarly
$\sum_{k=1}^{\widetilde{d}}g_{i}(x_{k},\widetilde{\Phi}^{0})g_{j}(x_{k},\widetilde{\Phi}^{0})=2\widetilde{d}\delta_{ij}$.

Following Eq. 38, the Hessian matrix at $\widetilde{\Phi}^{0}$ takes the form
in the theorem which completes the proof. ∎

To show that the Hessian matrix is bounded from below, we need the following
lemma regarding the perturbation analysis to the singular values of a matrix.

######  Lemma 25 ([7, Corollary 8.6.2]).

Given A, A+E $\in\mathbb{R}^{m\times n}$ with $m\geq n$, we have

$\left\lvert\sigma_{k}(A+E)-\sigma_{k}(A)\right\rvert\leq\sigma_{0}(E)=\left\lVert
E\right\rVert_{2},\quad\forall k\in[n].$

Then we bound the perturbation to the Jacobian matrix via that of the phase
factors.

######  Lemma 26.

Given any symmetric phase factors $\Phi$ of length $d+1$, let
$E:=A(\widetilde{\Phi}+\bm{\epsilon})-A(\widetilde{\Phi})$. Then $\left\lVert
E\right\rVert_{2}\leq
4\widetilde{d}^{\frac{3}{2}}\left\lVert\bm{\epsilon}\right\rVert_{2}$.

###### Proof.

Recall that in the proof of Theorem 24, we have shown that taking the
derivative w.r.t. $\phi_{i}$ yields terms by adding an $\mathrm{i}Z$
multiplication left to $e^{\mathrm{i}\phi_{i}Z}$. Note that
$\mathrm{i}Z=e^{\mathrm{i}\frac{\pi}{2}Z}$, it shifts the $i$-th site by
$\frac{\pi}{2}$ which essentially breaks the symmetry constraint. That means
we have to use asymmetric set of phase factors to characterize the derivative.
To capture this, we define the evaluation map of asymmetric phase factors as
follows. Given an asymmetric set of phase factors
$\Phi_{\text{as}}:=(\phi_{0},\phi_{1},\cdots,\phi_{d})$, we define

$U^{\text{(as)}}(x,\Phi_{\text{as}}):=e^{\mathrm{i}\phi_{0}Z}e^{\mathrm{i}\arccos(x)X}e^{\mathrm{i}\phi_{1}Z}e^{\mathrm{i}\arccos(x)X}\cdots
e^{\mathrm{i}\phi_{d-1}Z}e^{\mathrm{i}\arccos(x)X}e^{\mathrm{i}\phi_{d}Z}.$

We use a superscript (as) to distinguish it with $U(x,\widetilde{\Phi})$ which
takes the reduced phase factors $\widetilde{\Phi}$ of a symmetric set as
argument.

When $d$ is odd, we define an asymmetric set of phase factors as

$\Phi_{\text{as}}^{(i)}:=\left(\phi_{0},\cdots,\phi_{i-1},\phi_{i}+\frac{\pi}{2},\phi_{i+1},\cdots,\phi_{\widetilde{d}-1},\phi_{\widetilde{d}-1},\cdots,\phi_{i+1},\phi_{i},\phi_{i-1},\cdots,\phi_{0}\right)$

which unrolls the full set of phase factors and add a $\frac{\pi}{2}$ shift on
the $i$-th site. As a remark, the shifted site $i$ can be greater than
$\widetilde{d}$ which adds a phase shift on the right half.

Taking derivative yields that

$\begin{split}g_{i}(x,\widetilde{\Phi})&=\mathrm{Re}\left[\braket{0}{\frac{\partial}{\partial\phi_{i}}U(x,\widetilde{\Phi})}{0}\right]=\mathrm{Re}[\braket{0}{U^{\text{(as)}}(x,\Phi_{\text{as}}^{(i)})}{0}]+\mathrm{Re}[\braket{0}{U^{\text{(as)}}(x,\Phi_{\text{as}}^{(d-i)})}{0}]\\\
&=\mathrm{Re}[\braket{0}{U^{\text{(as)}}(x,\Phi_{\text{as}}^{(i)})}{0}]+\mathrm{Re}[\braket{0}{U^{\text{(as)}}(x,\Phi_{\text{as}}^{(i)})^{\top}}{0}]=2\mathrm{Re}[\braket{0}{U^{\text{(as)}}(x,\Phi_{\text{as}}^{(i)})}{0}].\end{split}$

Due to the unitarity, we have

$\left\lvert g_{i}(x,\widetilde{\Phi})\right\rvert\leq 2,\quad\forall
i\in[\widetilde{d}].$

To further take the second derivative, let us first introduce a double-shifted
asymmetric set of phase factors. When $0\leq i<j<\widetilde{d}$,

$\Phi_{\text{as}}^{(i,j)}:=\left(\phi_{0},\cdots,\phi_{i-1},\phi_{i}+\frac{\pi}{2},\phi_{i+1},\cdots,\phi_{j-1},\phi_{j}+\frac{\pi}{2},\phi_{j+1},\cdots,\phi_{\widetilde{d}-1},\phi_{\widetilde{d}-1},\cdots,\phi_{0}\right).$

When $0\leq i=j<\widetilde{d}$,

$\Phi_{\text{as}}^{(i,i)}:=\left(\phi_{0},\cdots,\phi_{i-1},\phi_{i}+\pi,\phi_{i+1},\cdots,\phi_{\widetilde{d}-1},\phi_{\widetilde{d}-1},\cdots,\phi_{i+1},\phi_{i},\phi_{i-1},\cdots,\phi_{0}\right).$

Similarly, when $i\text{ or }j\geq\widetilde{d}$, the phase shift is added on
the right half accordingly.

With the defined notation, we are able to compute the second derivative as

$\begin{split}g_{ij}(x,\widetilde{\Phi})&:=\frac{\partial}{\partial\phi_{j}}g_{i}(x,\widetilde{\Phi})=2\mathrm{Re}[\braket{0}{\frac{\partial}{\partial\phi_{j}}U^{\text{(as)}}(x,\Phi_{\text{as}}^{(i)})}{0}]\\\
&=2\mathrm{Re}[\braket{0}{U^{\text{(as)}}(x,\Phi_{\text{as}}^{(i,j)})}{0}]+2\mathrm{Re}[\braket{0}{U^{\text{(as)}}(x,\Phi_{\text{as}}^{(i,d-j)})}{0}].\end{split}$

Due to the unitarity, we have

$\left\lvert g_{ij}(x,\widetilde{\Phi})\right\rvert\leq 4,\quad\forall
i,j\in[\widetilde{d}].$

Therefore, for any given odd integer $d$ and symmetric
$\Phi\in\mathbb{R}^{d+1}$, we have

$\left\lVert\nabla
g_{i}(x,\widetilde{\Phi})\right\rVert_{2}=\sqrt{\sum_{j\in[\widetilde{d}]}\left\lvert
g_{ij}(x,\widetilde{\Phi})\right\rvert^{2}}\leq
4\sqrt{\widetilde{d}},\quad\forall i\in[\widetilde{d}].$ (81)

When $d$ is even, the full set of phase factors is
$(\phi_{0},...\phi_{\widetilde{d}-2},\phi_{\widetilde{d}-1},\phi_{\widetilde{d}-2},...,\phi_{0})$.
The only difference in the unrolled structure comparing to that in the odd
case is that there is a single unpaired phase factor $\phi_{\widetilde{d}-1}$
at the middle. Therefore, when the derivative only involves
$\phi_{0},\cdots,\phi_{\widetilde{d}-2}$, the analysis is exactly the same as
that when $d$ is odd.

When $i=\widetilde{d}-1$, taking derivative w.r.t. $\phi_{\widetilde{d}-1}$ is
equivalently shifting this site by $\frac{\pi}{2}$, namely,

$g_{\widetilde{d}-1}(x,\widetilde{\Phi})=\mathrm{Re}[\braket{0}{U(x,\widetilde{\Psi})}{0}]=g(x,\widetilde{\Psi}),$

where

$\Psi:=\Phi+\frac{\pi}{2}\mathbf{e}_{\widetilde{d}-1}=(\phi_{0},\cdots,\phi_{\widetilde{d}-2},\phi_{\widetilde{d}-1}+\frac{\pi}{2},\phi_{\widetilde{d}-2},\cdots,\phi_{0}).$

Note that $\Psi$ is a set of symmetric phase factors (and therefore we use a
different notation). Then, taking the second derivative on
$g(x,\widetilde{\Phi})$ is equivalently to taking the first derivative on
$g(x,\widetilde{\Psi})$ whose analysis is the same as that in the previous
case. Therefore,

$g_{i,\widetilde{d}-1}(x,\widetilde{\Phi})=g_{i}(x,\widetilde{\Psi})=2\mathrm{Re}\left[\braket{0}{U^{\text{(as)}}(x,\Psi_{\text{as}}^{(i)})}{0}\right],\
\forall i\in[\widetilde{d}-1],$

and

$g_{\widetilde{d}-1,\widetilde{d}-1}(x,\widetilde{\Phi})=g_{\widetilde{d}-1}(x,\widetilde{\Psi})=-g(x,\widetilde{\Phi}).$

These results lead to the inequality

$\left\lvert g_{i,\widetilde{d}-1}(x,\widetilde{\Phi})\right\rvert\leq
2,\quad\forall i\in[\widetilde{d}].$

Thus, the inequality in Eq. 81 also holds when $d$ is even.

To conclude, for any given integer $d$ and symmetric
$\Phi\in\mathbb{R}^{d+1}$, we obtain a uniform upper bound of the gradient of
$g_{j}(x,\widetilde{\Phi})$

$\left\lVert\nabla g_{j}(x,\widetilde{\Phi})\right\rVert_{2}\leq
4\sqrt{\widetilde{d}},\ \forall j\in[\widetilde{d}],\ \forall x\in[-1,1].$
(82)

Following the mean value theorem, there exists $t\in(0,1)$ so that for any
$j\in[\widetilde{d}]$ and $x\in[-1,1]$

$\left\lvert
g_{j}(x,\widetilde{\Phi}+\bm{\epsilon})-g_{j}(x,\widetilde{\Phi})\right\rvert\leq\left\lVert\nabla
g_{j}(x,\widetilde{\Phi}+t\bm{\epsilon})\right\rVert_{2}\left\lVert\bm{\epsilon}\right\rVert_{2}\leq
4\sqrt{\widetilde{d}}\left\lVert\bm{\epsilon}\right\rVert_{2}.$

Furthermore

$\left\lVert
A(\widetilde{\Phi}+\bm{\epsilon})-A(\widetilde{\Phi})\right\rVert_{F}=\sqrt{\sum_{j=0}^{\widetilde{d}-1}\sum_{k=1}^{\widetilde{d}}\left(g_{j}(x_{k},\widetilde{\Phi}+\bm{\epsilon})-g_{j}(x_{k},\widetilde{\Phi})\right)^{2}}\leq
4\widetilde{d}^{\frac{3}{2}}\left\lVert\bm{\epsilon}\right\rVert_{2}.$

Then, the theorem follows

$\left\lVert
A(\widetilde{\Phi}+\bm{\epsilon})-A(\widetilde{\Phi})\right\rVert_{2}\leq\left\lVert
A(\widetilde{\Phi}+\bm{\epsilon})-A(\widetilde{\Phi})\right\rVert_{F}\leq
4\widetilde{d}^{\frac{3}{2}}\left\lVert\bm{\epsilon}\right\rVert_{2}.$

∎

Then the perturbation analysis to the singular value of the Jacobian matrix is
obtained by directly applying Lemma 25.

######  Corollary 27.

If $\left\lVert\bm{\epsilon}\right\rVert_{2}\leq\frac{1}{20\widetilde{d}}$,
then

$\frac{4}{5}\sqrt{\widetilde{d}}\leq\sigma_{\min}(A(\widetilde{\Phi}^{0}+\bm{\epsilon}))\leq\sigma_{\max}(A(\widetilde{\Phi}^{0}+\bm{\epsilon}))\leq\left(\frac{1}{5}+\sqrt{2}\right)\sqrt{\widetilde{d}}.$

###### Proof.

Applying Lemma 25, we have

$\left\lvert\sigma_{k}(A(\widetilde{\Phi}^{0}+\bm{\epsilon}))-\sigma_{k}(A(\widetilde{\Phi}^{0}))\right\rvert\leq
4\widetilde{d}^{\frac{3}{2}}\left\lVert\bm{\epsilon}\right\rVert_{2}.$ (83)

As a consequence of Theorem 24,
$\sigma_{\max}(A(\widetilde{\Phi}^{0}))=\sqrt{2\widetilde{d}}$ and

$\sigma_{\text{min}}(A(\widetilde{\Phi}^{0})=\left\\{\begin{array}[]{ll}\sqrt{\widetilde{d}}&,\text{
when }d\text{ is even},\\\ \sqrt{2\widetilde{d}}&,\text{ when }d\text{ is
odd}.\end{array}\right.$

When $\left\lVert\bm{\epsilon}\right\rVert_{2}\leq\frac{1}{20\widetilde{d}}$,
we have

$\sigma_{\min}(A(\widetilde{\Phi}^{0}+\bm{\epsilon}))\geq\sigma_{\min}(A(\widetilde{\Phi}^{0}))-\left\lVert
E\right\rVert_{2}\geq\sqrt{\widetilde{d}}-4\widetilde{d}^{\frac{3}{2}}\left\lVert\bm{\epsilon}\right\rVert_{2}\geq\frac{4}{5}\sqrt{\widetilde{d}}.$
(84)

and

$\sigma_{\max}(A(\widetilde{\Phi}^{0}+\bm{\epsilon}))\leq\sigma_{\max}(A(\widetilde{\Phi}^{0}))+\left\lVert
E\right\rVert_{2}\leq\sqrt{2\widetilde{d}}+4\widetilde{d}^{\frac{3}{2}}\left\lVert\bm{\epsilon}\right\rVert_{2}\leq\left(\frac{1}{5}+\sqrt{2}\right)\sqrt{\widetilde{d}}.$
(85)

∎

Recall that the Hessian matrix at the maximal solution satisfies
$\text{Hess}(\widetilde{\Phi}^{*})=\frac{2}{\widetilde{d}}A(\widetilde{\Phi}^{*})^{\top}A(\widetilde{\Phi}^{*})$.
The positive definiteness of this matrix follows by combining the previous
results with Theorem 5.

######  Theorem 28.

Let $\Phi^{*}$ be the optimal phase factors corresponding to the maximal
solution. If $\left\lVert
f\right\rVert_{\infty}\leq\frac{\sqrt{3}}{20\pi\widetilde{d}}$, then
$\mathrm{Hess}(\widetilde{\Phi}^{*})$ is positive definite and

$\lambda_{\min}\left(\mathrm{Hess}(\widetilde{\Phi}^{*})\right)\geq\frac{32}{25}.$
(86)

###### Proof.

Theorem 5 implies that
$\left\lVert\widetilde{\Phi}^{*}-\widetilde{\Phi}^{0}\right\rVert_{2}\leq\frac{1}{20\widetilde{d}}$.
As a consequence of Corollary 27, we have
$\sigma_{\text{min}}(A(\widetilde{\Phi}^{*}))\geq\frac{4}{5}\sqrt{\widetilde{d}}$.
Therefore

$\lambda_{\text{min}}(\text{Hess}(\widetilde{\Phi}^{*}))=\lambda_{\text{min}}\left(\frac{2}{\widetilde{d}}A(\widetilde{\Phi}^{*})^{\top}A(\widetilde{\Phi}^{*})\right)=\frac{2}{\widetilde{d}}\sigma_{\text{min}}^{2}(A(\widetilde{\Phi}^{*}))\geq\frac{32}{25},$

which completes the proof. ∎

###  6.2 Lower bounding the distance to the initial point

In the last subsection, we show that the Hessian matrix is positive definite
at the maximal solution. To show the positive definiteness of the Hessian
matrix in a neighborhood of the optima, we need the following lower bound as
an intermediate step.

######  Theorem 29.

Given any symmetric phase factor $\widetilde{\Phi}$ satisfying
$\left\lVert\widetilde{\Phi}-\widetilde{\Phi}^{0}\right\rVert_{2}\leq\frac{1}{\sqrt{2}}$,
the following inequality holds

$\frac{1}{\sqrt{3}}\left\lVert
g(x,\widetilde{\Phi})\right\rVert_{T}\leq\left\lVert\widetilde{\Phi}-\widetilde{\Phi}^{0}\right\rVert_{2}.$
(87)

###### Proof.

We denote the expansion of the polynomial $g(x,\widetilde{\Phi})$ in Chebyshev
basis as

$g(x,\widetilde{\Phi})=\sum_{k=0}^{d}g_{k}(\widetilde{\Phi})T_{k}(x).$

For simplicity, we drop the $\widetilde{\Phi}$-dependence in the Chebyshev
coefficient $g_{k}$ (in particular, $g_{k}$ is not the derivative of
$g(x,\widetilde{\Phi})$ with respect to $\phi_{k}$ as in the previous
section).

Then we apply Eq. 26 to estimate $q_{d-1}$, the coefficient of $Q(x)$. Notice
that $\cos(\phi_{\ell})$ is positive for any $1\leq\ell\leq\widetilde{d}-1$,
which follows
$\left\lvert\phi_{\ell}\right\rvert\leq\left\lVert\widetilde{\Phi}-\widetilde{\Phi}^{0}\right\rVert_{2}\leq\frac{1}{\sqrt{2}}<\frac{\pi}{2}$.
We use Eq. 26 and the inequality
$\cos(\phi_{\widetilde{d}-1})\geq\cos(\phi_{\widetilde{d}-1})^{2}$ when $d$ is
even, and obtain

$q_{d-1}^{(0)}\geq\prod_{\ell=1}^{\widetilde{d}-1}\cos^{2}(\phi_{\ell})=\prod_{\ell=1}^{\widetilde{d}-1}\left(1-\sin^{2}(\phi_{\ell})\right)\geq
1-\sum_{\ell=1}^{\widetilde{d}-1}\phi_{\ell}^{2}\geq
1-\left\lVert\widetilde{\Phi}-\widetilde{\Phi}^{0}\right\rVert_{2}^{2}\geq\frac{1}{2}.$
(88)

Hence $\deg(Q)=d-1$ and the leading coefficient of $Q$ is positive. Now we are
able to apply Theorems 20 and 21 and have

$\sqrt{1-\left\lVert g(x,\widetilde{\Phi})\right\rVert_{\infty}^{2}}\leq
q^{(0)}_{d-1}\leq q^{(1)}_{d-3}\leq\cdots\leq
q^{(\ell)}_{d-1-2\ell}\leq\cdots\leq q^{(\hat{d})}_{d-1-2\hat{d}}\leq 1.$

According to the reduction procedure, the leading Chebyshev coefficients and
the phase factor are connected as follows

$\left\lvert\sin(2\phi_{0}-\frac{\pi}{2})\right\rvert^{2}=\left\lvert\frac{\mathrm{Re}[p_{d}]}{q_{d-1}}\right\rvert^{2}=\frac{g_{d}^{2}}{q_{d-1}^{2}}\geq
g_{d}^{2},$ (89)

where the last inequality follows $\left\lvert q_{d-1}\right\rvert\leq 1$.

Following Eq. 36, for any $1\leq\ell\leq\hat{d}$, we have

$\begin{split}\left\lvert\sin(2\phi_{\ell})\right\rvert^{2}&=\left\lvert\frac{\mathrm{Im}[p^{(\ell-1)}_{d-2\ell}e^{-2\mathrm{i}\phi_{\ell-1}}]}{q_{d-1-2\ell}^{(\ell)}}\right\rvert^{2}\geq\left\lvert\mathrm{Im}[p^{(\ell-1)}_{d-2\ell}e^{-2\mathrm{i}\phi_{\ell-1}}]\right\rvert^{2}\\\
&=4q^{(\ell-1)}_{d+1-2\ell}\left(q^{(\ell)}_{d-1-2\ell}-q^{(\ell-1)}_{d+1-2\ell}\right)\geq
4q_{d-1}^{(0)}\left(q^{(\ell)}_{d-1-2\ell}-q^{(\ell-1)}_{d+1-2\ell}\right),\\\
\end{split}$

where the last inequality follows the monotonicity in Theorem 20.

Now we are going to show

$\left\lVert\widetilde{\Phi}-\widetilde{\Phi}^{0}\right\rVert_{2}^{2}\geq\frac{1}{4}g_{d}^{2}+q_{d-1}^{(0)}\left(1-{q_{d-1}^{(0)}}\right),$
(90)

where $g_{d}$ denotes the leading Chebyshev coefficient of
$g(x,\widetilde{\Phi})$.

When $d$ is odd, we have

$\sum_{\ell=1}^{\hat{d}}\left\lvert\sin(2\phi_{\ell})\right\rvert^{2}\geq
4q_{d-1}^{(0)}\sum_{\ell=1}^{\hat{d}}\left(q^{(\ell)}_{d-1-2\ell}-q^{(\ell-1)}_{d+1-2\ell}\right)=4q_{d-1}^{(0)}\left(q_{0}^{(\hat{d})}-q_{d-1}^{(0)}\right)=4q_{d-1}^{(0)}\left(1-{q_{d-1}^{(0)}}\right).$

Here, we have used the fact that $q_{0}^{(\hat{d})}=1$ obtained from the proof
of Lemma 13. Hence

$\begin{split}\left\lVert\widetilde{\Phi}-\widetilde{\Phi}^{0}\right\rVert_{2}^{2}&=\left(\phi_{0}-\frac{\pi}{4}\right)^{2}+\sum_{\ell=1}^{\hat{d}}\phi_{\ell}^{2}\\\
&\geq\frac{1}{4}\left\lvert\sin(2\phi_{0}-\frac{\pi}{2})\right\rvert^{2}+\frac{1}{4}\sum_{\ell=1}^{\hat{d}}\left\lvert\sin(2\phi_{\ell})\right\rvert^{2}\\\
&\geq\frac{1}{4}g_{d}^{2}+q_{d-1}^{(0)}\left(1-{q_{d-1}^{(0)}}\right).\end{split}$

When $d$ is even, the dangling unpaired phase factor $\phi_{\widetilde{d}-1}$
has to be analyzed individually

$\left\lvert\sin\left(\phi_{\widetilde{d}-1}\right)\right\rvert^{2}=1-\left({q^{(\widetilde{d}-2)}_{1}}\right)^{2}.$

Furthermore, we have

$\sum_{\ell=1}^{\hat{d}}\left\lvert\sin(2\phi_{\ell})\right\rvert^{2}\geq
4q_{d-1}^{(0)}\sum_{\ell=1}^{\hat{d}}\left(q^{(\ell)}_{d-1-2\ell}-q^{(\ell-1)}_{d+1-2\ell}\right)=4q_{d-1}^{(0)}\left(q_{1}^{(\widetilde{d}-2)}-q_{d-1}^{(0)}\right).$

Unlike the odd case, the coefficient $q_{1}^{(\widetilde{d}-2)}$ is not
necessary equal to $1$ when $d$ is even. Then

$\begin{split}\left\lVert\widetilde{\Phi}-\widetilde{\Phi}^{0}\right\rVert_{2}^{2}&=\left(\phi_{0}-\frac{\pi}{4}\right)^{2}+\sum_{\ell=1}^{\hat{d}}\phi_{\ell}^{2}+\phi_{\widetilde{d}-1}^{2}\\\
&\geq\frac{1}{4}\left\lvert\sin(2\phi_{0})\right\rvert^{2}+\frac{1}{4}\sum_{\ell=1}^{\hat{d}}\left\lvert\sin(2\phi_{\ell})\right\rvert^{2}+\left\lvert\sin\left(\phi_{\widetilde{d}-1}\right)\right\rvert^{2}\\\
&\geq\frac{1}{4}g_{d}^{2}+q_{d-1}^{(0)}\left(q_{1}^{(\widetilde{d}-2)}-q_{d-1}^{(0)}\right)+1-\left(q^{(\widetilde{d}-2)}_{1}\right)^{2}\\\
&=\frac{1}{4}g_{d}^{2}+q_{d-1}^{(0)}\left(1-{q_{d-1}^{(0)}}\right)+\left(1-{q_{d-1}^{(0)}}\right)\left(1-{q_{1}^{(\widetilde{d}-2)}}\right)+q_{1}^{(\widetilde{d}-2)}\left(1-{q_{1}^{(\widetilde{d}-2)}}\right)\\\
&\geq\frac{1}{4}g_{d}^{2}+q_{d-1}^{(0)}\left(1-{q_{d-1}^{(0)}}\right).\end{split}$

According to Eqs. 88 and 21, we have $\frac{1}{2}\leq
q_{d-1}^{(0)}\leq\sqrt{1-\left\lVert
g(x,\widetilde{\Phi})\right\rVert_{T}^{2}+\frac{1}{2}g_{d}^{2}}$. It implies
that

$\begin{split}\left\lVert\widetilde{\Phi}-\widetilde{\Phi}^{0}\right\rVert_{2}^{2}&\geq\frac{1}{4}g_{d}^{2}+\sqrt{1-\left\lVert
g(x,\widetilde{\Phi})\right\rVert_{T}^{2}+\frac{1}{2}g_{d}^{2}}-1+\left\lVert
g(x,\widetilde{\Phi})\right\rVert_{T}^{2}-\frac{1}{2}g_{d}^{2}\\\
&\geq\left\lVert
g(x,\widetilde{\Phi})\right\rVert_{T}^{2}-\frac{1}{4}g_{d}^{2}-\frac{\left\lVert
g(x,\widetilde{\Phi})\right\rVert_{T}^{2}-\frac{1}{2}g_{d}^{2}}{\sqrt{1-\left\lVert
g(x,\widetilde{\Phi})\right\rVert_{T}^{2}+\frac{1}{2}g_{d}^{2}}+1}\\\
&\geq\left\lVert
g(x,\widetilde{\Phi})\right\rVert_{T}^{2}-\frac{1}{4}g_{d}^{2}-\frac{2\left\lVert
g(x,\widetilde{\Phi})\right\rVert_{T}^{2}-g_{d}^{2}}{3}\\\
&\geq\frac{1}{3}\left\lVert
g(x,\widetilde{\Phi})\right\rVert_{T}^{2}.\end{split}$

∎

###  6.3 Local strong convexity of the optimization problem

Now we are going to present the proof of Theorem 6, where we use Corollary 27
and the following lemma regarding the perturbative analysis to the eigenvalues
of a matrix.

######  Lemma 30 ([7, Theorem 8.1.4] Wielandt–Hoffman).

If A, A+E $\in\mathbb{R}^{n\times n}$ are symmetric matrices, then

$\sum_{k\in[n]}\left(\lambda_{k}(A+E)-\lambda_{k}(A)\right)^{2}\leq\left\lVert
E\right\rVert_{F}^{2}.$

###### Proof of Theorem 6.

Recall that

$F(\widetilde{\Phi})=\frac{2}{\widetilde{d}}\left(A^{\top}(\widetilde{\Phi})A(\widetilde{\Phi})\right)-R(\widetilde{\Phi}),$
(91)

where
$R_{ij}:=\frac{2}{\widetilde{d}}\sum_{k=1}^{\widetilde{d}}\left(g(x_{k},\widetilde{\Phi})-f(x_{k})\right)g_{ij}(x_{k},\widetilde{\Phi})$.
The residual term is estimated as follows

$\begin{split}\left\lVert
R\right\rVert_{F}^{2}&=\sum_{i=0}^{\widetilde{d}-1}\sum_{j=0}^{\widetilde{d}-1}\left(\frac{2}{\widetilde{d}}\sum_{k=1}^{\widetilde{d}}\left(g(x_{k},\widetilde{\Phi})-f(x_{k})\right)g_{ij}(x_{k},\widetilde{\Phi})\right)^{2}\\\
&\leq\frac{4}{\widetilde{d}^{2}}\sum_{i=0}^{\widetilde{d}-1}\sum_{j=0}^{\widetilde{d}-1}\left(\sum_{k=1}^{\widetilde{d}}\left(g(x_{k},\widetilde{\Phi})-f(x_{k})\right)^{2}\right)\left(\sum_{k=1}^{\widetilde{d}}g_{ij}(x_{k},\widetilde{\Phi})^{2}\right)\\\
&\leq
8\left(\sum_{k=1}^{\widetilde{d}}g(x_{k},\widetilde{\Phi})^{2}+\sum_{k=1}^{\widetilde{d}}f(x_{k})^{2}\right)\max_{i,j}\left(\sum_{k=1}^{\widetilde{d}}g_{ij}(x_{k},\widetilde{\Phi})^{2}\right)\\\
&=8\widetilde{d}^{2}\left(\left\lVert
g(x,\widetilde{\Phi})\right\rVert_{T}^{2}+\left\lVert
f(x)\right\rVert_{T}^{2}\right)\max_{i,j}\left\lVert
g_{ij}(x,\widetilde{\Phi})\right\rVert_{T}^{2}.\end{split}$ (92)

Here, the last line uses the discrete orthogonality of Chebyshev polynomial.

Then, we have the following inequality

$\left\lVert R\right\rVert_{F}^{2}\leq
128\widetilde{d}^{2}\left(3\left\lVert\widetilde{\Phi}-\widetilde{\Phi}^{0}\right\rVert_{2}^{2}+\left\lVert
f\right\rVert_{\infty}^{2}\right)\leq\frac{24}{25}\left(1+\frac{1}{\pi^{2}}\right),$
(93)

where the first inequality uses Theorem 29, the fact $\left\lVert
f\right\rVert_{T}\leq\left\lVert f\right\rVert_{\infty}$, and its consequence
$\left\lVert g_{ij}(x,\widetilde{\Phi})\right\rVert_{T}\leq\left\lVert
g_{ij}(x,\widetilde{\Phi})\right\rVert_{\infty}\leq 4$ derived in the proof of
Lemma 26.

Then, applying Lemma 30, one has

$\begin{split}&\lambda_{\min}(\mathrm{Hess}(\widetilde{\Phi}))\geq\lambda_{\min}(\frac{2}{\widetilde{d}}A(\widetilde{\Phi})^{\top}A(\widetilde{\Phi}))-\left\lVert
R\right\rVert_{F}\geq\frac{32}{25}-\sqrt{\frac{24}{25}\left(1+\frac{1}{\pi^{2}}\right)}\geq\frac{1}{4},\\\
&\lambda_{\max}(\mathrm{Hess}(\widetilde{\Phi}))\leq\lambda_{\max}(\frac{2}{\widetilde{d}}A(\widetilde{\Phi})^{\top}A(\widetilde{\Phi}))+\left\lVert
R\right\rVert_{F}\leq
2\left(\frac{1}{5}+\sqrt{2}\right)^{2}+\sqrt{\frac{24}{25}\left(1+\frac{1}{\pi^{2}}\right)}\leq\frac{25}{4}.\end{split}$
(94)

∎

######  Remark 31.

According to Theorem 6, we know that if $\left\lVert
f\right\rVert_{\infty}\leq\frac{\sqrt{3}}{20\pi\widetilde{d}}$, the
optimization problem Eq. 5 is strongly convex in the following region

$\mathcal{S}:=\left\\{\widetilde{\Phi}\in\mathbb{R}^{\widetilde{d}}:\left\lVert\widetilde{\Phi}-\widetilde{\Phi}^{0}\right\rVert_{2}\leq\frac{1}{20\widetilde{d}}\right\\}.$
(95)

Furthermore, the maximal solution is included in that region
$\widetilde{\Phi}^{*}\in\mathcal{S}$. For convenience, the lower bound of the
Hessian matrix is referred to as $\sigma=\frac{1}{4}$, while the upper bound
is referred to as $L=\frac{25}{4}$.

###  6.4 Projected gradient descend algorithm for symmetric phase-factor
evaluation

Algorithm 1 Projected gradient descend algorithm for optimizing symmetric
phase factors (with performance guarantee)

Input: A target real polynomial $f$ of degree $d$ so that it has definite
parity and $\left\lVert
f\right\rVert_{\infty}\leq\frac{\sqrt{3}}{20\pi\widetilde{d}}$, an initial
guess
$\widetilde{\Phi}^{0}=\left(\frac{\pi}{4},0,\cdots,0\right)\in\mathcal{S}$ and
the error tolerance $\epsilon$.

Output: A set of reduced phase factors $\widetilde{\Phi}$ so that
$F(\widetilde{\Phi})\leq\epsilon$ implying $\left\lVert
g(x,\widetilde{\Phi})-f(x)\right\rVert_{\infty}\leq d\sqrt{\epsilon}$ [5,
Theorem 4].

Evaluate the objective function $F(\widetilde{\Phi})$ at
$\widetilde{\Phi}^{0}$ according to Eq. 5.

while $F(\widetilde{\Phi})>\epsilon$ do

Evaluate the gradient $\nabla F(\widetilde{\Phi})$ of the objective.

Update $\widetilde{\Phi}_{i}$ by gradient descent
$\widetilde{\Phi}\leftarrow\widetilde{\Phi}-\frac{1}{L}\nabla
F(\widetilde{\Phi})$.

if
$\left\lVert\widetilde{\Phi}-\widetilde{\Phi}^{0}\right\rVert_{2}>\frac{1}{20\widetilde{d}}$
then

Project the point onto the ball
$\widetilde{\Phi}\leftarrow\widetilde{\Phi}^{0}-\frac{\widetilde{\Phi}-\widetilde{\Phi}^{0}}{20\widetilde{d}\left\lVert\widetilde{\Phi}-\widetilde{\Phi}^{0}\right\rVert_{2}}$.

end if

end while

return $\widetilde{\Phi}$

We list a practical algorithm for symmetric phase-factor evaluation in
Algorithm 1 which is based on the projected gradient descent method [2]. In
Corollary 7, we give the convergence analysis of this algorithm and its proof
can be obtained by directly applying [2, Theorem 3.10]. Let
$\widetilde{\Phi}^{\ell}$ be the set of reduced phase factors in the $\ell$-th
iteration step. Then, $\widetilde{\Phi}^{\ell}\in\mathcal{S}$ for all $\ell$
follows the construction of the algorithm. Note that
$F(\widetilde{\Phi}^{*})=0,\ \nabla F(\widetilde{\Phi}^{*})=0$ and
$\widetilde{\Phi}^{*}\in\mathcal{S}$. Following Taylor’s theorem, there exists
$s\in(0,1)$ so that

$\begin{split}F(\widetilde{\Phi}^{\ell})&=F(\widetilde{\Phi}^{*})+\nabla
F(\widetilde{\Phi}^{*})^{\top}\left(\widetilde{\Phi}^{\ell}-\widetilde{\Phi}^{*}\right)+\frac{1}{2}\left(\widetilde{\Phi}^{\ell}-\widetilde{\Phi}^{*}\right)^{\top}\text{Hess}\left(s\widetilde{\Phi}^{\ell}+(1-s)\widetilde{\Phi}^{*}\right)\left(\widetilde{\Phi}^{\ell}-\widetilde{\Phi}^{*}\right)\\\
&=\frac{1}{2}\left(\widetilde{\Phi}^{\ell}-\widetilde{\Phi}^{*}\right)^{\top}\text{Hess}\left(s\widetilde{\Phi}^{\ell}+(1-s)\widetilde{\Phi}^{*}\right)\left(\widetilde{\Phi}^{\ell}-\widetilde{\Phi}^{*}\right).\end{split}$

Because $\widetilde{\Phi}^{\ell},\widetilde{\Phi}^{*}\in\mathcal{S}$ and
$\mathcal{S}$ is convex, the intermediate point also lies in the convex set
$s\widetilde{\Phi}^{\ell}+(1-s)\widetilde{\Phi}^{*}\in\mathcal{S}$. Then, the
Hessian is upper bounded at the intermediate point by Theorem 6. Thus

$F(\widetilde{\Phi}^{\ell})\leq\frac{1}{2}\lambda_{\text{max}}\left(\text{Hess}\left(s\widetilde{\Phi}^{\ell}+(1-s)\widetilde{\Phi}^{*}\right)\right)\left\lVert\widetilde{\Phi}^{\ell}-\widetilde{\Phi}^{*}\right\rVert_{2}^{2}\leq\frac{L}{2}\left\lVert\widetilde{\Phi}^{\ell}-\widetilde{\Phi}^{*}\right\rVert_{2}^{2}.$

Furthermore, using Corollary 7, we have

$F(\widetilde{\Phi}^{\ell})\leq\frac{L}{2}\left\lVert\widetilde{\Phi}^{\ell}-\widetilde{\Phi}^{*}\right\rVert_{2}^{2}\leq\frac{L}{2}e^{-\frac{\sigma}{L}\ell}\left\lVert\widetilde{\Phi}^{0}-\widetilde{\Phi}^{*}\right\rVert_{2}^{2}.$

To obtain $F(\widetilde{\Phi}^{\ell})\leq\epsilon$, it suffices to take
$\ell=\mathcal{O}\left(\log\left(\frac{1}{\epsilon}\right)\right)$ iterations.

Remarkably, Algorithm 1 does not require a large number of bits. To see the
claim, we consider the model of finite precision arithmetic (see standard
axioms in [9]). The parameter $u$ is referred to the unit roundoff. In each
iteration step, the rounding error when computing the function $F$ and its
gradient is $\mathcal{O}(d^{2}u)$. Furthermore, the projection is needed to be
executed when the proposed new $\Phi$ lies outside $\mathcal{S}$. Under this
circumstance, the possible error amplification by taking square roots and
division does not contribute significantly and the total rounding error
accumulated in each iteration step remains $\mathcal{O}(d^{2}u)$. To conclude,
the numerical error of the algorithm is upper bounded by invoking triangle
inequality. Then, there exists some constant $C$ so that

$Cd^{2}u\log\frac{1}{\epsilon}\leq\epsilon\ \Rightarrow\
\frac{1}{u}:=2^{n_{\text{bits}}}\geq\frac{1}{C}\frac{d^{2}\log(1/\epsilon)}{\epsilon}\
\Rightarrow\
n_{\text{bits}}=\mathcal{O}\left(\log\left(\frac{d\log(1/\epsilon)}{\epsilon}\right)\right).$

Therefore, the number of bits required by Algorithm 1 is
$\mathcal{O}\left(\log\left(\frac{d\log(1/\epsilon)}{\epsilon}\right)\right)$.
Therefore the algorithm is numerically stable, and it is sufficient to use the
standard double precision arithmetic operations.

##  7 Additional numerical results

In this section, we present some additional numerical results demonstrating
some of the features of the energy landscape of the optimization problem for
symmetric phase factors in Eq. 5.

###  7.1 Hessian matrix without symmetry constraint is singular at global
optima

When the symmetry constraint is not imposed, the Hessian matrix can be derived
from Eq. 38. Note that the Jacobian matrix
$A\in\mathbb{R}^{\widetilde{d}\times d}$ does not have full column rank
because $d>\widetilde{d}$. Then, following that $\ker(A^{\top}A)=\ker(A)$, we
have $\text{rank}(A^{\top}A)=\text{rank}(A)\leq\widetilde{d}<d$. Note that at
the global optima, the Hessian matrix is precisely
$\text{Hess}(\Phi^{*})=(A^{*})^{\top}A^{*}$ due to the vanishing residual term
$R=0$. Here, $A^{*}$ refers to the Jacobian matrix at the optima. Thus, the
Hessian matrix at optima is singular. Because the residual term $R$ is small
near the optima by continuity, the Hessian matrix is almost singular in a
neighborhood of the optima. This behavior is demonstrated by the numerical
result displayed in Fig. 1 (a). On the other hand, when the symmetry
constraint is imposed, Fig. 1 (b) shows that the Hessian matrix is positive
definite in a neighborhood of the maximal solution and it agrees our
theoretical results about the local strong convexity of the Hessian matrix.

Figure 1: The smallest singular value ($\sigma_{\text{min}}$) or eigenvalue
($\lambda_{\text{min}}$) of the Hessian matrix evaluated at 100 randomly
sampled $\Phi$ near the optimizer $\Phi^{*}$ of a given polynomial. Blue dots
correspond to an odd polynomial of degree $61$ and orange crosses correspond
to an even polynomial of degree $80$. (a) The symmetry constraint is not
imposed. The plotted data concentrates around the machine precision $\sim
10^{-17}$, which implies a large basin near the optimizer on which the Hessian
matrix is singular if the symmetry constraint is not imposed. (b) The symmetry
constraint is imposed. The numerical results show that the Hessian matrix is
well conditioned near the optimizer, and it is positive definite which agrees
our result of local strong convexity.

###  7.2 Existence of local minima

To visualize the full energy landscape of the optimization problem, we first
consider the simplest scenario so that only two phase factors suffice to solve
the optimization problem exactly. We choose two target polynomials
$f(x)=x^{2}-\frac{1}{2}$ and
$f(x)=\frac{1}{\sqrt{3}}x^{3}-\frac{2}{\sqrt{3}}x$. The energy landscape on
the irreducible domain is displayed in Fig. 2. The global minima derived from
the proposed method are annotated in the domain. In this special case, there
is not any local minimum in the landscape. However, this is a rare exception
rather than the norm. Existence of local minima can be observed by slightly
increasing the degree of the target polynomial
$f(x)=0.2103T_{4}(x)+0.1231T_{2}(x)+0.1666$. The local minimizer is
numerically searched by randomly initiating the stochastic gradient descent
algorithm to solve Eq. 5. The value of the objective function at
$\Phi^{\text{loc}}$ is $0.0218$. The Hessian matrix has eigenvalues
$0.1359,4.5815,7.7510$, which indicates that $\Phi^{\text{loc}}$ is indeed a
local minimum. To visualize the landscape, we plot a two-dimensional section
of the optimization landscape spanned by the affine plane spanned by the
eigenvectors of the least two eigenvalues of the Hessian matrix (Fig. 2).

Figure 2: The optimization landscape of the modified cost function
$F(\widetilde{\Phi})^{1/3}$ on the irreducible domain $D_{d}$. Here, the cube
root is taken to signify the structure of the landscape near optima. The
annotated values are all optima exactly computed from the proposed method. (a)
The target function is set to $f(x)=x^{2}-\frac{1}{2}$ which yields 8
inequivalent optima. (b) The target function is set to
$f(x)=\frac{1}{\sqrt{3}}x^{3}-\frac{2}{\sqrt{3}}x$ which yields 4 inequivalent
optima. Figure 3: A two-dimensional section of the optimization landscape
$F(\widetilde{\Phi}^{\mathrm{loc}}+xv_{1}+yv_{2})$, where
$\Phi^{\mathrm{loc}}$ is a local minimum of a given degree-$5$ target
polynomial, and $v_{1},v_{2}$ are unit eigenvectors corresponding to the least
two eigenvalues of the Hessian matrix evaluated at $\Phi^{\mathrm{loc}}$. The
objective value of the local minimum is annotated on the figure.

###  7.3 Numerical behavior near the maximal solution

Here we demonstrate that among all global minima, the maximal solution is
particularly desirable from the perspective of numerical optimization. Given
an arbitrarily chosen polynomial
$f(x)=\frac{1}{4}T_{6}(x)+\frac{5}{4}T_{4}(x)+\frac{1}{8}T_{2}(x)-T_{0}(x)$,
we consider a sequence of scaled polynomials $f^{(k)}(x):=10^{-k}f(x)$, so
that $\lim_{k\to\infty}f^{(k)}(x)=0$. All exact solutions to the optimization
problem associated to $f^{(k)}(x)$ can be constructed explicitly via Remark
16. This procedure yields several sequences of phase factors, which converge
to different limiting points as $k\to\infty$. For the specific degree-$6$
polynomial above, there are four distinct limits which is referred to as
$\Phi_{0}^{\text{class}}$,

  * •

class 1: $\Phi^{1}_{0}=\left(\frac{\pi}{4},0,0,0,0,0,\frac{\pi}{4}\right),$

  * •

class 2:
$\Phi^{2}_{0}=\left(\frac{\pi}{4},0,\frac{\pi}{4},-\frac{\pi}{2},\frac{\pi}{4},0,\frac{\pi}{4}\right),$

  * •

class 3:
$\Phi^{3}_{0}=\left(\frac{\pi}{4},0,-\frac{\pi}{4},\frac{\pi}{2},-\frac{\pi}{4},0,\frac{\pi}{4}\right),$

  * •

class 4:
$\Phi^{4}_{0}=\left(\frac{\pi}{4},\frac{\pi}{4},0,-\frac{\pi}{2},0,\frac{\pi}{4},\frac{\pi}{4}\right).$

Specifically, the first class is associated with the maximal solution. The
sequence of the set of phase factors corresponding to the scaled polynomial
$f^{(k)}(x)$ and a given class of construction is referred to as
$\Phi^{\text{class}}_{*}(k)$. Fig. 4 shows that the class associated with the
maximal solution distinguishes from other classes, in the sense that the
convergence rate of $\Phi^{\text{1}}_{*}(k)$ towards $\Phi^{\text{1}}_{0}$ is
much faster. This suggests that the convergence basin near
$\Phi^{\text{1}}_{0}$ is much flatter, which justifies the choice of using
$\Phi^{\text{1}}_{0}$ as the initial guess of the optimization.

To further test the performance of numerical optimization associated with
different initial guesses, we use the limit point $\Phi_{0}^{\text{class}}$ as
the initial guess and find the optimal phase factors
$\Phi_{*}^{\text{class}}(k)$ generating the scaled polynomial $f^{(k)}(x)$ by
running QSPPACK [5]. The trajectory of the optimization process is displayed
in Fig. 5. It shows that the optimization starting from $\Phi^{\text{1}}_{0}$
also has the fastest convergence rate.

Figure 4: The convergence to the limit point of different class as the target
polynomial being scaled down. Given a polynomial $f(x)$, the optimal set of
phase factors $\Phi_{*}^{\text{class}}(k)$ parameterizes the scaled polynomial
$f^{(k)}(x):=10^{-k}f(x)$. Here, $\Phi^{\mathrm{class}}_{0}$ is the limit
point of the corresponding class. Figure 5: The objective value in each
iteration step by running QSPPACK. Initiating the optimization from the limit
point of a class $\Phi_{0}^{\text{class}}$, $\Phi_{t}^{\text{class}}$ is the
set of phase factors in the $t$-th optimization step. Different symbols
correspond to the optimization trajectories starting from distinct classes of
initial guess.

## References

  * [1] Dimitri P Bertsekas.  On the Goldstein-Levitin-Polyak gradient projection method.  IEEE Transactions on automatic control, 21(2):174–184, 1976.
  * [2] Sébastien Bubeck.  Convex optimization: Algorithms and complexity, 2015.
  * [3] Rui Chao, Dawei Ding, Andras Gilyen, Cupjin Huang, and Mario Szegedy.  Finding angles for quantum signal processing with machine precision, 2020\.
  * [4] Andrew M. Childs, Dmitri Maslov, Yunseong Nam, Neil J. Ross, and Yuan Su.  Toward the first quantum simulation with quantum speedup.  Proc. Nat. Acad. Sci., 115:9456–9461, 2018.
  * [5] Yulong Dong, Xiang Meng, K. Birgitta Whaley, and Lin Lin.  Efficient phase factor evaluation in quantum signal processing, 2020.
  * [6] András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe.  Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics.  In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, pages 193–204. ACM, 2019.
  * [7] Gene H. Golub and Charles F. Van Loan.  Matrix Computations.  The Johns Hopkins University Press, third edition, 1996.
  * [8] Jeongwan Haah.  Product decomposition of periodic functions in quantum signal processing.  Quantum, 3:190, 2019.
  * [9] Nicholas J. Higham.  Accuracy and Stability of Numerical Algorithms.  Society for Industrial and Applied Mathematics, second edition, 2002.
  * [10] J. L. W. V. Jensen.  Sur un nouvel et important théorème de la théorie des fonctions.  Acta Mathematica, 22(none):359 – 364, 1900.
  * [11] C. T. Kelley.  Iterative methods for optimization, volume 18.  SIAM, 1999.
  * [12] Lin Lin and Yu Tong.  Near-optimal ground state preparation.  Quantum, 4:372, 2020.
  * [13] Lin Lin and Yu Tong.  Optimal quantum eigenstate filtering with application to solving quantum linear systems.  Quantum, 4:361, 2020.
  * [14] Guang Hao Low and Isaac L Chuang.  Optimal hamiltonian simulation by quantum signal processing.  Physical review letters, 118(1):010501, 2017.
  * [15] Kurt Mahler.  On some inequalities for polynomials in several variables.  Journal of The London Mathematical Society-second Series, pages 341–344, 1962.
  * [16] John M Martyn, Zane M Rossi, Andrew K Tan, and Isaac L Chuang.  A grand unification of quantum algorithms.  arXiv:2105.02859, 2021.
  * [17] Michael A Nielsen and Isaac Chuang.  Quantum computation and quantum information, 2000.
  * [18] J. Nocedal and S. J. Wright.  Numerical optimization.  Springer Verlag, 1999.

##  Appendix A Proof of Lemma 10

###### Proof.

One can check that

$\begin{split}&P(x)=e^{\mathrm{i}\left(\phi_{0}+\phi_{d}\right)}\bra{0}W(x)\prod_{j=1}^{d-1}\left[e^{\mathrm{i}\phi_{j}Z}W(x)\right]\ket{0},\\\
&\mathrm{i}\sqrt{1-x^{2}}Q(x)=e^{\mathrm{i}\left(\phi_{0}-\phi_{d}\right)}\bra{0}W(x)\prod_{j=1}^{d-1}\left[e^{\mathrm{i}\phi_{j}Z}W(x)\right]\ket{1}.\end{split}$
(96)

Then, we expand the product of matrices from left to right by using
$e^{\mathrm{i}\phi Z}=\cos(\phi)+\mathrm{i}\sin(\phi)Z$ and
$W(x)Z=ZW(x)^{-1}$. We want to show that the term containing $\sin(\phi_{j})$
does not contribute to $p_{d}$ and $q_{d-1}$. We use $l.o.$ to denote any
polynomial or polynomial-valued matrix of degree $<d$. Here, we relax the
definition of polynomial and classify the entry of the matrix of our interest
as polynomial by considering $\sqrt{1-x^{2}}$ as another variable. Note that

$\begin{split}&W(x)\prod_{j=1}^{d-1}\left[e^{\mathrm{i}\phi_{j}Z}W(x)\right]\\\
&=\cos\left(\phi_{1}\right)W(x)^{2}\prod_{j=2}^{d-1}\left[e^{\mathrm{i}\phi_{j}Z}W(x)\right]+\mathrm{i}\sin\left(\phi_{1}\right)W(x)ZW(x)\prod_{j=2}^{d-1}\left[e^{\mathrm{i}\phi_{j}Z}W(x)\right]\\\
&=\cos\left(\phi_{1}\right)W(x)^{2}\prod_{j=2}^{d-1}\left[e^{\mathrm{i}\phi_{j}Z}W(x)\right]+\mathrm{i}\sin\left(\phi_{1}\right)Z\prod_{j=2}^{d-1}\left[e^{\mathrm{i}\phi_{j}Z}W(x)\right]\\\
&=\cos\left(\phi_{1}\right)W(x)^{2}\prod_{j=2}^{d-1}\left[e^{\mathrm{i}\phi_{j}Z}W(x)\right]+l.o.,\end{split}$
(97)

where the product of the second term is indeed a matrix whose elements are
polynomials of degree $\leq d-2$.

Using the fact that

$[W(x)]^{d}=\begin{pmatrix}T_{d}(x)&\mathrm{i}\sqrt{1-x^{2}}U_{d-1}(x)\\\
\mathrm{i}\sqrt{1-x^{2}}U_{d-1}(x)&T_{d}(x)\end{pmatrix},$ (98)

we may prove inductively,

$\begin{split}&P(x)=e^{\mathrm{i}\left(\phi_{0}+\phi_{d}\right)}\prod_{j=1}^{d-1}\cos\left(\phi_{j}\right)T_{d}(x)+l.o.\\\
&\mathrm{i}\sqrt{1-x^{2}}Q(x)=\mathrm{i}\sqrt{1-x^{2}}e^{\mathrm{i}\left(\phi_{0}-\phi_{d}\right)}\prod_{j=1}^{d-1}\cos\left(\phi_{j}\right)U_{k-1}(x)+l.o.,\end{split}$
(99)

which proves the lemma. ∎

##  Appendix B Proof of Lemma 11

###### Proof.

Now that $d$ is odd and $\Phi$ is symmetric, Lemma 10 indicates

$q_{d-1}=\prod_{j=1}^{\widetilde{d}-1}\cos^{2}\left(\phi_{j}\right)\geq 0.$

If $\left\lvert\phi_{j}\right\rvert\neq\frac{\pi}{2}$ for all $1\leq
j\leq\widetilde{d}-1$, then the leading Chebyshev coefficient of $Q$ is
$q_{d-1}$, which is positive. Otherwise, there exists positive integer $j_{0}$
such that $\phi_{j_{0}}=\phi_{d-j_{0}}=\pm\frac{\pi}{2}$.

If $j_{0}<\widetilde{d}-1$, using the identity

$W(x)e^{\pm\mathrm{i}\frac{\pi}{2}Z}W(x)=e^{\pm\mathrm{i}\frac{\pi}{2}Z},$

the product of matrices

$e^{\mathrm{i}\phi_{j_{0}-1}Z}W(x)e^{\mathrm{i}\phi_{j_{0}}Z}W(x)e^{\mathrm{i}\phi_{j_{0}+1}Z}$

can be replaced by

$e^{\mathrm{i}\left(\phi_{j_{0}-1}+\phi_{j_{0}}+\phi_{j_{0}+1}\right)Z}.$

We can perform a similar reduction for for

$e^{\mathrm{i}\phi_{d-j_{0}-1}Z}W(x)e^{\mathrm{i}\phi_{d-j_{0}}Z}W(x)e^{\mathrm{i}\phi_{d-j_{0}+1}Z}.$

Therefore, $U(x,\Phi)$ can be constructed by a new set of symmetric phase
factors whose full length is $d-3$ and $\deg(Q)\leq d-5$. Apply Lemma 10 again
and one has $q_{d-5}\geq 0$.

If $j_{0}=\widetilde{d}-1$, the product of matrices

$e^{\mathrm{i}\phi_{\widetilde{d}-2}Z}W(x)e^{\mathrm{i}\phi_{\widetilde{d}-1}Z}W(x)e^{\mathrm{i}\phi_{\widetilde{d}}Z}W(x)e^{\mathrm{i}\phi_{\widetilde{d}+1}Z}$

can be replaced by

$-e^{\mathrm{i}\phi_{\widetilde{d}-2}Z}W(x)e^{\mathrm{i}\phi_{\widetilde{d}+1}Z}.$

Therefore, $U(x,\Phi)$ can be constructed by a new set of symmetric phase
factors whose full length is $d-1$ and $\deg(Q)\leq d-3$. Apply Lemma 10 again
and one has $q_{d-3}\leq 0$.

Repeat the procedure above until the length of the new set of symmetric phase
factors reaches $\text{deg}(Q)+2$, and it follows that the leading Chebyshev
coefficient of $Q$ has the same sign as
$\left(-1\right)^{\frac{d-1-\deg(Q)}{2}}$. ∎

##  Appendix C Proof of Lemma 12

First, we consider the expansion of $P(x)$ and $Q(x)$ in the monomial basis,