Gradient Descent Learns Linear Dynamical Systems
Moritz Hardt∗

Tengyu Ma†

Benjamin Recht‡

arXiv:1609.05191v1 [cs.LG] 16 Sep 2016

September 19, 2016

Abstract
We prove that gradient descent efficiently converges to the global optimizer of the maximum
likelihood objective of an unknown linear time-invariant dynamical system from a sequence of
noisy observations generated by the system. Even though the objective function is non-convex,
we provide polynomial running time and sample complexity bounds under strong but natural
assumptions. Linear systems identification has been studied for many decades, yet, to the best
of our knowledge, these are the first polynomial guarantees for the problem we consider.

1

Introduction

Many learning problems are by their nature sequence problems where the goal is to fit a model that
maps a sequence of input words x1 , . . . , xT to a corresponding sequence of observations y1 , . . . , yT .
Text translation, speech recognition, time series prediction, video captioning and question answering
systems, to name a few, are all sequence to sequence learning problems. For a sequence model to
be both expressive and parsimonious in its parameterization, it is crucial to equip the model with
memory thus allowing its prediction at time t to depend on previously seen inputs.
Recurrent neural networks form an expressive class of non-linear sequence models. Through
their many variants, such as long-short-term-memory [HS97], recurrent neural networks have seen
remarkable empirical success in a broad range of domains. At the core, neural networks are typically
trained using some form of (stochastic) gradient descent. Even though the training objective is
non-convex, it is widely observed in practice that gradient descent quickly approaches a good set of
model parameters. Understanding the effectiveness of gradient descent for non-convex objectives
on theoretical grounds is a major open problem in this area.
If we remove all non-linear state transitions from a recurrent neural network, we are left with the
state transition representation of a linear dynamical system. Notwithstanding, the natural training
objective for linear systems remains non-convex due to the composition of multiple linear operators
in the system. If there is any hope of eventually understanding recurrent neural networks, it will
be inevitable to develop a solid understanding of this special case first.
To be sure, linear dynamical systems are important in their own right and have been studied
for many decades independently of machine learning within the control theory community. Control
theory provides a rich set techniques for identifying and manipulating linear systems. The learning problem in this context corresponds to “linear dynamical system identification”. Maximum
∗

Google. Email: m@mrtz.org
Princeton University. Email: tengyu@cs.princeton.edu. Work performed at Google. Partially supported by
Simons Award in Theoretical Computer Science, IBM PhD Fellowship, Dodds Fellowship and Siebel Scholarship.
‡
UC Berkeley. Email: brecht@berkeley.edu. Work performed at Google.
†

1

likelihood estimation with gradient descent is a popular heuristic for dynamical system identification [Lju98]. In the context of machine learning, linear systems play an important role in
numerous tasks. For example, their estimation arises as subroutines of reinforcement learning in
robotics [LK13], location and mapping estimation in robotic systems [DWB06], and estimation of
pose from video [RRD05].
In this work, we show that gradient descent efficiently minimizes the maximum likelihood objective of an unknown linear system given noisy observations generated by the system. More formally,
we receive noisy observations generated by the following time-invariant linear system:
ht+1 = Aht + Bxt

(1.1)

yt = Cht + Dxt + ξt
Here, A, B, C, D are linear transformations with compatible dimensions and we denote by Θ =
(A, B, C, D) the parameters of the system. The vector ht represents the hidden state of the model
at time t. Its dimension n is called the order of the system. The stochastic noise variables {ξt }
perturb the output of the system which is why the model is called an output error model in control
theory. We assume the variables are drawn i.i.d. from a distribution with mean 0 and variance σ 2 .
Throughout the paper we focus on controllable and externally stable systems. A linear system
is externally stable (or equivalently bounded-input bounded-output stable) if and only if the spectral
radius of A, denoted ρ(A), is strictly bounded by 1. Controllability is a mild non-degeneracy
assumption that we formally define later. Without loss of generality, we further assume that the
transformations B, C and D have bounded Frobenius norm. This can be achieved by a rescaling
of the output variables. We assume we have N pairs of sequences (x, y) as training examples,
n
o
S = (x(1) , y (1) ), . . . , (x(N ) , y (N ) ) .
Each input sequence x ∈ RT of length T is drawn from a distribution and y is the corresponding
output of the system above generated from an unknown initial state h. We allow the unknown
initial state to vary from one input sequence to the next. This only makes the learning problem
more challenging.
Our goal is to fit a linear system to the observations. We parameterize our model in exactly
the same way as (1.1). That is, for linear mappings (Â, B̂, Ĉ, D̂), the trained model is defined as:
ĥt+1 = Âĥt + B̂xt ,

ŷt = Ĉ ĥt + D̂xt

(1.2)

The (population) risk of the model is obtained by feeding the learned system with the correct
initial states and comparing its predictions with the ground truth in expectation over inputs and
errors. Denoting by ŷt the t-th prediction of the trained model starting from the correct initial state
b as a short hand for (Â, B̂, Ĉ, D̂), we formally define population risk
that generated yt , and using Θ
as:
"
#
T
X
1
2
b =
f (Θ)
kŷt − yt k
(1.3)
E
{xt },{ξt } T
t=1

Note that even though the prediction ŷt is generated from the correct initial state, the learning
algorithm does not have access to the correct initial state for its training sequences.
While the squared loss objective turns out to be non-convex, it has many appealing properties.
Assuming the inputs xt and errors ξt are drawn independently from a Gaussian distribution, the
corresponding population objective corresponds to maximum likelihood estimation. In this work,
we make the weaker assumption that the inputs and errors are drawn independently from possibly
2

different distributions. The independence assumption is certainly idealized for some learning applications. However, in control applications the inputs can often be chosen by the controller rather
than by nature. Moreover, the outputs of the system at various time steps are correlated through
the unknown hidden state and therefore not independent even if the inputs are.

1.1

Our results

We show that we can efficiently minimize the population risk using projected stochastic gradient
descent. The bulk of our work applies to single-input single-output (SISO) systems meaning that
inputs and outputs are scalars xt , yt ∈ R. However, the hidden state can have arbitrary dimension n.
Every controllable SISO admits a convenient canonical form called controllable canonical form that
we formally introduce later. In this canonical form, the transition matrix A is governed by n
parameters a1 , . . . , an which coincide with the coefficients of the characteristic polynomial of A.
The minimal assumption under which we might hope to learn the system is that the spectral radius
of A is smaller than 1. However, the set of such matrices is non-convex and does not have enough
structure for our analysis. We will therefore make additional assumptions. The assumptions we
need differ between the case where we are trying to learn A with n parameter system, and the case
where we allow ourselves to over-specify the trained model with n0 > n parameters. The former
is sometimes called proper learning, while the latter is called improper learning. In the improper
case, we are essentially able to learn any system with spectral radius less than 1 under a mild
separation condition on the roots of the characteristic polynomial. Our assumption in the proper
case is stronger and we introduce it next.

1.2

Proper learning

Suppose that the state transition matrix A is given by parameters a1 , . . . , an and consider the
polynomial q(z) = 1 + a1 z + a2 z 2 + · · · + an z n over the complex numbers C.
We will require that the image of the unit circle on the
complex plane under the polynomial q is contained in the
Complex plane
cone of complex numbers whose real part is larger than
their absolute imaginary part. Formally, for all z ∈ C
1
such that |z| = 1, we require that <(q(z)) > |=(q(z))|.
Here, <(z) and =(z) denote the real and imaginary part
of z, respectively. We illustrate this condition in the figure
on the right for a degree 4 system.
0
Our assumption has three important implications.
First, it implies (via Rouché’s theorem) that the spectral radius of A is smaller than 1 and therefore ensures
1
the stability of the system. Second, the vectors satisfying
n
our assumption form a convex set in R . Finally, it en1
0
1
sures that the objective function is weakly quasi-convex, a
condition we introduce later when we show that it enables
stochastic gradient descent to make sufficient progress.
√
We note in passing that our assumption can be satisfied via the `1 -norm constraint kak1 ≤
√ 2/2.
Moreover, if we pick random Gaussian coefficients with expected norm bounded by o(1/ log n),
then the resulting vector will satisfy our assumption with probability 1 − o(1). Roughly speaking,
the assumption requires the roots of the characteristic polynomial p(z) = z n + a1 z n−1 + · · · + an
are relatively dispersed inside the unit circle. (For comparison, on the other end of the spectrum,

3

the polynomial p(z) = (z − 0.99)n have all its roots colliding at the same point and doesn’t satisfy
the assumption.)
Theorem 1.1 (Informal). Under our assumption, projected stochastic gradient descent, when
b with population risk
given N sample sequence of length T , returns parameters Θ
!
r
5 + σ 2 n3
n
b ≤ f (Θ) + O
.
f (Θ)
TN
Recall that f (Θ) is the population risk of the optimal system, and σ 2 refers to the variance
of the noise variables. We also assume that the inputs xt are drawn from a pairwise independent
distribution with mean 0 and variance 1. Note, however, that this does not imply independence
of the outputs as these are correlated by a common hidden state. The stated version of our
result glosses over the fact that we need our assumption to hold with a small amount of slack;
a precise version follows in Section 4. Our theorem establishes a polynomial convergence rate for
stochastic gradient descent. Since each iteration of the algorithm only requires a sequence of matrix
operations and an efficient projection step, the running time is polynomial, as well. Likewise, the
sample requirements are polynomial since each iteration requires only a single fresh example. An
important feature of this result is that the error decreases with both the length T and the number
of samples N . Computationally, the projection step can be a bottleneck, although it is unlikely to
be required in practice and may be an artifact of our analysis.

1.3

The power of over-parameterization

Endowing the model with additional parameters compared to the ground truth turns out to be
surprisingly powerful. We show that we can essentially remove the assumption we previously made
in proper learning. The idea is simple. If p is the characteristic polynomial of A of degree n. We
can find a system of order n0 > n such that the characteristic polynomial of its transition matrix
becomes p · p0 for some polynomial p0 of order n0 − n. This means that to apply our result we only
need the polynomial p · p0 to satisfy our assumption. At this point, we can choose p0 to be an
approximation of the inverse p−1 . For sufficiently good approximation, the resulting polynomial
p · p0 is close to 1 and therefore satisfies our assumption. Such an approximation exists generically
for n0 = O(n) under mild non-degeneracy assumptions on the roots of p. In particular, any small
random perturbation of the roots would suffice.
Theorem 1.2 (Informal). Under a mild non-degeneracy assumption, stochastic gradient descent
b corresponding to a system of order n0 = O(n) with population risk
returns parameters Θ
!
r
5 + σ 2 n3
n
b ≤ f (Θ) + O
f (Θ)
,
TN
when given N sample sequences of length T .
We remark that the idea we sketched also shows that, in the extreme, improper learning of
linear dynamical systems becomes easy in the sense that the problem can be solved using linear
regression against the outputs of the system. However, our general result interpolates between the
proper case and the regime where linear regression works. We discuss in more details in Section 6.3.

4

1.4

Multi-input multi-output systems

Both results we saw immediately extend to single-input multi-output (SIMO) systems as the dimensionality of C and D are irrelevant for us. The case of multi-input multi-output (MIMO) systems
is more delicate. Specifically, our results carry over to a broad family of multi-input multi-output
systems. However, in general MIMO systems no longer enjoy canonical forms like SISO systems. In
Section 7, we introduce a natural generalization of controllable canonical form for MIMO systems
and extend our results to this case.

1.5

Related work

System identification is a core problem in dynamical systems and has been studied in depth for
many years. The most popular reference on this topic is the text by Ljung [Lju98]. Nonetheless,
the list of non-asymptotic results on identifying linear systems from noisy data is surprisingly
short. Several authors have recently tried to estimate the sample complexity of dynamical system
identification using machine learning tools [VK08, CW02, WC99]. All of these result are rather
pessimistic with sample complexity bounds that are exponential in the degree of the linear system
and other relevant quantities. Contrastingly, we prove that gradient descent has an associated
polynomial sample complexity in all of these parameters. Moreover, all of these papers only focus
on how well empirical risk approximates the true population risk and do not provide guarantees
about any algorithmic schemes for minimizing the empirical risk.
The only result to our knowledge which provides polynomial sample complexity for identifying
linear dynamical systems is in Shah et al [SBTR12]. Here, the authors show that if certain frequency
domain information about the linear dynamical system is observed, then the linear system can be
identified by solving a second-order cone programming problem. This result is about improper
learning only, and the size of the resulting system may be quite large, scaling as the (1 − ρ(A))−2 .
As we describe in this work, very simple algorithms work in the improper setting when the system
degree is allowed to be polynomial in (1 − ρ(A))−1 . Moreover, it is not immediately clear how to
translate the frequency-domain results to the time-domain identification problem discussed above.
Our main assumption about the image of the polynomial q(z) is an appeal to the theory of
passive systems. A system is passive if the dot product between the input sequence ut and output
sequence yt are strictly positive. Physically, this notion corresponds to systems that cannot create
energy. For example, a circuit made solely of resistors, capacitors, and inductors would be a
passive electrical system. If one added an amplifier to the internals of the system, then it would
no longer be passive. The set of passive systems is a subset of the set of stable systems, and the
subclass is somewhat easier to work with mathematically. Indeed, Megretski used tools from passive
systems to provide a relaxation technique for a family of identification problems in dynamical
systems [Meg08]. His approach is to lower bound a nonlinear least-squares cost with a convex
functional. However, he does not prove that his technique can identify any of the systems, even
asymptotically. Bazanella et al use a passivity condition to prove quasi-convexity of a cost function
that arises in control design [BGMA08]. Building on this work, Eckhard and Bazanella prove
a weaker version of Lemma 3.3 in the context of system identification [EB11], but no sample
complexity or global convergence proofs are provided.

1.6

Proof overview

The first important step in our proof is to develop population risk in Fourier domain where it is
closely related to what we call idealized risk. Idealized risk essentially captures the `2 -difference
between the transfer function of the learned system and the ground truth. The transfer function
5

is a fundamental object in control theory. Any linear system is completely characterized by its
transfer function G(z) = C(zI − A)−1 B. In the case of a SISO, the transfer function is a rational
function of degree n over the complex numbers and can be written as G(z) = s(z)/p(z). In the
canonical form introduced in Section 1.7, the coefficients of p(z) are precisely the parameters that
specify A. Moreover, z n p(1/z) = 1 + a1 z + a2 z 2 + · · · + an z n is the polynomial we encountered in the
introduction. Under the assumption illustrated earlier, we show in Section 3 that the idealized risk
is weakly quasi-convex (Lemma 3.3). Quasi-convexity implies that gradients cannot vanish except
at the optimum of the objective function; we review this (mostly known) material in Section 2. In
particular, this lemma implies that in principle we can hope to show that gradient descent converges
to a global optimum. However, there are several important issues that we need to address. First,
the result only applies to idealized risk, not our actual population risk objective. Therefore it is not
clear how to obtain unbiased gradients of the idealized risk objective. Second, there is a subtlety in
even defining a suitable empirical risk objective. The reason is that risk is defined with respect to
the correct initial state of the system which we do not have access to during training. We overcome
both of these problems. In particular, we design an almost unbiased estimator of the gradient of
the idealized risk in Lemma 5.4 and give variance bounds of the gradient estimator (Lemma 5.5).
Our results on improper learning in Section 6 rely on a surprisingly simple but powerful insight. We can extend the degree of the transfer function G(z) by extending both numerator and
denominator with a polynomial u(z) such that G(z) = s(z)u(z)/p(z)u(z). While this results in an
equivalent system in terms of input-output behavior, it can dramatically change the geometry of
the optimization landscape. In particular, we can see that only p(z)u(z) has to satisfy the assumption of our proper learning algorithm. This allows us, for example, to put u(z) ≈ p(z)−1 so that
p(z)u(z) ≈ 1, hence trivially satisfying our assumption. A suitable inverse approximation exists
under light assumptions and requires degree no more than d = O(n). Algorithmically, there is
almost no change. We simply run stochastic gradient descent with n + d model parameters rather
than n model parameters.

1.7

Preliminaries

For complex matrix (or vector, number) C, we use <(C) to denote the real part and =(C) the
imaginary part, and C̄ the conjugate and C ∗ = C̄ > its conjugate transpose . We use | · | to denote
the absolute value of a complex number
c. For complex vector u and v, we use hu, vi = u∗ v to
√
denote the inner product and kuk = u∗ u is the norm of u. For complexpmatrix A and B with same
dimension, hA, Bi = tr(A∗ B) defines an inner product, and kAkF = tr(A∗ A) is the Frobenius
norm. For a square matrix A, we use ρ(A) to denote the spectral radius of A, that is, the largest
absolute value of the elements in the spectrum of A. We use In to denote the identity matrix with
dimension n × n, and we drop the subscript when it’s clear from the context.We let ei denote the
i-th standard basis vector.
A SISO of order n is in controllable canonical form if A and B have the following form


 
0
1
0
···
0
0
 0
 0 

0
1
·
·
·
0


 

 .. 
..
..
.. 
..
A =  ...
B
=
(1.4)

 . 
.
.
.
. 

 
 0
 0 
0
0
···
1 
1
−an −an−1 −an−2 · · · −a1
We will parameterize Â, B̂, Ĉ, D̂ accordingly. We will write A = CC(a) for brevity, where a
is used to denote the unknown last row [−an , . . . , −a1 ] of matrix A. We will use â to denote
6

the corresponding training variables for a. Since here B is known, so B̂ is no longer a trainable
parameter, and is forced to be equal to B. Moreover, C is a row vector and we use [c1 , · · · , cn ] to
denote its coordinates (and similarly for Ĉ).
A SISO is controllable if and only if the matrix [B | AB | A2 B | · · · | An−1 B] has rank n.
This statement corresponds to the condition that all hidden states should be reachable from some
initial condition and input trajectory. Any controllable system admits a controllable canonical
form [HRS07]. For vector a = [an , . . . , a1 ], let pa (z) denote the polynomial
pa (z) = z n + a1 z n−1 + · · · + an .

(1.5)

When a defines the matrix A that appears in controllable canonical form, then pa is precisely the
characteristic polynomial of A. That is, pa (z) = det(zI − A).

2

Gradient descent and quasi-convexity

It is known that under certain mild conditions (stochastic) gradient descent converges even on nonconvex functions to local minimum [GHJY15, LSJR16]. Though usually for concrete problems the
challenge is to prove that there is no spurious local minimum other than the target solution. Here
we introduce a condition similar to the quasi-convexity notion in [HLS15], which ensures that any
point with vanishing gradient is the optimal solution . Roughly speaking, the condition says that at
any point θ the negative of the gradient −∇f (θ) should be positively correlated with direction θ∗ −θ
pointing towards the optimum. Our condition is slightly weaker than that in [HLS15] since we only
require quasi-convexity and smoothness with respect to the optimum, and this (simple) extension
will be necessary for our analysis.
Definition 2.1 (Weak quasi-convexity). We say an objective function f is τ -weakly-quasi-convex
(τ -WQC) over a domain B with respect to global minimum θ∗ if there is a positive constant τ > 0
such that for all θ ∈ B,
∇f (θ)> (θ − θ∗ ) ≥ τ (f (θ) − f (θ∗ )) .
(2.1)
We further say f is Γ-weakly-smooth if for for any point θ, k∇f (θ)k2 ≤ Γ(f (θ) − f (θ∗ )).
Note that indeed any Γ-smooth convex function in the usual sense is O(Γ)-weakly-smooth. For
a random vector X ∈ Rn , we define it’s variance to be Var[X] = E[kX − EXk2 ].
Definition 2.2. We call r(θ) an unbiased estimator of ∇f (θ) with variance V if it satisfies E[r(θ)] =
∇f (θ) and Var[r(θ)] ≤ V .
Projected stochastic gradient descent over some closed convex set B with learning rate η > 0
refers to the following algorithm in which ΠB denotes the Euclidean projection onto B:
for k = 0 to K − 1 :
wk+1 = θk − ηr(θk )
θk+1 = ΠB (wk+1 )
return θj with j uniformly picked from 1, . . . , K

(2.2)

The following Proposition is well known for convex objective functions (corresponding to 1weakly-quasi-convex functions). We extend it (straightforwardly) to the case when τ -WQC holds
with any positive constant τ .

7

Proposition 2.3. Suppose the objective function f is τ -weakly-quasi-convex and Γ-weakly-smooth,
and r(·) is an unbiased estimator for ∇f (θ) with variance V . Moreover, suppose the global minimum
θ∗ belongs to B, and the initial point θ0 satisfies kθ0 −θ∗ k ≤ R. Then projected gradient descent (2.2)
with a proper learning rate returns θK in K iterations with expected error
(
√ )!
ΓR2 R V
f (θK ) ≤ O max
, √
.
τ 2K τ K
Remark 1. It’s straightforward to see (from the proof) that the algorithm tolerates inverse exponential
in√the

n bias
ogradient estimator. Technically, suppose E[r(θ)] = ∇f (θ) ± ζ then f (θK ) ≤
ΓR2 R√ V
O max τ 2 K , τ K
+ poly(K) · ζ. Throughout the paper, we assume that the error that we are
shooting for is inverse polynomial and therefore the effect of inverse exponential bias is negligible.
We defer the proof of Proposition 2.4 to Appendix A which is a simple variation of the standard
convergence analysis of stochastic gradient descent (see, for example, [Bot98]). Finally, we note
that the sum of two quasi-convex functions may no longer be quasi-convex. However, if a sequence
functions is τ -WQC with respect to a common point θ∗ , then their sum is also τ -WQC. This follows
from the linearity of gradient operation.
Proposition 2.4. Suppose functions f1 , . . . , fn are individually τ -weakly-quasi-convex in B with
∗
respect
w1 , . . . , wn the linear combination
Pnto a common global minimum θ , then for non-negative
f = i=1 wi fi is also τ -weakly-quasi-convex with respect to θ∗ in B.

3

Population risk in frequency domain

We next establish conditions under which risk is weakly-quasi-convex. Our strategy is to first
b by what we call idealized risk. This approximation of the
approximate the risk functional f (Θ)
objective function is fairly straightforward; we justify it toward the end of the section. We can
show that

b ≈ kD − D̂k2 + P∞ Ĉ Âk B − CAk B 2 .
(3.1)
f (Θ)
k=0
The leading term kD − D̂k2 is convex in D̂ which appears nowhere else in the objective. It
therefore doesn’t affect the convergence of the algorithm (up to lower order terms) by virtue of
Proposition 2.4, and we restrict our attention to the remaining terms.
Definition 3.1 (Idealized risk). We define the idealized risk as
g(Â, Ĉ) =

∞ 
X

Ĉ Âk B − CAk B

2

.

(3.2)

k=0

We now use basic concepts from control theory (see [HRS07, Hes09] for more background) to
express the idealized risk (3.2) in Fourier domain. The transfer function of the linear system is
G(z) = C(zI − A)−1 B .

(3.3)

Note that G(z) is a rational function over the complex numbers of degree n and hence we
s(z)
can find polynomials s(z) and p(z) such that G(z) = p(z)
, with the convention that the leading

8

coefficient of p(z) is 1. In controllable canonical form (1.4), the coefficients of p will correspond to
the last row of the A, while the coefficients of s correspond to the entries of C. Also note that
G(z) =

∞
X

z −t CAt−1 B =

t=1

∞
X

z −t rt−1

t=1

The sequence r = (r0 , r1 , . . . , rt , . . .) = (CB, CAB, . . . , CAt B, . . .) is called the impulse response of
the linear system. The behavior of a linear system is uniquely determined by the impulse response
and therefore by G(z). Analogously, we denote the transfer function of the learned system by
b
G(z)
= Ĉ(zI − Â)−1 B = ŝ(z)/p̂(z) . The idealized risk (3.2) is only a function of the impulse
b
response r̂ of the learned system, and therefore it is only a function of G(z).
For future reference, we note that by some elementary calculation (see Lemma B.1), we have
G(z) = C(zI − A)−1 B =

zn

c1 + · · · + cn z n−1
,
+ a1 z n−1 + · · · + an

(3.4)

which implies that s(z) = c1 + · · · + cn z n−1 and p(z) = z n + a1 z n−1 + · · · + an .
With these definitions in mind, we are ready to express idealized risk in Fourier domain.
Proposition 3.2. Suppose pâ (z) has all its roots inside unit circle, then the idealized risk g(â, Ĉ)
can be written in the Fourier domain as
Z 2π
2
b iθ ) − G(eiθ ) dθ .
g(Â, Ĉ) =
G(e
0

b iθ ) the Fourier
Proof. Note that G(eiθ ) is the Fourier transform of the sequence {rk } and so is G(e
1
transform of r̂k . Therefore by Parseval’ Theorem, we have that
g(Â, Ĉ) =

∞
X

Z

2

kr̂k − rk k =

b iθ ) − G(eiθ )|2 dθ .
|G(e

0

k=0

3.1

2π

Quasi-convexity of the idealized risk

Now that we have a convenient expression for risk in Fourier domain, we can prove that the idealized
risk g(Â, Ĉ) is weakly-quasi-convex when â is not so far from the true a in the sense that pa (z) and
p̂a (z) have an angle less than π/2 for every z on the unit circle. We will use the convention that a
and â refer to the parameters that specify A and Â, respectively.
Lemma 3.3. For τ > 0 and every Ĉ, the idealized risk g(Â, Ĉ) is τ -weakly-quasi-convex over the
domain




pa (z)
n
Nτ (a) = â ∈ R : <
≥ τ /2, ∀ z ∈ C, s.t. |z| = 1 .
(3.5)
pâ (z)
b
b
Proof. We first analyze a single term h = |G(z)
− G(z)|2 . Recall that G(z)
= ŝ(z)/p̂(z) where
n
n−1
p̂(z) = pâ (z) = z + â1 z
+ · · · + ân . Note that z is fixed and h is a function of â and Ĉ. Then
it is straightforward to see that


 
∂h
1
ŝ(z) s(z) ∗
= 2<
−
.
(3.6)
∂ŝ(z)
p̂(z) p̂(z) p(z)
1

The Fourier transform exists since krk k2 = kĈ Âk B̂k2 ≤ kĈkkÂk kkB̂k ≤ cρ(Â)k where c doesn’t depend on k
and ρ(Â) < 1.

9

and
∂h
= −2<
∂ p̂(z)



ŝ(z)
p̂(z)2



ŝ(z) s(z)
−
p̂(z) p(z)

∗ 
.

(3.7)

Since ŝ(z) and p̂(z) are linear in Ĉ and â respectively, by chain rule we have that
h

∂h
∂h
∂h ∂ p̂(z)
∂h ∂ŝ(z)
, â − ai + h
, Ĉ − Ci =
h
, â − ai +
h
, Ĉ − Ci
∂â
∂ p̂(z) ∂â
∂ŝ(z) ∂ Ĉ
∂ Ĉ
∂h
∂h
=
(p̂(z) − p(z)) +
(ŝ(z) − s(z)) .
∂ p̂(z)
∂ŝ(z)

Plugging the formulas (3.6) and (3.7) for

∂h
∂ŝ(z)

and

∂h
∂ p̂(z)

into the equation above, we obtain that



−ŝ(z)(p̂(z) − p(z)) + p̂(z)(ŝ(z) − s(z))
p̂(z)2


 
ŝ(z)p(z) − s(z)p̂(z) ŝ(z) s(z) ∗
= 2<
−
p̂(z)2
p̂(z) p(z)
)
(
p(z) ŝ(z) s(z) 2
−
= 2<
p̂(z) p̂(z) p(z)


2
p(z)
b
G(z)
− G(z) .
= 2<
p̂(z)

∂h
∂h
h , â − ai + h
, Ĉ − Ci = 2<
∂â
∂ Ĉ



ŝ(z) s(z)
−
p̂(z) p(z)

∗ 

o
n
b
Hence h = |G(z)
− G(z)|2 is τ -weakly-quasi-convex with τ = 2 min|z|=γ < p(z)
p̂(z) . This implies
our claim, since by Proposition 3.2, the idealized risk g is convex combination of functions of the
b
form |G(z)
− G(z)|2 for |z| = 1. Moreover, Proposition 2.4 shows convex combination preserves
quasi-convexity.
For future reference, we also prove that the idealized risk is O(n2 /τ14 )-weakly smooth.
Lemma 3.4. The idealized risk g(Â, Ĉ) is Γ-weakly smooth with Γ = O(n2 /τ14 ).
Proof. By equation (3.7) and the chain rule we get that
∂g
∂ Ĉ

Z
=
T

b
∂|G(z)
− G(z)|2 ∂p(z)
·
dz =
∂p(z)
∂ Ĉ



Z
2<
T

1
p̂(z)



ŝ(z) s(z)
−
p̂(z) p(z)

∗ 

· [1, . . . , z n−1 ]dz .

therefore we can bound the norm of the gradient by
! Z

Z
∂g 2
ŝ(z) s(z) 2
1 2
n−1 2
≤
−
dz ·
4k[1, . . . , z
]k · |
| dz ≤ O(n/τ12 ) · g(Â, Ĉ) .
p(z)
p(z)
∂ Ĉ
T p̂(z)
T
Similarly, we could show that

3.2

∂g
∂â

2

≤ O(n2 /τ12 ) · g(Â, Ĉ).

Justifying idealized risk

We need to justify the approximation we made in Equation (3.1).

10

Lemma 3.5. Assume that ξt and xt are drawn i.i.d. from an arbitrary distribution with mean 0
b can be written as,
and variance 1. Then the population risk f (Θ)
2

b = (D̂ − D) +
f (Θ)

T
−1 
X
k=1

k
1−
T



Ĉ Âk−1 B − CAk−1 B

2

+ σ2 .

(3.8)

Proof of Lemma 3.5. Under the distributional assumptions on ξt and xt , we can calculate the objective functions above analytically. We write out yt , ŷt in terms of the inputs,
yt = Dxt +

t−1
X

t−k−1

CA

t−1

Bxk + CA

h0 + ξt ,

ŷt = D̂xt +

k=1

t−1
X

Ĉ Ât−k−1 B̂xk + CAt−1 h0 .

k=1

Therefore, using the fact that xt ’s are independent and with mean 0 and covariance I, the expectation of the error can be calculated (formally by Claim B.2),


Pt−1
2
2
E kŷt − yt k = kD̂ − DkF + k=1 Ĉ Ât−k−1 B̂ − CAt−k−1 B

2
F

+ E[kξt k2 ] .

2
F

+ σ2 .

(3.9)

Using E[kξt k2 ] = σ 2 , it follows that
b = kD̂ − Dk2F + PT −1 1 −
f (Θ)
k=1

k
T



Ĉ Âk−1 B̂ − CAk−1 B

(3.10)

Recall that under the controllable canonical form (1.4), B = en is known and therefore B̂ = B is
no longer a variable. We use â for the training variable corresponding to a. Then the expected
objective function (3.10) simplifies to
b = (D̂ − D)2 + PT −1 1 −
f (Θ)
k=1

k
T



Ĉ Âk−1 B − CAk−1 B

2

.

The previous lemma does not yet control higher order contributions present in the idealized
risk. This requires additional structure that we introduce in the next section.

4

Effective relaxations of spectral radius

The previous section showed quasi-convexity of the idealized risk. However, several steps are missing
towards showing finite sample guarantees for stochastic gradient descent. In particular, we will need
to control the variance of the stochastic gradient at any system that we encounter in the training.
For this purpose we formally introduce our main assumption now and show that it serves as an
effective relaxation of spectral radius. This results below will be used for proving convergence of
stochastic gradient descent in Section 5.
Consider the following convex region C in the complex plane,
C = {z : <z ≥ (1 + τ0 )|=z|} ∩ {z : τ1 < <z < τ2 } .

(4.1)

where τ0 , τ1 , τ2 > 0 are constants that are considered as fixed constant throughout the paper. Our
bounds will have polynomial dependency on these parameters.
Definition 4.1. We say a polynomial p(z) is α-acquiescent if {p(z)/z n : |z| = α} ⊆ C. A linear
system with transfer function G(z) = s(z)/p(z) is α-acquiescent if the denominator p(z) is.

11

The set of coefficients a ∈ Rn defining acquiescent systems form a convex set. Formally, for a
positive α > 0, define the convex set Bα ⊆ Rn as

Bα = a ∈ Rn : {pa (z)/z n : |z| = α} ⊆ C .
(4.2)

We note that definition (4.2) is equivalent to the definition Bα = a ∈ Rn : {z n p(1/z) : |z| =
1/α} ⊆ C , which is the version that we used in introduction for simplicity. Indeed, we can verify the
convexity of Bα by definition and the convexity of C: a, b ∈ Bα implies that pa (z)/z n , pb (z)/z n ∈ C
and therefore, p(a+b)/2 (z)/z n = 12 (pa (z)/z n + pb (z)/z n ) ∈ C. We also note that the parameter α
in the definition of acquiescence corresponds to the spectral radius of the companion matrix. In
particular, an acquiescent system is stable for α < 1.
Lemma 4.2. Suppose a ∈ Bα , then the roots of polynomial pa (z) have magnitudes bounded by α.
Therefore the controllable canonical form A = CC(a) defined by a has spectral radius ρ(A) < α.
Proof. Define holomorphic function f (z) = z n and g(z) = pa (z) = z n + a1 z n−1 + · · · + an . We apply
the symmetric form of Rouche’s theorem [Est62] on the circle K = {z : |z| = α}. For any point z
on K, we have that |f (z)| = αn , and that |f (z) − g(z)| = αn · |1 − pa (z)/z n |. Since a ∈ Bα , we have
that pa (z)/z n ∈ C for any z with |z| = α. Observe that for any c ∈ C we have that |1 − c| < 1 + |c|,
therefore we have that
|f (z) − g(z)| = αn |1 − pa (z)/z n | < αn (1 + |pa (z)|/|z n |) = |f (z)| + |pa (z)| = |f (z)| + |g(z)| .
Hence, using Rouche’s Theorem, we conclude that f and g have same number of roots inside
circle K. Note that function f = z n has exactly n roots in K and therefore g have all its n roots
inside circle K.
The following lemma establishes the fact that Bα is a monotone family of sets in α. The
proof follows from the maximum modulo principle of the harmonic functions <(z n p(1/z)) and
=(z n p(1/z)). We defer the short proof to Section C.1. We remark that there are larger convex sets
than Bα that ensure bounded spectral radius. However, in order to also guarantee monotonicity
and the no blow-up property below, we have to restrict our attention to Bα .
Lemma 4.3 (Monotonicity of Bα ). For any 0 < α < β, we have that Bα ⊂ Bβ .
Our next lemma entails that acquiescent systems have well behaved impulse responses.
Lemma 4.4 (No blow-up property). Suppose a ∈ Bα for some α ≤ 1. Then the companion matrix
A = CC(a) satisfies
∞
X
kα−k Ak Bk2 ≤ 2πnα−2n /τ12 .
(4.3)
k=0

Moreover, it holds that for any k ≥ 0,
kAk Bk2 ≤ min{2πn/τ12 , 2πnα2k−2n /τ12 } .
P
iλk α−k Ak B be the Fourier transform of the series α−k Ak B.
Proof of Lemma 4.4. Let fλ = ∞
k=0 e
Then using Parseval’s Theorem, we have
Z 2π
Z 2π
∞
X
−k k
2
2
kα A Bk =
|fλ | dλ =
|(I − α−1 eiλ A)−1 B|2 dλ
0

k=0

Z
=
0

0
2π

Pn
|pa

2j
j=1 α
dλ
(αe−iλ )|2

12

Z
≤
0

2π

n
dλ.
|pa (αe−iλ )|2

(4.4)

where at the last step we used the fact that (I − wA)−1 B =

1
[w−1 , w−2 . . . , z −n ]> (see
pa (w−1 )
|qa (α−1 eiλ )| ≥ τ1 , and therefore

Lemma B.1), and that α ≤ 1. Since a ∈ Bα , we have that
pa (αe−iλ ) = αn e−inλ q(eiλ /α) has magnitude at least τ1 αn . Plugging in this into equation (4.4), we
conclude that
Z 2π
∞
X
n
−k k
2
kα A Bk ≤
dλ ≤ 2πnα−2n /τ12 .
|pa (αe−iλ )|2
0
k=0

Finally we establish the bound for kAk Bk2 . By Lemma 4.3, we have Bα ⊂ B1 for α ≤ 1. Therefore
we can pick α = 1 in equation (4.3) and it still holds. That is, we have that
∞
X

kAk Bk2 ≤ 2πn/τ12 .

k=0

This also implies that kAk Bk2 ≤ 2πn/τ12 .

4.1

Efficiently computing the projection

In our algorithm, we require a projection onto Bα . However, the only requirement of the projection
step is that it projects onto a set contained inside Bα that also contains the true linear system.
So a variety of subroutines can be used to compute this projection or an approximation. First,
the explicit projection onto Bα is representable by a semidefinite program. This is because each
of the three constrains can be checked by testing if a trigonometric polynomial is non-negative.
A simple inner approximation can be constructed by requiring the constraints to hold on an a
finite grid of size O(n). One can check that this provides a tight, polyhedral approximation to the
set Bα , following an argument similar to Appendix C of Bhaskar et al [BTR13]. See Section F
for more detailed discussion on why projection on a polytope suffices. Furthermore, sometimes
we can replace the constraint by an `1 or `2 -constraint if we know that the system satisfies the
corresponding assumption. Removing the projection step entirely is an interesting open problem.

5

Learning acquiescent systems

Next we show that we can learn acquiescent systems.
Theorem 5.1. Suppose the true system Θ is α-acquiescent and satisfies kCk ≤ 1. Then with N
samples of length T ≥ Ω(n + 1/(1 − α)), stochastic gradient descent (Algorithm 1) with projection
b = (Â, B̂, Ĉ, D̂) with population risk
set Bα returns parameters Θ
!
r
2
5 + σ 2 n3
n
n
b ≤ f (Θ) + O
f (Θ)
+
,
(5.1)
N
TN
where O(·)-notation hides polynomial dependencies on 1/(1 − α), 1/τ0 , 1/τ1 , τ2 , and R = kak.
Recall that T is the length of the sequence and N is the number of samples. The first term
in the bound (5.1) comes from the smoothness of the population risk and the second comes from
the variance of the gradient estimator of population risk (which will be described in detail below).
3

Note that ỹt is different from ŷt defined in equation (1.2) which is used to define the population risk: here ŷt is
obtained from the (wrong) initial state h0 = 0 while ŷt is obtained from the correct initial state.
3
See Algorithm Box 3 for a detailed back-propagation algorithm that computes the gradient.

13

Algorithm 1 Projected stochastic gradient descent with partial loss
For i = 0 to N :
1. Take a fresh sample ((x1 , . . . , xT ), (y1 , . . . , yT )). Let ỹt be the simulated outputs2 of system
b on inputs x and initial states h0 = 0.
Θ
b =
2. Let TP
Run stochastic gradient descent3 on loss function `((x, y), Θ)
1 = T /4.
1
∂`
∂`
∂`
2
, and , GD =
, we update
t>T1 kỹt − yt k . Concretely, let GA = ∂â , GC =
T −T1
∂ Ĉ

∂ D̂

[â, Ĉ, D̂] → [â, Ĉ, D̂] − η[GA , GC , GD ] .
b = (â, Ĉ, D̂) to the set Bα ⊗ Rn ⊗ R.
3. Project Θ

An important (but not surprising) feature here is the variance
scale in 1/T and therefore for long
√
sequence actually we got 1/N convergence instead of 1/ N (for relatively small N ).
We can further balance the variance of the estimator with the number of samples by breaking
each long sequence of length T into Θ(T /n) short sequences of length Θ(n), and then run backpropagation (1) on these T N/n shorter sequences. This leads us to the following bound which gives
the right dependency in T and N as we expected: T N should be counted as the true number of
samples for the sequence-to-sequence model.
b with
Corollary 5.2. Under the assumption of Theorem 5.1, Algorithm 2 returns parameters Θ
population risk
!
r
5 + σ 2 n3
n
b ≤ f (Θ) + O
f (Θ)
,
TN
where O(·)-notation hides polynomial dependencies on 1/(1 − α), 1/τ0 , 1/τ1 , τ2 , and R = kak.
Algorithm 2 Projected stochastic gradient descent for long sequences
Input: N samples sequences of length T
b
Output: Learned system Θ
1. Divide each sample of length T into T /(βn) samples of length βn where β is a large enough
b
constant. Then run algorithm 1 with the new samples and obtain Θ.
We remark the the gradient computation procedure takes time linear in T n since one can use
chain-rule (also called back-propagation) to compute the gradient efficiently . For completeness,
Algorithm 3 gives a detailed implementation. Finally and importantly, we remark that although we
defined the population risk as the expected error with respected to sequence of length T , actually our
error bound generalizes to any longer (or shorter) sequences of length T 0  max{n, 1/(1 − α)}. By
b (Lemma 3.5) and the fact that kCAk Bk decays exponentially for k  n
the explicit formula for f (Θ)
(Lemma 4.4), we can bound the population risk on sequences of different lengths. Concretely, let
b denote the population risk on sequence of length T 0 , we have for all T 0  max{n, 1/(1 − α)},
fT 0 (Θ)
!
r
5 + σ 2 n3
n
b ≤ 1.1f (Θ)
b + exp(−(1 − α) min{T, T 0 }) ≤ O
fT 0 (Θ)
.
TN
14

We note that generalization to longer sequence does deserve attention. Indeed in practice, it’s
usually difficult to train non-linear recurrent networks that generalize to longer sequences than the
training data.
Our proof of Theorem 5.1 simply consists of three parts: a) showing the idealized risk is quasiconvex in the convex set Bα (Lemma 5.3); b) designing an (almost) unbiased estimator of the gradient of the idealized risk (Lemma 5.4); c) variance bounds of the gradient estimator (Lemma 5.5).
First of all, using the theory developed in Section 3 (Lemma 3.3 and Lemma 3.4), it is straightforward to verify that in the convex set Bα ⊗ Rn , the idealized risk is both weakly-quasi-convex and
weakly-smooth.
Lemma 5.3. Under the condition of Theorem 5.1, the idealized risk (3.2) is τ -weakly-quasi-convex
in the convex set Bα ⊗ Rn and Γ-weakly smooth, where τ = Ω(τ0 τ1 /τ2 ) and Γ = O(n2 /τ14 ).
Proof of Lemma 5.3. It suffices to show that for all â, a ∈ Bα , it satisfies â ∈ Nτ (a) for τ =
Ω(τ0 τ1 /τ2 ). Indeed, by the monotonicity of the family of sets Bα (Lemma 4.3), we have that
â, a ∈ B1 , which by definition means for every z on unit circle, pa (z)/z n , pâ (z)/z n ∈ C. By definition
of C, for any point w, ŵ ∈ C, the angle φ between w and ŵ is at most π − Ω(τ0 ) and ratio
of the magnitude is at least τ1 /τ2 , which implies that <(w/ŵ) = |w|/|ŵ| · cos(φ) ≥ Ω(τ0 τ1 /τ2 ).
Therefore <(pa (z)/pâ (z)) ≥ Ω(τ0 τ1 /τ2 ), and we conclude that â ∈ Nτ (a). The smoothness bound
was established in Lemma 3.4.
Towards designing an unbiased estimator of the gradient, we note that there is a small caveat
here that prevents us to just use the gradient of the empirical risk, as commonly done for other
(static) problems. Recall that the population risk is defined as the expected risk with known initial
state h0 , while in the training we don’t have access to the initial states and therefore using the
naive approach we couldn’t even estimate population risk from samples without knowing the initial
states.
We argue that being able to handle the missing initial states is indeed desired: in most of the
interesting applications h0 is unknown (or even to be learned). Moreover, the ability of handling
unknown h0 allows us to break a very long sequence into shorter sequences, which helps us to
obtain Corollary 5.2. Here the difficulty is essentially that we have a supervised learning problem
with missing data h0 . We get around it by simply ignoring first T1 = Ω(T ) outputs of the system
and setting the corresponding errors to 0. Since the influence of h0 to any outputs later than time
k ≥ T1  max{n, 1/(1 − α)} is inverse exponentially small, we could safely assume h0 = 0 when
the error earlier than time T1 is not taken into account.
This small trick also makes our algorithm suitable to the cases when these early outputs are
actually not observed. This is indeed an interesting setting, since in many sequence-to-sequence
model [SVL14], there is no output in the first half fraction of iterations (of course these models
have non-linear operation that we cannot handle).
The proof of the correctness of the estimator is almost trivial and deferred to Section C.
Lemma 5.4. Under the assumption of Theorem 5.1, suppose â, a ∈ Bα . Then in Algorithm 1,
at each iteration, GA , GC are unbiased estimators of the gradient of the idealized risk (3.2) in the
sense that:


∂g ∂g
,
± exp(−Ω((1 − α)T )) .
E [GA , GC ] =
∂â ∂ Ĉ
(5.2)
Finally, we control the variance of the gradient estimator.
15

Lemma 5.5. The (almost) unbiased estimator (GA , GC ) of the gradient of g(Â, Ĉ) has variance
bounded by

O n3 Λ2 /τ16 + σ 2 n2 Λ/τ14
Var [GA ] + Var [GC ] ≤
.
T
where Λ = O(max{n, 1/(1 − α) log 1/(1 − α)}).
Note that Lemma 5.5 does not directly follow from the Γ-weakly-smoothness of the population
b is also Γ-smooth for every sample.
risk, since it’s not clear whether the loss function `((x, y), Θ)
Moreover, even if it could work out, from smoothness the variance bound can be only as small as Γ2 ,
while the true variance scales linearly in 1/T . Here the discrepancy comes from that smoothness
implies an upper bound of the expected squared norm of the gradient, which is equal to the variance
plus the expected squared mean. Though typically for many other problems variance is on the same
order as the squared mean, here for our sequence-to-sequence model, actually the variance decreases
in length of the data, and therefore the bound of variance from smoothness is very pessimistic.
We bound directly the variance instead. It’s tedious but very simple in spirit. We mainly
need Lemma 4.4 to control various difference sums that shows up from calculating the expectation.
The only tricky part is to obtain the 1/T dependency which corresponds to the cancellation of
the contribution from the cross terms. In the proof we will basically write out the variance as a
(complicated) function of Â, Ĉ which consists of sums of terms involving (Ĉ Âk B − CAk B) and
Âk B. We control these sums using Lemma 4.4. The proof is deferred to Section C.
Finally we are ready to prove Theorem 5.1. We essentially just combine Lemma 5.3, Lemma 5.4
and Lemma 5.5 with the generic convergence Proposition 2.3. This will give us low error in idealized
risk and then we relate the idealized risk to the population risk.
Proof of Theorem 5.1. We consider g 0 (Â, Ĉ, D̂) = (D̂ − D)2 + g(Â, Ĉ), an extended version of
the idealized risk which takes the contribution of D̂ into account. By Lemma 5.4 we have that
Algorithm 1 computes GA , GC which are almost unbiased estimators of the gradients of g 0 up
to negligible error exp(−Ω((1 − α)T )), and by Lemma C.2 we have GD is an unbiased estimator
of g 0 with respect to D̂. Moreover by Lemma 5.5, these unbiased estimator has total variance
O(n5 +σ 2 n3 )
V =
where O(·) hides dependency on τ1 and (1 − α). Applying Proposition 2.3 (which
T
only requires an unbiased estimator of the
of g 0
), we obtain that after T iterations, we
 gradient
q
2
5
2 n3
b ≤
converge to a point with g 0 (â, Ĉ, D̂) ≤ O nN + n +σ
. Then, by Lemma 3.5 we have f (Θ)
TN


q
n2
n5 +σ 2 n3
0
2
0
+ f (Θ) which completes the proof.
g (â, Ĉ, D̂) + σ = g (â, Ĉ, D̂) + f (Θ) ≤ O N +
TN

6

The power of improper learning

We observe an interesting and important fact about the theory in Section 5: it solely requires
a condition on the characteristic function p(z). This suggests that the geometry of the training
objective function depends mostly on the denominator of the transfer function, even though the
system is uniquely determined by the transfer function G(z) = s(z)/p(z). This might seem to be an
undesirable discrepancy between the behavior of the system and our analysis of the optimization
problem.
However, we can actually exploit this discrepancy to design improper learning algorithms that
succeed under much weaker assumptions. We rely on the following simple observation about the
s(z)
. For an arbitrary polynomial u(z) of leading coefficient 1, we
invariance of a system G(z) = p(z)

16

can write G(z) as
G(z) =

s̃(z)
s(z)u(z)
=
,
p(z)u(z)
p̃(z)

where s̃ = su and p̃ = pu. Therefore the system s̃(z)/p̃(z) has identical behavior as G. Although
this is a redundant representation of G(z), it should counted as an acceptable solution. After all,
learning the minimum representation4 of linear system is impossible in general. In fact, we will
encounter an example in Section 6.1.
While not changing the behavior of the system, the extension from p(z) to p̃(z), does affect the
geometry of the optimization problem. In particular, if p̃(z) is now an α-acquiescent characteristic
polynomial as defined in Definition 4.1, then we could find it simply using stochastic gradient
descent as shown in Section 5. Observe that we don’t require knowledge of u(z) but only its
existence. Denoting by d the degree of u, the algorithm itself is simply stochastic gradient descent
with n + d model parameters instead of n.
Our discussion motivates the following definition.
Definition 6.1. A polynomial p(z) of degree n is α-acquiescent by extension of degree d if there
exists a polynomial u(z) of degree d and leading coefficient 1 such that p(z)u(z) is α-acquiescent.
For a transfer function G(z), we define it’s H2 norm as
kGk2H2

1
=
2π

Z

2π

|G(eiθ )|2 dθ .

0

We assume (with loss of generality) that the true transfer function G(z) has bounded H2 norm,
that is, kGkH2 ≤ 1. This can be achieve by a rescaling5 of the matrix C.
Theorem 6.2. Suppose the true system has transfer function G(z) = s(z)/p(z) with a characteristic function p(z) that is α-acquiescent by extension of degree d, and kGkH2 ≤ 1, then projected
stochastic gradient descent with m = n + d states (that is, Algorithm 2 with m states) returns a
b with population risk
system Θ
!
r
5 + σ 2 m3
m
b ≤O
f (Θ)
.
TK
where the O(·) notation hides polynomial dependencies on τ0 , τ1 , τ2 , 1/(1 − α).
The theorem follows directly from Corollary 5.2 (with some additional care about the scaling.
Proof of Theorem 6.2. Let p̃(z) = p(z)u(z) be the acquiescent extension of p(z). Since τ2 ≥
|u(z)p(z)| = |p̃(z)| ≥ τ0 on the unit circle, we have that |s̃(z)| = |s(z)||u(z)| = s(z) · Oτ (1/p(z)).
Therefore we have that s̃(z) satisfies that ks̃kH2 = Oτ (ks(z)/p(z)kH2 ) = Oτ (kG(z)kH2 ) ≤ Oτ (1).
That means that the vector C that determines the coefficients of s̃ satisfies that kCk ≤ Oτ (1),
since for a polynomial h(z) = b0 + · · · + bn−1 z n−1 , we have khkH2 = kbk. Therefore we can apply
Corollary 5.2 to complete the proof.
In the rest of this section, we discuss in subsection 6.1 the instability of the minimum representation in subsection, and in subsection 6.2 we show several examples where the characteristic
function p(z) is not α-acquiescent but is α-acquiescent by extension with small degree d.
4

The minimum representation of a transfer function G(z) is defined as the representation G(z) = s(z)/p(z) with
p(z) having minimum degree.
5
In fact, this is a natural scaling that makes comparing error easier. Recall that the population risk is essentially
kĜ − GkH2 , therefore rescaling C so that kGkH2 = 1 implies that when error  1 we achieve non-trivial performance.

17

As a final remark, the examples illustrated in the following sub-sections may be far from optimally analyzed. It is beyond the scope of this paper to understand the optimal condition under
which p(z) is acquiescent by extension.

6.1

Instability of the minimum representation

We begin by constructing a contrived example where the minimum representation of G(z) is not
stable at all and as a consequence one can’t hope to recover the minimum representation of G(z).
s(z)
s0 (z)
z n −0.8−n
1
0
Consider G(z) = p(z)
:= (z−0.1)(z
n −0.9−n ) and G (z) = p0 (z) := z−0.1 . Clearly these are the
minimum representations of the G(z) and G0 (z), which also both satisfy acquiescence. On the
one hand, the characteristic polynomial p(z) and p0 (z) are very different. On the other hand, the
transfer functions G(z) and G0 (z) have almost the same values on unit circle up to exponentially
small error,
0.8−n − 0.9−n
|G(z) − G0 (z)| ≤
≤ exp(−Ω(n)) .
(z − 0.1)(z − 0.9−n )
Moreover, the transfer functions G(z) and Ĝ(z) are on the order of Θ(1) on unit circle. These suggest
that from an (inverse polynomially accurate) approximation of the transfer function G(z), we cannot
hope to recover the minimum representation in any sense, even if the minimum representation
satisfies acquiescence.

6.2

Power of improper learning in various cases

We illustrate the use of improper learning through various examples below.
6.2.1

Example: artificial construction

We consider a simple contrived example where improper learning can help us learn the transfer
function dramatically. We will show an example of characteristic function which is not 1-acquiescent
but (α + 1)/2-(α + 1)/2-acquiescent by extension of degree 3.
Let n be a large enough
integer and α be a constant. Let J = {1, n − 1, n} and ω = e2πi/n , and
Q
3
j
then define p(z) = z
j∈[n],j ∈J
/ (z − αω ). Therefore we have that
p(z)/z n = z 3

Y

(1 − αω j /z) =

j∈[n],j∈J

1 − αn /z n
(1 − ω/z)(1 − ω −1 /z)(1 − 1/z)

(6.1)

Taking z = e−iπ/2 we have that p(z)/z n has argument (phase) roughly −3π/4, and therefore
it’s not in C, which implies that p(z) is not 1-acquiescent. On the other hand, picking u(z) =
(z−ω)(z−1)(z−ω −1 ) as the helper function, from equation (6.1) we have p(z)u(z)/z n+3 = 1−αn /z n
takes values inverse exponentially close to 1 on the circle with radius (α + 1)/2. Therefore p(z)u(z)
is (α + 1)/2-acquiescent.
6.2.2

Example: characteristic function with separated roots

A characteristic polynomial with well separated roots will be acquiescent by extension. Our bound
will depend on the following quantity of p that characterizes the separateness of the roots.

18

Definition 6.3. For a polynomial h(z) of degree n with roots λ1 , . . . , λn inside unit circle, define
the quantity Γ(·) of the polynomial h as:

Γ(h) :=

X
Q
j∈[n]

λnj
.
i6=j (λi − λj )

Lemma 6.4. Suppose p(z) is a polynomial of degree n with distinct roots inside circle with radius
√
α. Let Γ = Γ(p), then p(z) is α-acquiescent by extension of degree d = O(max{(1 − α)−1 log( nΓ ·
kpkH2 ), 0}).
Our main idea to extend p(z) by multiplying some polynomial u that approximates p−1 (in a
relatively weak sense) and therefore pu will always take values in the set C. We believe the following
lemma should be known though for completeness we provide the proof in Section D.
Lemma 6.5 (Approximation of inverse of a polynomial). Suppose p(z) is a polynomial of degree
n and leading coefficient 1 with distinct roots inside circle with radius α, and Γ = Γ(p). Then for
1
Γ
d = O(max{( 1−α
log (1−α)ζ
, 0}), there exists a polynomial h(z) of degree d and leading coefficient
1 such that for all z on unit circle,
z n+d
− h(z) ≤ ζ .
p(z)
Proof of Lemma 6.4. Let γ = 1 − α. Using Lemma 6.5 with ζ = 0.5kpk−1
H∞ , we have that there
1
exists polynomial u of degree d = O(max{ 1−α log(ΓkpkH∞ ), 0}) such that
z n+d
− u(z) ≤ ζ .
p(z)
Then we have that
p(z)u(z)/z n+d − 1 ≤ ζ|p(z)| < 0.5 .
Therefore p(z)u(z)/z n+d ∈ Cτ0 ,τ1 ,τ2 for constant τ0 , τ1 , τ2 . Finally noting that for degree n polyno√
mial we have khkH∞ ≤ n · khkH2 , which completes the proof.
6.2.3

Example: Characteristic polynomial with random roots

We consider the following generative model for characteristic polynomial of degree 2n. We generate
n complex numbers λ1 , . . . , λn uniformly randomly on circle with radius α < 1, and take λi , λ̄i for
i = 1, . . . , n as the roots of p(z). That is, p(z) = (z − λ1 )(z − λ̄1 ) . . . (z − λn )(z − λ̄n ). We show that
with good probability (over the randomness of λi ’s), polynomial p(z) will satisfy the condition in
subsection 6.2.2 so that it can be learned efficiently by our improper learning algorithm.
Theorem 6.6. Suppose p(z) with random roots inside circle of radius α is generated from the
process described above. Then with high probability over the choice of p, we have that Γ(p) ≤
e √n)) and kpkH ≤ exp(Õ(√n)). As a corollary, p(z) is α-acquiescent by extension of degree
exp(O(
2
e
O((1
− α)−1 n).
Towards proving Theorem 6.6, we need the following lemma about the expected distance of two
random points with radius ρ and r in log-space.
19

Lemma 6.7. Let x ∈ C be a fixed point with |x| = ρ, and λ uniformly drawn on the circle with
radius r. Then E [ln |x − λ|] = ln max{ρ, r} .
Proof. When r 6= ρ, let N be an integer and ω = e2iπ/N . Then we have that
N
1 X
ln |x − rω k |
N →∞ N

E[ln |x − λ| | r] = lim

(6.2)

k=1

P
k
The right hand of equation (6.2) can be computed easily by observing that N1 N
k=1 ln |x − rω | =
Q
P
N
N
1
1
1
k
N
N
k=1 (x − rω ) = N ln |x − r |. Therefore, when ρ > r, we have limN →∞ N
k=1 ln |x −
N ln
rω k | = limN →∞ ρ + N1 ln |(x/ρ)N − (r/ρ)N | = ln ρ. On the other hand, when ρ < r, we have
P
k
that limN →∞ N1 N
k=1 ln |x − rω | = ln r. Therefore we have that E[ln |x − λ| | r] = ln(max ρ, r).
For ρ = r, similarly proof (with more careful concern of regularity condition) we can show that
E[ln |x − λ| | r] = ln r.
Now we are ready to prove Theorem 6.6.
Proof of Theorem 6.6. Fixing index i, and the choice of λi , we consider the random variable Yi =
2n
P
ln( Q |λ −λ|λi| |Q |λ −λ¯ | )n ln |λi |− j6=i ln |λi −λj |. By Lemma 6.7, we have that E[Yi ] = n ln |λi |−
P j6=i i j j6=i i j
j6=i E[ln |λi − λj |] = ln(1 − δ). Let Zj = ln |λi − λj |. Then we have that Zj are random variable
with mean 0 and ψ1 -Orlicz norm bounded by 1 since E[eln |λi −λj | − 1] ≤ 1. Therefore by Bernstein
inequality for sub-exponential tail random variable (for example, [LT13, Theorem 6.21]), we have
P
e √
e
that with high probability (1 − n−10 ), it holds that
j6=i Zj ≤ O( n) where O hides logarithmic
e √n).
factors. Therefore, with high probability, we have |Yi | ≤ O(
Finally we take union bound overPall i ∈ [n], and obtain that with high probability, for ∀i ∈
e √n)). With similar technique, we can
e √n), which implies that n exp(Yi ) ≤ exp(O(
[n], |Yi | ≤ O(
i=1
√
prove that kpkH2 ≤ exp(Õ( n).
6.2.4

Example: Passive systems

We will show that with improper learning we can learn almost all passive systems, an important
class of stable linear dynamical system as we discussed earlier. We start off with the definition of
a strict-input passive system.
Definition 6.8 (Passive System, c.f [KA10]). A SISO linear system is strict-input passive if and
only if for some τ0 > 0 and any z on unit circle, <(G(z)) ≥ τ0 .
In order to learn the passive system, we need to add assumptions in the definition of strict
passivity. To make it precise, we define the following subsets of complex plane: For positive
constant τ0 , τ1 , τ2 , define
Cτ+0 ,τ1 ,τ2 = {z ∈ C : |z| ≤ τ2 , <(z) ≥ τ1 , <(z) ≥ τ0 |=(z)| } .

(6.3)

We say a transfer function G(z) = s(z)/p(z) is (τ0 , τ1 , τ2 )-strict input passive if for any z on
unit circle we have G(z) ∈ Cτ+0 ,τ1 ,τ2 . Note that for small constant τ0 , τ1 and large constant τ2 , this
basically means the system is strict-input passive.
Now we are ready to state our main theorem in this subsection. We will prove that passive
systems could be learned improperly with a constant factor more states (dimensions), assuming
s(z) has all its roots strictly inside unit circles and Γ(s) ≤ exp(O(n)).
20

Theorem 6.9. Suppose G(z) = s(z)/p(z) is (τ0 , τ1 , τ2 )-strict-input passive. Moreover, suppose
the roots of s(z) have magnitudes inside circle with radius α and Γ = Γ(s) ≤ exp(O(n)) and
kpkH2 ≤ exp(O(n)). Then p(z) is α-acquiescent by extension of degree d = Oτ,α (n), and as a
consequence we can learn G(z) with n + d states in polynomial time.
Moreover, suppose in addition we assume that G(z) ∈ Cτ0 ,τ1 ,τ2 for every z on unit circle. Then
p(z) is α-acquiescent by extension of degree d = Oτ,α (n).
The proof of Theorem 6.9 is similar in spirit to that of Lemma 6.4, and is deferred to Section D.

6.3

Improper learning using linear regression

In this subsection, we show that under stronger assumption than α-acquiescent by extension, we
can improperly learn a linear dynamical system with linear regression, up to some fixed bias.
The basic idea is to fit a linear function that maps [xk−` , . . . , xk ] to yk . This is equivalent to a
dynamical system with ` hidden states and with the companion matrix A in (1.4) being chosen as
a` = 1 and a`−1 = · · · = a1 = 0. In this case, the hidden states exactly memorize all the previous `
inputs, and the output is a linear combination of the hidden states.
Equivalently, in the frequency space, this corresponds to fitting the transfer function G(z) =
`−1 +···+c
1
s(z)/p(z) with a rational function of the form c1 z z `−1
= c1 z −(`−1) + · · · + cn . The following is
a sufficient condition on the characteristic polynomial p(x) that guarantees the existence of such
fitting,
Definition 6.10. A polynomial p(z) of degree n is extremely-acquiescent by extension of degree d
with bias ε if there exists a polynomial u(z) of degree d and leading coefficient 1 such that for all
z on unit circle,
p(z)u(z)/z n+d − 1 ≤ ε

(6.4)

We remark that if p(z) is 1-acquiescent by extension of degree d, then there exists u(z) such
that p(z)u(z)/z n+d ∈ C. Therefore, equation (6.4) above is a much stronger requirement than
acquiescence by extension.6
When p(z) is extremely-acquiescent, we see that the transfer function G(z) = s(z)/p(z) can be
approximated by s(z)u(z)/z n+d up to bias ε. Let ` = n + d + 1 and s(z)u(z) = c1 z `−1 + · · · + c` .
Then we have that G(z) can be approximated by the following dynamical system of ` hidden states
with ε bias: we choose A = CC(a) with a` = 1 and a`−1 = · · · = a1 = 0, and C = [c1 , . . . , c` ]. As
we have argued previously, such a dynamical system simply memorizes all the previous ` inputs,
and therefore it is equivalent to linear regression from the feature [xk−` , . . . , xk ] to output yk .
Proposition 6.11 (Informal). If the true system G(z) = s(z)/p(z) satisfies that p(z) is extremelyacquiescent by extension of degree d. Then using linear regression we can learn mapping from
[xk−` , . . . , xk ] to yk with bias ε and polynomial sampling complexity.
We remark that with linear regression the bias ε will only go to zero as we increase the length `
of the feature, but not as we increase the number of samples. Moreover, linear regression requires
a stronger assumption than the improper learning results in previous subsections do. The latter
can be viewed as an interpolation between the proper case and the regime where linear regression
works.
6

We need (1−δ)-acquiescence by extension in previous subsections for small δ > 0, though this is merely additional
technicality needed for the sample complexity. We ignore this difference between 1−δ-acquiescence and 1-acquiescence
and for the purpose of this subsection

21

7

Learning multi-input multi-output (MIMO) systems

We consider multi-input multi-output systems with the transfer functions that have a common
denominator p(z),
1
G(z) =
· S(z)
(7.1)
p(z)
where S(z) is an `in × `out matrix with each entry being a polynomial with real coefficients of degree
at most n and p(z) = z n + a1 z n−1 + · · · + an . Note that here we use `in to denote the dimension of
the inputs of the system and `out the dimension of the outputs.
Although a special case of a general MIMO system, this class of systems still contains many
interesting cases, such as the transfer
functions studied in [FHB01, FHB04], where G(z) is assumed
P
Ri
to take the form G(z) = R0 + ni=1 z−λ
, for λ1 , . . . , λn ∈ C with conjugate symmetry and Ri ∈
i
`
×`
C out in satisfies that Ri = R̄j whenever λi = λ̄j .
In order to learn the system G(z), we parametrize p(z) by its coefficients a1 , . . . , an and S(z)
by the coefficients of its entries. Note that each entry of S(z) depends on n + 1 real coefficients and
therefore the collection of coefficients forms a third order tensor of dimension `out × `in × (n + 1). It
will be convenient to collect the leading coefficients of the entries of S(z) into a matrix of dimension
`out × `in , named D, and the rest of the coefficients into a matrix of dimension `out × `in n, denoted
by C. This will be particularly intuitive when a state-space representation is used to learn the
system with samples as discussed later. We parameterize the training transfer function Ĝ(z) by â,
Ĉ and D̂ using the same way.
Let’s define the risk function in the frequency domain as,
Z 2π
2
g(Â, Ĉ, D̂) =
G(eiθ ) − Ĝ(eiθ ) dθ .
(7.2)
F

0

The following lemma is an analog of Lemma 3.3 for the MIMO case. Itss proof actually follows
from a straightforward extension of the proof of Lemma 3.3 by observing that matrix S(z) (or Ŝ(z))
commute with scalar p(z) and p̂(z), and that Ŝ(z), p̂(z) are linear in â, Ĉ.
Lemma 7.1. The risk function g(â, Ĉ) defined in (7.2) is τ -weakly-quasi-convex in the domain




pa (z)
0
n
Nτ (a) = â ∈ R : <
≥ τ /2, ∀ z ∈ C, s.t. |z| = 1 ⊗ R`in ×`out ×n
pâ (z)
Finally, as alluded before, we use a particular state space representation for learning the system
in time domain with example sequences. It is known that any transfer function of the form (7.1)
can be realized uniquely by the state space system of the following special case of Brunovsky normal
form [Bru70],




0
I`in
0
···
0
0


0
0
I`in
···
0

 .. 

 . 


.
.
.
..
.
..
..
..
..
(7.3)
A = 
,
B
=
,


.


 0 


0
0
0
···
I`in
I`in
−an I`in −an−1 I`in −an−2 I`in · · · −a1 I`in
and,
C ∈ R`out ×n`in ,

D ∈ R`out ×`in .

The following Theorem is a straightforward extension of Corollary 5.2 and Theorem 6.2 to the
MIMO case.
22

Theorem 7.2. Suppose transfer function G(z) of a MIMO system takes form (7.1), and has norm
kGkH2 ≤ 1. If the common denominator p(z) is α-acquiescent by extension of degree d then
b with
projected stochastic gradient descent over the state space representation (7.3) will return Θ
risk
−1
b ≤ poly(n + d, σ, τ, (1 − α) ) .
f (Θ)
TN
We note that since A and B are simply the tensor product of I`in with CC(a) and en , the no
blow-up property (Lemma 4.4) for Ak B still remains true. Therefore to prove Theorem 7.2, we
essentially only need to run the proof of Lemma 5.5 with matrix notation and matrix norm. We
defer the proof to the full version.

Acknowledgments
We thank Amir Globerson, Alexandre Megretski, Pablo A. Parrilo, Yoram Singer and Ruixiang
Zhang for helpful discussions. We are indebted to Mark Tobenkin for pointers to relevant prior
work. We also thank Alexandre Megretski for helpful feedback, insights into passive systems and
suggestions on how to organize Section 3.

References
[BGMA08] Alexandre S. Bazanella, Michel Gevers, Ljubisa Miskovic, and Brian D.O. Anderson.
Iterative minimization of h2 control performance criteria. Automatica, 44:2549–2559,
2008.
[Bot98]

Léon Bottou. On-line learning in neural networks. chapter On-line Learning and
Stochastic Approximations, pages 9–42. Cambridge University Press, New York, NY,
USA, 1998.

[Bru70]

Pavol Brunovsky. A classification of linear controllable systems.
06(3):(173)–188, 1970.

[BTR13]

Badri Narayan Bhaskar, Gongguo Tang, and Benjamin Recht. Atomic norm denoising
with applications to line spectral estimation. IEEE Transactions on Signal Processing,
61(23):5987–5999, 2013.

[CW02]

M. C. Campi and Erik Weyer. Finite sample properties of system identification methods. IEEE Transactions on Automatic Control, 47(8):1329–1334, 2002.

[DWB06]

Hugh Durrant-Whyte and Tim Bailey. Simultaneous localization and mapping: part i.
Robotics & Automation Magazine, IEEE, 13(2):99–110, 2006.

[EB11]

Diego Eckhard and Alexandre Sanfelice Bazanella. On the global convergence of identification of output error models. In Proc. 18th IFAC World congress, 2011.

[Est62]

T. Estermann. Complex numbers and functions. Athlone Press, 1962.

[FHB01]

Maryam Fazel, Haitham Hindi, and Stephen P Boyd. A rank minimization heuristic
with application to minimum order system approximation. In Proc. American Control
Conference, volume 6, pages 4734–4739. IEEE, 2001.

23

Kybernetika,

[FHB04]

Maryam Fazel, Haitham Hindi, and S Boyd. Rank minimization and applications in
system theory. In Proc. American Control Conference, volume 4, pages 3273–3278.
IEEE, 2004.

[GHJY15] Rong Ge, Furong Huang, Chi Jin, and Yang Yuan. Escaping from saddle points - online
stochastic gradient for tensor decomposition. In Proc. 28th COLT, pages 797–842, 2015.
[Hes09]

Joao P Hespanha. Linear systems theory. Princeton university press, 2009.

[HLS15]

E. Hazan, K. Y. Levy, and S. Shalev-Shwartz. Beyond Convexity: Stochastic QuasiConvex Optimization. ArXiv e-prints, July 2015.

[HRS07]

Christiaan Heij, André Ran, and Freek van Schagen. Introduction to mathematical
systems theory : linear systems, identification and control. Birkhäuser, Basel, Boston,
Berlin, 2007.

[HS97]

Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Neural Computation, 9(8):1735–1780, 1997.

[KA10]

Nicholas Kottenstette and Panos J Antsaklis. Relationships between positive real,
passive dissipative, & positive systems. In American Control Conference (ACC), 2010,
pages 409–416. IEEE, 2010.

[Lju98]

Lennart Ljung. System Identification. Theory for the user. Prentice Hall, Upper Saddle
River, NJ, 2nd edition, 1998.

[LK13]

Sergey Levine and Vladlen Koltun. Guided policy search. In Proceedings of The 30th
International Conference on Machine Learning, pages 1–9, 2013.

[LSJR16]

J. D. Lee, M. Simchowitz, M. I. Jordan, and B. Recht. Gradient Descent Converges to
Minimizers. ArXiv e-prints, February 2016.

[LT13]

Michel Ledoux and Michel Talagrand. Probability in Banach Spaces: isoperimetry and
processes, volume 23. Springer Science & Business Media, 2013.

[Meg08]

Alexandre Megretski. Convex optimization in robust identification of nonlinear feedback. In Proceedings of the 47th Conference on Decision and Control, 2008.

[RRD05]

Ali Rahimi, Ben Recht, and Trevor Darrell. Learning appearance manifolds from video.
In Proc. IEEE CVPR, 2005.

[SBTR12] Parikshit Shah, Badri Narayan Bhaskar, Gongguo Tang, and Benjamin Recht. Linear system identification via atomic norm regularization. In Proceedings of the 51st
Conference on Decision and Control, 2012.
[Sch41]

A. C. Schaeffer. Inequalities of a. markoff and s. bernstein for polynomials and related
functions. Bull. Amer. Math. Soc., 47(8):565–579, 08 1941.

[SVL14]

Ilya Sutskever, Oriol Vinyals, and Quoc V. Le. Sequence to sequence learning with
neural networks. In Proc. 27th NIPS, pages 3104–3112, 2014.

[VK08]

M. Vidyasagar and Rajeeva L. Karandikar. A learning theory approach to system
identification and stochastic adaptive control. Journal of Process Control, 18(3):421–
430, 2008.
24

[WC99]

A

Erik Weyer and M. C. Campi. Finite sample properties of system identification methods. In Proceedings of the 38th Conference on Decision and Control, 1999.

Background on optimization

The proof below uses the standard analysis of gradient descent for non-smooth objectives and
demonstrates that the argument still works for weakly-quasi-convex functions.
Proof of Proposition 2.3. We start by using the weakly-quasi-convex condition and then the rest
follows a variant of the standard analysis of non-smooth projected sub-gradient descent7 . We
conditioned on θk , and have that
τ (f (θk ) − f (θ∗ )) ≤ ∇f (θk )> (θk − θ∗ ) = E[r(θk )> (θk − θ∗ ) | θk ]


1
∗
= E (θk − wk+1 )(θk − θ ) | θk
η




1
∗ 2
2
∗ 2
=
E kθk − wk+1 k | θk + kθk − θ k − E kwk+1 − θ k | θk
η

 1


= η E kr(θk )k2 +
kθk − θ∗ k2 − E kwk+1 − θ∗ k2 | θk
η

(A.1)

where the first inequality uses weakly-quasi-convex and the rest of lines are simply algebraic manipulations. Since θk+1 is the projection of wk+1 to B and θ∗ belongs to B, we have kwk+1 − θ∗ k ≥
kθk+1 − θ∗ k. Together with (A.1), and


2
∗
2
E kr(θk )k = k∇f (θk )k + Var[r(θk )] ≤ Γ(f (θk ) − f (θ )) + V,
we obtain that
τ (f (θk ) − f (θ∗ )) ≤ ηΓ(f (θk ) − f (θ∗ )) + ηV +



1
kθk − θ∗ k2 − E kθk+1 − θ∗ k2 | θk .
η

Taking expectation over all the randomness and summing over k we obtain that




K−1
X
1
1
1
1 2
∗
∗ 2
ηKV + kθ0 − θ k ≤
ηKV + R
.
E [f (θk ) − f (θ )] ≤
τ − ηΓ
η
τ − ηΓ
η
k=0

2 2

Γ
where we use the assumption that kθ0 − θ∗ k ≤ R. Suppose K ≥ 4R
, then we take η = √VRK .
V τ2
Therefore we have that τ − ηΓ ≥ τ /2 and therefore
√ √
K−1
X
4R V K
∗
.
(A.2)
E [f (θk ) − f (θ )] ≤
τ
k=0

4R2 Γ2
,
V τ2

τ
we pick η = 2Γ
and obtain that


K−1
X
2 τ KV
2ΓR2
8ΓR2
∗
+
≤
.
E [f (θk ) − f (θ )] ≤
τ
2Γ
τ
τ2

On the other hand, if K ≤

(A.3)

k=0

Therefore using equation (A.3) and (A.2) we obtain that when choosing η properly according
to K as above,
(
√ )
2 4R V
8ΓR
∗
, √
.
E [f (θk ) − f (θ )] ≤ max
τ 2K τ K
k∈[K]
7

Although we used weak smoothness to get a slightly better bound

25

B

Toolbox

Lemma B.1. Let B = en ∈ Rn×1 and λ ∈ [0, 2π], w ∈ C. Suppose A with ρ(A) · |w| < 1 has the
controllable canonical form A = CC(a). Then
 −1 
w
−2 

1
w 
−1
(I − wA) B =
 . 
pa (w−1 )  .. 
w−n
where pa (x) = xn + a1 xn−1 + · · · + an is the characteristic polynomial of A.
Proof. let v = (I − wA)−1 B then we have (I − wA)v = B.
form

1
−w
0
 0
1
−w

 ..
.
..
.
I − wA =  .
.
.

 0
0
0
an w an−1 w an−2 w

Note that B = en , and I − wA is of the

···
0

···
0


.
..
..
(B.1)

.


···
−w
· · · 1 + a1 w

Therefore we obtain that vk = wvk+1 for 1 ≤ k ≤ n − 1. That is, vk = v0 w−k for v0 = v1 w1 .
Using the fact that ((I − wA)v)n = 1, we obtain that v0 = pa (w−1 )−1 where pa (·) is the polynomial
−1
−n
nw
pa (x) = xn + a1 xn−1 + · · · + an . Then we have that u(I − wA)−1 B = u1 w pa+···+u
(w−1 )
Claim B.2. Suppose x1 , . . . , xn are independent variables with mean 0 and covariance matrices
and I, U1 , . . . , Ud are fixed matrices, then
 Pn
 Pn
E k k=1 Uk xk k2 = k=1 kUk k2F .
Proof. We have that
 Pn
 Pn
Pn
Pn
2
> >
>
E k k=1 Uk xk k2F = k,` tr(Uk xk x>
k=1 kUk kF
k tr(Uk xk xk Uk ) =
` U` ) =

C
C.1

Missing proofs in Sections 4 and 5
Monotonicity of acquiescence: Proof of Lemma 4.3

Lemma C.1 (Lemma 4.3 restated). For any 0 < α < β, we have that Bα ⊂ Bβ .
Proof. Let qa (z) = 1 + a1 z + · · · + an z n . Note that q(z −1 ) = pa (z)/z n . Therefore we note that
Bα = {a : qa (z) ∈ C, ∀|z| = 1/α}. Suppose a ∈ Bα , then <(qa (z)) ≥ τ1 for any z with |z| =
1/α. Since <(qa (z)) is the real part of the holomorphic function qa (z), its a harmonic function.
By maximum (minimum) principle of the harmonic functions, we have that for any |z| ≤ 1/α,
<(qa (z)) ≥ inf |z|=1/α <(qa (z)) ≥ τ1 . In particular, it holds that for |z| = 1/β < 1/α, <(qa (z)) ≥ τ1 .
Similarly we can prove that for z with |z| = 1/β, <(qa (z)) ≥ (1 + τ0 )=(qa (z)), and other conditions
for a being in Bβ .

26

C.2

Proof of Lemma 5.4

Lemma 5.4 follows directly from the following general Lemma which also handles the multi-input
multi-output case. It can be seen simply from calculation similar to the proof of Lemma 3.5. We
mainly need to control the tail of the series using the no-blow up property (Lemma 4.4) and argue
that the wrong value of the initial states h0 won’t cause any trouble to the partial loss function
b (defined in Algorithm 1). This is simply because after time T1 = T /4, the influence of
`((x, y), Θ)
the initial state is already washed out.
Lemma C.2. In algorithm 3 the values of GA , GC , GD are equal to the gradients of g(Â, Ĉ) + (D̂ −
D)2 with respect to Â, Ĉ and D̂ up to inverse exponentially small error.
b has expecProof of Lemma C.2. We first show that the partial empirical loss function `((x, y), Θ)
tation almost equal to the idealized risk (up to the term for D̂ and exponential small error),
b = g(Â, Ĉ) + (D̂ − D)2 ± exp(−Ω((1 − α)T )).
E[`((x, y), Θ)]
This can be seen simply from similar calculation to the proof of Lemma 3.5. Note that
yt = Dxt +

t−1
X

CAt−k−1 Bxk + CAt−1 h0 + ξt

and ỹt = D̂xt +

t−1
X

Ĉ Ât−k−1 B̂xk .

(C.1)

k=1

k=1

Therefore noting that when t ≥ T1 ≥ Ω(T ), we have that kCAt−1 h0 k ≤ exp(−Ω((1 − α)T ) and
therefore the effect of h0 is negligible. Then we have that


T
X
1
b =
kyt − yt k2  ± exp(−Ω((1 − α)T ))
E[`((x, y), Θ)]
E
T − T1
t>T1

= kD̂ − Dk2 +

= kD̂ − Dk2 +

1
T − T1
T1
X

X

X

T ≥t>T1 0≤j≤t−1

kĈ Âj B − CAj Bk2 +

= kD̂ − Dk2 +

X
T ≥j≥T1

j=0
∞
X

kĈ Âj B − CAj Bk2 ± exp(−Ω((1 − α)T ))
T −j
kĈ Âj B − CAj Bk2 ± exp(−Ω((1 − α)T ))
T − T1

kĈ Âj B − CAj Bk2 ± exp(−Ω((1 − α)T )) .

j=0

where the first line use the fact that kCAt−1 h0 k ≤ exp(−Ω((1−α)T ), the second uses equation (3.9)
and the last line uses the no-blowing up property of Ak B (Lemma 4.4).
b is also close to the gradient of
Similarly, we can prove that the gradient of E[`((x, y), Θ)]
2
g(Â, Ĉ) + (D̂ − D) up to inverse exponential error.

C.3

Proof of Lemma 5.5

Proof of Lemma 5.5. Both GA and GC can be written in the form of a quadratic form (with vector
coefficients) of x1 , . . . , xT and ξ1 , . . . , ξT . That is, we will write
X
X
X
X
0
GA =
xs xt ust +
xs ξt u0st and GC =
xs xt vst +
xs ξt vst
.
s,t

s,t

s,t

27

s,t

where ust and vst are vectors that will be calculated later. By Claim C.3, we have that
#
"
X
X
X
X
0
ku0st k2 .
(C.2)
kust k2 + O(σ 2 )
xs ξt ust ≤ O(1)
Var
xs xs ust +
s,t

s,t

s,t

s,t

P
P 0 2
Therefore in order to bound from above Var [GA ], it suffices to bound
kust k2 and
kust k ,
and similarly for GC .
We begin by writing out ust for fixed s, t ∈ [T ] and bounding its norm. We use the same set
of notations as int the proof of Lemma 5.4. Recall that we set rk = CAk B and r̂k = Ĉ Âk B, and
∆rk = r̂k − rk . Moreover, let zk = Âk B. We note that the sums of kzk k2 and rk2 can be controlled.
By the assumption of the Lemma, we have that
∞
X

kzk k2 ≤ 2πnτ1−2 ,

k=t
∞
X

∆rk2 ≤ 4πnτ1−2 ,

kzk k2 ≤ 2πnα2k−2n τ1−2 .

(C.3)

k∆rk k2 ≤ 4πnα2k−2n τ1−2 .

(C.4)

k=t

which will be used many times in the proof that follows.
We calculate the explicit form of GA using the explicit back-propagation Algorithm 3. We have
that in Algorithm 3,
P
P
(C.5)
ĥk = kj=1 Âk−j Bxj = kj=1 zk−j xj
and
T
X

∆hk =

(Â> )j−k Ĉ > ∆yj =

αj−k (Â> )j−k Ĉ > 1j>T1

ξj +

> ∆h

j
X

!
∆rj−` x`

(C.6)

`=1

j=k

j=k

P

T
X

0>
k hk−1

and equation (C.5) and equation (C.6) above, we have that
P

P
j−k B 1
k−t−1 B
ust = Tk=2
∆r
Ĉ
Â
j−s
k≥t+1 · Â
j≥max{k,s,T1 +1}

P
P
(C.7)
∆r
r̂
= Tk=2
j−s j−k 1k≥t+1 · zk−t−1 .
j≥max{k,s,T1 +1}

Then using GA =

k≥2 B

and that,
u0st =

T
X

0
zk−1−s · 1k≥s+1 · r̂t−k
· 1t>max{T1 ,k} =

k=2

X

0
zk−1−s · r̂t−k
· 1t>max{T1 }

(C.8)

s+1≤k≤t



1
Towards bounding kust k, we consider four different cases. Let Λ = Ω {max{n, (1 − α)−1 log( 1−α
)}
be a threshold.
When 0 ≤ s − t ≤ Λ, we rewrite ust by rearranging equation (C.7),
X
X
X
X
zk−t−1
∆rj−s r̂j−k +
zk−t−1
∆rj−s r̂j−k
ust =

Case 1:

T ≥k≥s

=

X
`≥0,`≥T1 +1−s

t<k<s

j≥max{k,T1 +1}

∆r`

X

r̂`+s−k zk−t−1 +

s≤k≤l+s,k≤T

j≥max{s,T1 +1}

X
`≥0,`≥T1 +1−s

28

∆r`

X
s>k>t

r̂`+s−k zk−t−1

where at the second line, we did the change of variables ` = j − s. Then by Cauchy-Schartz
inequality, we have,


2
X
X
X
kust k2 ≤ 2 
∆r`2  
r̂`+s−k zk−t−1 
`≥0,`≥T1 +1−s

`≥0,`≥T1 +1−s s≤k≤l+s,k≤T

{z

|



X

+ 2

}

T1

2

X

∆r`2  

`≥0,`≥T1 +1−s

X

r̂`+s−k zk−t−1


.

(C.9)

`≥0,`≥T1 +1−s s>k>t

{z

|

T2

}

We could bound the contribution from ∆rk2 ssing equation (C.4), and it remains to bound terms
T1 and T2 . Using the tail bounds for kzk k (equation (C.3)) and the fact that |r̂k | = |Ĉ Âk B| ≤
kÂk Bk = kzk k , we have that

2
2
X
X
X
X

T1 =
r̂`+s−k zk−t−1 ≤
|r̂`+s−k |kzk−t−1 k .
(C.10)
`≥0,`≥T1 +1−s s≤k≤l+s,k≤T

`≥0

s≤k≤`+s

We bound the inner sum of RHS of (C.10) using the fact that kzk k2 ≤ O(nα2k−2n /τ12 ) and
obtain that,
X
X
|r̂`+s−k |kzk−t−1 k ≤
O(nα(`+s−t−1)−2n /τ12 )
s≤k≤`+s

s≤k≤`+s

≤ O(`nα(`+s−t−1)−2n /τ12 ) .

(C.11)

Note that equation (C.11) is particular effective when ` > Λ. When ` ≤ Λ, we can refine the bound
using equation (C.3) and obtain that

1/2 
1/2
X
X
X
|r̂`+s−k |kzk−t−1 k ≤ 
|r̂`+s−k |2  
kzk−t−1 k2 
s≤k≤`+s

s≤k≤`+s

s≤k≤`+s

√
√
≤ O( n/τ1 ) · O( n/τ1 ) = O(n/τ12 ) .

(C.12)

Plugging equation (C.12) and (C.11) into equation (C.10), we have that

2
X
X
X
X

|r̂`+s−k |kzk−t−1 k ≤
O(n2 /τ14 ) +
O(`2 n2 α2(`+s−t−1)−4n /τ14 )
`≥0

s≤k≤`+s

Λ≥`≥0

≤ O(n

2

Λ/τ14 )

`>Λ
2

+ O(n /τ14 ) = O(n2 Λ/τ14 ) .

(C.13)

For the second term in equation (C.9), we bound similarly,
2

T2 ≤

X

X

r̂`+s−k zk−t−1

≤ O(n2 Λ/τ14 ) .

(C.14)

`≥0,`≥T1 +1−s s>k>t

Therefore using the bounds for T1 and T2 we obtain that,
kust k2 ≤ O(n3 Λ/τ16 )
29

(C.15)

Case 2:

When s − t > Λ, we tighten equation (C.13) by observing that,

2
X
X
X

O(`2 n2 α2` /τ14 )
T1 ≤
|r̂`+s−k |kzk−t−1 k ≤ α2(s−t−1)−4n
`≥0

`≥0

s≤k≤`+s

2

(C.16)

kust k2 ≤ O(n3 /((1 − α)3 τ16 )) · αs−t−1 .

(C.17)

≤α

s−t−1

· O(n /(τ14 (1 − α)3 )) .

where we used equation (C.11). Similarly we can prove that
T2 ≤ αs−t−1 · O(n2 /(τ14 (1 − α)3 )) .
Therefore, we have when s − t ≥ Λ,

Case 3: When −Λ ≤ s − t ≤ 0, we can rewrite us t and use the Cauchy-Schwartz inequality and
obtain that
X
X
X
X
r̂`+s−k zk−t−1 .
ust =
zk−t−1
∆rj−s r̂j−k =
∆r`
T ≥k≥t+1

`≥0,`≥T1 +1−s

j≥max{k,T1 +1}

t+1≤k≤l+s,k≤T

and,
kust k2 ≤ 

2




X

X

∆r`2  

X

r̂`+s−k zk−t−1

.

`≥0,`≥T1 +1−s t+1≤k≤l+s,k≤T

`≥0,`≥T1 +1−s

Using almost the same arguments as in equation (C.11) and (C.12), we that
X
|r̂`+s−k | · kzk−t−1 k ≤ O(`nα(`+s−t−1)−2n /τ12 )
t+1≤k≤`+s

X

and

√
√
|r̂`+s−k | · kzk−t−1 k ≤ O( n/τ1 ) · O( n/τ1 ) = O(n/τ12 ) .

t+1≤k≤`+s

Then using a same type of argument as equation (C.13), we can have that
2

X

X

0
0
r̂`+s−k
zk−t−1

≤ O(n2 Λ/τ14 ) + O(n2 /τ14 )

`≥0,`≥T1 +1−s t+1≤k≤l+s,k≤T

= O(n2 Λ/τ14 ) .
It follows that in this case kust k can be bounded with the same bound in (C.15).
Case 4: When s − t ≤ −Λ, we use a different simplification of ust from above. First of all, it
follows (C.7) that


T
X
X
0

kust k ≤
k∆rj−s r̂j−k
zk−t−1 k1k≥t+1 
k=2

j≥max{k,s,T1 +1}

(C.18)
≤

X
k≥t+1

0
kzk−t−1
k

X
j≥max{k,T1 +1}

30

0
|∆rj−s r̂j−k
|.

Since j − s ≥ k − s > 4n and it follows that
X
X
0
|∆rj−s r̂j−k
|≤

√
√
O( n/τ1 · αj−s−n ) · O( n/τ1 · αj−k−n )

j≥max{k,T1 +1}

j≥max{k,T1 +1}

≤ O(n/(τ12 (1 − α)) · αk−s−n )
Then we have that
X

kust k2 ≤

X

0
kzk−t−1
k

k≥t+1

j≥max{k,T1 +1}


≤


X

X


k≥t+1

O(n/τ12 )

· O(n

2

/(τ14 (1

2 


X

0
kzk−t−1
k2  

k≥t+1

≤

0
|∆rj−s r̂j−k
|

0
|∆rj−s r̂j−k
| 

j≥max{k,T1 +1}

3

− α) )αt−s ) = O(n3 /(τ16 δ 3 )αt−s )

Therefore, using the bound for kust k2 obtained in the four cases above, taking sum over s, t, we
obtain that
X
X
X
kust k2 ≤
O(n3 Λ/τ16 ) +
O(n3 /(τ16 (1 − α)3 )α|t−s|−1 )
1≤s,t≤T

s,t∈[T ]:|s−t|≤Λ
3

≤ O(T n Λ

2

/τ16 )

s,t:|s−t|≥Λ

+ O(n

3

/τ16 )

= O(T n3 Λ2 /τ16 ) .

(C.19)

We finished the bounds for kust k and now we turn to bound ku0st k2 . Using the formula for u0st
(equation C.8), we have that for t ≤ s + 1, u0st = 0. For s + Λ ≥ t ≥ s + 2, we have that by
Cauchy-Schwartz inequality,
1/2

1/2 
X
X
0
|2  ≤ O(n/τ12 ) ≤ O(n/τ12 ) .
ku0st k ≤ 
kzk−1−s k2  
|r̂t−k
s+1≤k≤t

s+1≤k≤t

On the other hand, for t > s + Λ, by the bound that |r̂k0 |2 ≤ kzk0 k2 ≤ O(nα2k−2n /τ12 ), we have,
ku0st k

≤

T
X

kzk−1−s k ·

0
|
|r̂t−k

nαt−s−1 /τ12

s+1≤k≤t−1

s+1≤k≤t−1

≤ O(n(t − s)α

≤

T
X

t−s−1

/τ12 ) .

Therefore taking sum over s, t, similarly to equation (C.19),
X
ku0st k2 ≤ O(T n2 Λ/τ14 ) .
s,t∈[T ]

Then using equation (C.2) and equation (C.19) and (C.20), we obtain that

Var[kGA k2 ] ≤ O T n3 Λ2 /τ16 + σ 2 T n2 Λ/τ14 .
Hence, it follows that

O n3 Λ2 /τ16 + σ 2 n2 Λ/τ14
1
Var[GA ] ≤
.
Var[GA ] ≤
(T − T1 )2
T
We can prove the bound for GC similarly.
31

(C.20)

Claim C.3. Let x1 , . . . , xT be independent random variables with mean 0 and variance 1 and 4-th
moment bounded by O(1), and uij be vectors for i, j ∈ [T ]. Moreover, let ξ1 , . . . , ξT be independent
random variables with mean 0 and variance σ 2 and u0ij be vectors for i, j ∈ [T ]. Then,
Var

hP

i,j xi xj uij +

P

i
P
P
0
2
2
0 2
x
ξ
u
i,j i j ij ≤ O(1)
i,j kuij k + O(σ )
i,j kuij k .

Proof. Note that the two sums in the target are independent with mean 0, therefore we only need
to bound the variance of both sums individually. The proof follows the linearity of expectation and
the independence of xi ’s:

 XX h
i
2
P
>
=
x
x
x
x
u
u
x
x
u
E
E
i j k ` ij k`
i,j i j ij
i,j

=

X

≤

X

k,`

X

>
4
E[uii uii xi ] +

i

X

+ O(1)

X

h
i
2 2 >
>
x
x
(u
u
+
u
u
)
E i j ij ij
ij ji

i,j

i6=j

u>
ii ujj

i,j

=k

>
2 2
E[uii ujj xi xj ] +
2

kuij + uji k

i,j
2
i uii k + O(1)

P

X

kuij k2

i,j
α
where at second line we usedP
the fact that for
Pany monomial x with an odd degree on one of the
α
xi ’s, E[x ] = 0. Note that E[ i,j xi xj uij ] = i uii . Therefore,

Var

hP

i
h P
i
P
P
2 −k [
x
x
u
k
k
x
x
u
=
E i,j xi xj uij ]k2 ≤ O(1) i,j kuij k2
E
i,j i j ij
i,j i j ij

Similarly, we can control Var

D

(C.21)

hP

i
P
0 2
0
2
x
ξ
u
i,j kuij k .
i,j i j ij by O(σ )

Missing proofs in Section 6

D.1

Proof of Lemma 6.5

Towards proving Lemma 6.5, we use the following lemma to express the inverse of a polynomial as
a sum of inverses of degree-1 polynomials.
Lemma D.1. Let p(z) = (z − λ1 ) . . . (z − λn ) where λj ’s are distinct. Then we have that
n

X tj
1
=
,
p(z)
z − λj

where

tj =

Q

i6=j (λj

−1
− λi )
.

(D.1)

j=1

Proof of Lemma D.1. By interpolating constant function at points λ1 , . . . , λn using Lagrange interpolating formula, we have that
n Q
X
(x − λi )
Q i6=j
1=
·1
(D.2)
i6=j (λj − λi )
j=1

Dividing p(z) on both sides we obtain equation (D.1).
The following lemma computes the Fourier transform of function 1/(z − λ).
32

Lemma D.2. Let m ∈ Z, and K be the unit circle in complex plane, and λ ∈ C inside the K. Then
we have that

Z
zm
2πiλm for m ≥ 0
dz =
0
o.w.
K z−λ
Proof of Lemma D.2. For m ≥ 0, since z m is a holomorphic function, by Cauchy’s integral formula,
we have that
Z
zm
dz = 2πiλm .
z
−
λ
K
For m < 0, by changing of variable y = z −1 we have that
Z
Z −m−1
zm
y
dz =
dy .
K z−λ
K 1 − λy
since |λy| = |λ| < 1, then we by Taylor expansion we have,
Z
K

y −m−1
dy =
1 − λy

Z

y −m−1

K

∞
X

!
(λy)k

dy .

k=0

Since the series λy is dominated by |λ|k which converges, we can switch the integral with the sum.
Note that y −m−1 is holomorphic for m < 0, and therefore we conclude that
Z −m−1
y
dy = 0 .
K 1 − λy
Now we are ready to prove Lemma 6.5.
Proof of Lemma 6.5. Let m = n + d. We compute the Fourier transform of z m /p(z). That is, we
write
∞
X
eimθ
βk eikθ .
=
p(eiθ )
k=−∞

where
1
βk =
2π

Z
0

2π

ei(m−k)θ
1
dθ =
iθ
2πi
p(e )

By Lemma D.1, we write

Z
K

z m−k−1
dz
p(z)

n

X tj
1
=
.
p(z)
z − λj
j=1

Then it follows that
n

1 X
βk =
tj
2πi

Z

j=1

Using Lemma D.2, we obtain that
 Pn
m−k−1
j=1 tj λj
βk =
0
We claim that

n
X
j=1

tj λn−1
= 1,
j

and

K

z m−k−1
dz
z − λj

if − ∞ ≤ k ≤ m − 1
o.w.
n
X
j=1

33

tj λsj = 0 ,

0 ≤ s < n − 1.

(D.3)

Indeed these can be obtained by writing out the lagrange interpolation for polynomial f (x) = xs
with s ≤ n − 1 and compare the leading coefficient. Therefore, we further simplify βk to
 Pn
m−k−1
if − ∞ < k < m − n

j=1 tj λj
βk =
(D.4)
1
if k = m − n

0
o.w.
P
Let h(z) = k≥0 βk z k . Then we have that h(z) is a polynomial with degree d = m − n and leading
term 1. Moreover, for our choice of d,
X
X
X
zm
− h(z) =
βk z k ≤
|βk | ≤ max |tj |(1 − λj )n
(1 − γ)d−k−1
j
p(z)
k<0

k<0

k<0

d

≤ Γ(1 − γ) /γ < ζ .

D.2

Proof of Theorem 6.9

Theorem 6.9 follows directly from a combination of Lemma D.3 and Lemma D.4 below. Lemma D.3
shows that the denominator of a function (under the stated assumptions) can be extended to a
polynomial that takes values in C + on unit circle. Lemma D.4 shows that it can be further extended
to another polynomial that takes values in C.
Lemma D.3. Suppose the roots of s are inside circle with radius α < 1, and Γ = Γ(s). If transfer
function G(z) = s(z)/p(z) satisfies that G(z) ∈ Cτ0 ,τ1 ,τ2 (or G(z) ∈ Cτ+0 ,τ1 ,τ2 ) for any z on unit circle,
√

nΓ·kpk

H2
1
then there exists u(z) of degree d = Oτ (max{( 1−α
log
, 0}) such that p(z)u(z)/z n+d ∈
1−α
+
0
n+d
∈ Cτ 0 ,τ 0 ,τ 0 respectively) for τ = Θτ (1) , where Oτ (·), Θτ (·) hide the
Cτ 0 ,τ 0 ,τ 0 (or p(z)u(z)/z
0 1 2
0 1 2
polynomial dependencies on τ0 , τ1 , τ2 .

Proof of Lemma D.3. By the fact that G(z) = s(z)/p(z) ∈ Cτ0 ,τ1 ,τ2 , we have that p(z)/s(z) ∈
Cτ00 ,τ10 ,τ20 for some τ 0 that polynomially depend on τ . Using Lemma 6.5, there exists u(z) of degree
d such that
z n+d
− u(z) ≤ ζ .
s(z)
where we set ζ  min{τ00 , τ10 }/τ20 · kpk−1
H∞ . Then we have that
p(z)u(z)/z n+d −

p(z)
≤ |p(z)|ζ  min{τ00 , τ10 }.
s(z)

(D.5)

It follows from equation (D.5) implies that that p(z)u(z)/z n+d ∈ Cτ000 ,τ100 ,τ200 , where τ 00 polynomially
depends on τ . The same proof still works when we replace C by C + .
Lemma D.4. Suppose p(z) of degree n and leading coefficient 1 satisfies that p(z) ∈ Cτ+0 ,τ1 ,τ2 for
any z on unit circle. Then there exists u(z) of degree d such that p(z)u(z)/z n+d ∈ Cτ00 ,τ10 ,τ20 for any
z on unit circle with d = Oτ (n) and τ00 , τ10 , τ20 = Θτ (1), where Oτ (·), Θτ (·) hide the dependencies on
τ0 , τ1 , τ2 .
p
Proof of Lemma D.4. We fix z on unit circle first. Let’s defined p(z)/z n be the square root
of p(z)/z n with principle value. Let’s write p(z)/z n = τ2 (1 + ( τp(z)
− 1)) and we take Taylor
P 2 z n

−1/2
−1/2
p(z)
p(z)
∞
1
−1/2 = τ
k . Note that since
expansion for √
=
τ
(1
+
(
−
1))
(
−
1)
n
n
2
2
k=0
τ
z
τ
z
n
2
2
p(z)/z

34

τ1 < |p(z)| < τ2 , we have that | τp(z)
n − 1| < 1 − τ1 /τ2 . Therefore truncating the Taylor series at
2z
k = Oτ (1) we obtain a polynomial a rational function h(z) of the form
h(z) =

p(z)
j≥0 ( τ2 z n

Pk

− 1)j ,

1
p(z)/z n

with precision ζ  min{τ0 , τ1 }/τ2 , that is, √ 1 n − h(z) ≤ ζ .
p(z)/z
p
p(z)h(z)
n
n
− p(z)/z ≤ ζ|p(z)/z | ≤ ζτ2 . Note that since p(z)/z n ∈
Therefore, we obtain that
zn
p
Cτ+0 ,τ1 ,τ2 , we have that p(z)/z n ∈ Cτ00 ,τ10 ,τ20 for some constants τ00 , τ10 , τ20 . Therefore p(z)h(z)
∈ Cτ00 ,τ10 ,τ20 .
zn
nk
Note that h(z) is not a polynomial yet. Let u(z) = z h(z) and then u(z) is polynomial of degree
at most nk and p(z)u(z)/z (n+1)k ∈ Cτ00 ,τ10 ,τ20 for any z on unit circle.
which approximates √

E

Back-propagation implementation

In this section we give a detailed implementation of using back-propagation to compute the gradient
of the loss function. The algorithm is for general MIMO case with the parameterization (7.3). To
obtain the SISO sub-case, simply take `in = `out = 1.
Algorithm 3 Back-propagation
Parameters: â ∈ Rn , Ĉ ∈ R`in ×n`out , and D̂ ∈ R`in ×`out . Let Â = MCC(â) = CC(â) ⊗ I`in and
B̂ = en ⊗ I`in .
Input: samples ((x(1) , y 1 ), . . . , x(N ) , y (N ) ) and projection set Bα .
for each sample (x(j) , y j ) = ((x1 , . . . , xT ), (y1 , . . . , yT )) do
Feed-forward pass:
h0 = 0 ∈ Rn`in .
for k = 1 to T
ĥk ← Âhk−1 + B̂xk , ŷt ← Ĉhk + D̂xk and ĥk ← Âhk−1 + B̂xk .
end for
Back-propagation:
∆hT +1 ← 0, GA ← 0, GC ← 0. GD ← 0
T1 ← T /4
for k = T to 1
if k > T1 , ∆yk ← ŷk − yk , o.w. ∆yk ← 0. Let ∆hk ← Ĉ > ∆yk + Â> ∆hk+1 .
1
1
update GC ← GC + T −T
∆yk ĥk , GA ← GA − T −T
B > ∆hk ĥ>
k−1 , and GD ← GD +
1
1
1
T −T1 ∆yk xk .
end for
Gradient update: Â ← Â − η · GA , Ĉ ← Ĉ − η · GC , D̂ ← D̂ − η · GD .
Projection step: Obtain â from Â and set â ← ΠB (â), and Â = MCC(â)
end for

F

Projection to the set Bα

In order to have a fast projection algorithm to the convex set Bα , we consider a grid GM of size M
over the circle with radius α. We will show that M = Oτ (n) will be enough to approximate the set
Bα in the sense that projecting to the approximating set suffices for the convergence.
35

0
Let Bα,τ
= {a : pa (z)/z n ∈ Cτ0 ,τ1 ,τ2 , ∀z ∈ GM } and Bα,τ0 ,τ1 ,τ2 = {a : pa (z)/z n ∈ Cτ0 ,τ1 ,τ2 , ∀|z| =
0 ,τ1 ,τ2
α}. Here Cτ0 ,τ1 ,τ2 is defined the same as before though we used the subscript to emphasize the dependency on τi ’s,

Cτ0 ,τ1 ,τ2 = {z : <z ≥ (1 + τ0 )|=z|} ∩ {z : τ1 < <z < τ2 } .

(F.1)

0
We will first show that with M = Oτ (n), we can make Bα,τ
to be sandwiched within to
1 ,τ2 ,τ3
two sets Bα,τ0 ,τ1 ,τ2 and Bα,τ00 ,τ10 ,τ20 .

Lemma F.1. For any τ0 > τ00 , τ1 > τ10 , τ2 < τ20 , we have that for M = Oτ (n), there exists κ0 , κ1 , κ2
0
⊂ Bα,τ00 ,τ10 ,τ20
that polynomially depend on τi , τi0 ’s such that Bα,τ0 ,τ1 ,τ2 ⊂ Bα,κ
0 ,κ1 ,κ2
Before proving the lemma, we demonstrate how to use the lemma in our algorithm: We will
pick τ00 = τ0 /2, τ10 = τ1 /2 and τ20 = 2τ2 , and find κi ’s guaranteed in the lemma above. Then
0
we use Bα,κ
as the projection set in the algorithm (instead of Bα,τ0 ,τ1 ,τ2 )). First of all, the
0 ,κ1 ,κ2
0
0
ground-truth solution Θ is in the set Bα,κ
. Moreover, since Bα,κ
⊂ Bα,τ00 ,τ10 ,τ20 , we will
0 ,κ1 ,κ2
0 ,κ1 ,κ2
b will remain in the set Bα,τ 0 ,τ 0 ,τ 0 and therefore the quasi-convexity of
guarantee that the iterates Θ
0

1

2

the objective function still holds8 .
0
Note that the set Bα,κ
contains O(n) linear constraints and therefore we can use linear
0 ,κ1 ,κ2
programming to solve the projection problem. Moreover, since the points on the grid forms a Fourier
basis and therefore Fast Fourier transform can be potentially used to speed up the projection.
Finally, we will prove Lemma F.1. We need S. Bernstein’s inequality for polynomials.
Theorem F.2 (Bernstein’s inequality, see, for example, [Sch41]). Let p(z) be any polynomial of
degree n with complex coefficients. Then,
sup |p0 (z)| ≤ n sup |p(z)|.
|z|≤1

|z|≤1

We will use the following corollary of Bernstein’s inequality.
Corollary F.3. Let p(z) be any polynomial of degree n with complex coefficients. Then, for m =
20n,
sup |p0 (z)| ≤ 2n sup |p(e2ikπ/m )|.
|z|≤1

k∈[m]

Proof. For simplicity let τ = supk∈[m] |p(e2ikπ/m )|, and let τ 0 = supk∈[m] |p(e2ikπ/m )|. If τ 0 ≤ 2τ
then we are done by Bernstein’s inequality. Now let’s assume that τ 0 > 2τ . Suppose p(z) = τ 0 .
Then there exists k such that |z − e2πik/m | ≤ 4/m and |p(e2πik/m )| ≤ τ . Therefore by Cauchy
mean-value theorem we have that there exists ξ that lies between z and e2πik/m such that p0 (ξ) ≥
m(τ 0 − τ )/4 ≥ 1.1nτ 0 , which contradicts Bernstein’s inequality.
Lemma F.4. Suppose a polynomial of degree n satisfies that |p(w)| ≤ τ for every w = αe2iπk/m for
some m ≥ 20n. Then for every z with |z| = α there exists w = αe2iπk/m such that |p(z) − p(w)| ≤
O(nατ /m).
Proof. Let g(z) = p(αz) by a polynomial of degree at most n. Therefore we have g 0 (z) = αp(z).
Let w = αe2iπk/m such that |z − w| ≤ O(α/m). Then we have
1
|p(z) − p(w)| = |g(z/α) − p(w/α)| ≤ sup |g 0 (x)| · |z − w| (By Cauchy’s mean-value Theorem)
α
|x|≤1
≤ sup |p0 (x)| · |z − w| ≤ nτ |z − w| .

(Corallary F.3)

|x|≤1

≤ O(αnτ /m) .
8

with a slightly worse parameter up to constant factor since τi ’s are different from τi ’s up to constant factors

36

Now we are ready to prove Lemma F.1.
Proof of Lemma F.1. We choose κi = 12 (τi + τi0 ).The first inequality is trivial. We prove the second
0
−1 )z n .
one. Consider a such that a ∈ Bα,κ0 ,κ1 ,κ2 . We wil show that a ∈ Bα,τ
0 ,τ 0 ,τ 0 . Let qa (z) = p(z
0

1

2

By Lemma F.4, for every z with |z| = 1/α, we have that there exists w = α−1 e2πik/M for some
integer k such that |qa (z) − qa (w)| ≤ O(τ2 n/(αM )). Therefore let M = cn for sufficiently large c
(which depends on τi ’s), we have that for every z with |z| = 1/α, qa (z) ∈ Cτ00 ,τ10 ,τ20 . This completes
the proof.

37

