Gradient Descent Learns Linear Dynamical Systems

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 $x_{1},\dots,x_{T}$ to a corresponding sequence of observations $y_{1},\dots,y_{T}.$ 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 Hochreiter and Schmidhuber (1997), 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 likelihood estimation with gradient descent is a popular heuristic for dynamical system identification Ljung (1998). 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 Levine and Koltun (2013), location and mapping estimation in robotic systems Durrant-Whyte and Bailey (2006), and estimation of pose from video Rahimi et al. (2005).

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:

Here, $A,B,C,D$ are linear transformations with compatible dimensions and we denote by $\Theta=(A,B,C,D)$ the parameters of the system. The vector $h_{t}$ represents the hidden state of the model at time $t.$ Its dimension $n$ is called the order of the system. The stochastic noise variables $\{\xi_{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 and variance $\sigma^{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 $\rho(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,

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 $(\hat{A},\hat{B},\hat{C},\hat{D})$, the trained model is defined as:

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 $\hat{y}_{t}$ the $t$-th prediction of the trained model starting from the correct initial state that generated $y_{t}$, and using $\widehat{\Theta}$ as a short hand for $(\hat{A},\hat{B},\hat{C},\hat{D})$, we formally define population risk as:

Note that even though the prediction $\hat{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 $x_{t}$ and errors $\xi_{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 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.

2 Proper learning

Under our assumption, projected stochastic gradient descent, when given $N$ sample sequence of length $T$, returns parameters $\widehat{\Theta}$ with population risk

Here the expectation on LHS is with respect to the randomness of the algorithm. Recall that $f(\Theta)$ is the population risk of the optimal system, and $\sigma^{2}$ refers to the variance of the noise variables. We also assume that the inputs $x_{t}$ are drawn from a pairwise independent distribution with mean 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$. The dependency on the dimension $n$, on the other hand, is likely to be quite loose, and tighter bounds are left for future work.

The algorithm requires a (polynomial-time) projection step to a convex set at every iteration (formally defined in Section 4 and Algorithm 1). Computationally, it can be a bottleneck, although it is unlikely to be required in practice and may be an artifact of our analysis.

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 $n^{\prime}>n$ such that the characteristic polynomial of its transition matrix becomes $p\cdot p^{\prime}$ for some polynomial $p^{\prime}$ of order $n^{\prime}-n.$ This means that to apply our result we only need the polynomial $p\cdot p^{\prime}$ to satisfy our assumption. At this point, we can choose $p^{\prime}$ to be an approximation of the inverse $p^{-1}$. For sufficiently good approximation, the resulting polynomial $p\cdot p^{\prime}$ is close to $1$ and therefore satisfies our assumption. Such an approximation exists generically for $n^{\prime}=O(n)$ under mild non-degeneracy assumptions on the roots of $p$. In particular, any small random perturbation of the roots would suffice.

Under a mild non-degeneracy assumption, stochastic gradient descent returns parameters $\widehat{\Theta}$ corresponding to a system of order $n^{\prime}=O(n)$ with population risk

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 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.

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 Ljung (1998). 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 Vidyasagar and Karandikar (2008); Campi and Weyer (2002); Weyer and Campi (1999). 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 Shah et al. (2012). 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-\rho(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-\rho(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 $u_{t}$ and output sequence $y_{t}$ 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 Megretski (2008). 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. Söderström and Stoica (1982); Stoica and Söderström (1982, 1984) and later Bazanella et al. (2008); Eckhard and Bazanella (2011) prove the quasi-convexity of a cost function under a passivity condition in the context of system identification, but no sample complexity or global convergence proofs are provided.

6 Proof overview

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)\approx p(z)^{-1}$ so that $p(z)u(z)\approx 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.

7 Preliminaries

A SISO of order $n$ is in controllable canonical form if $A$ and $B$ have the following form

We will parameterize $\hat{A},\hat{B},\hat{C},\hat{D}$ accordingly. We will write $A=\textup{CC}(a)$ for brevity, where $a$ is used to denote the unknown last row $[-a_{n},\dots,-a_{1}]$ of matrix $A$. We will use $\hat{a}$ to denote the corresponding training variables for $a$. Since here $B$ is known, so $\hat{B}$ is no longer a trainable parameter, and is forced to be equal to $B$. Moreover, $C$ is a row vector and we use $[c_{1},\cdots,c_{n}]$ to denote its coordinates (and similarly for $\hat{C}$).

A SISO is controllable if and only if the matrix $[B\mid AB\mid A^{2}B\mid\cdots\mid A^{n-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 Heij et al. (2007). For vector $a=[a_{n},\dots,a_{1}]$, let $p_{a}(z)$ denote the polynomial

Gradient descent and quasi-convexity

It is known that under certain mild conditions (stochastic) gradient descent converges even on non-convex functions to local minimum Ge et al. (2015); Lee et al. (2016). 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 Hazan et al. (2015), which ensures that any point with vanishing gradient is the optimal solution . Roughly speaking, the condition says that at any point $\theta$ the negative of the gradient $-\nabla f(\theta)$ should be positively correlated with direction $\theta^{*}-\theta$ pointing towards the optimum. Our condition is slightly weaker than that in Hazan et al. (2015) since we only require quasi-convexity and smoothness with respect to the optimum, and this (simple) extension will be necessary for our analysis.

We say an objective function $f$ is $\tau$-weakly-quasi-convex ($\tau$-WQC) over a domain $\mathcal{B}$ with respect to global minimum $\theta^{*}$ if there is a positive constant $\tau>0$ such that for all $\theta\in\mathcal{B}$,

We further say $f$ is $\Gamma$-weakly-smooth if for for any point $\theta$, $\|\nabla f(\theta)\|^{2}\leq\Gamma(f(\theta)-f(\theta^{*}))$.

Projected stochastic gradient descent over some closed convex set $\mathcal{B}$ with learning rate $\eta>0$ refers to the following algorithm in which $\Pi_{\mathcal{B}}$ denotes the Euclidean projection onto $\mathcal{B}$:

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

Suppose the objective function $f$ is $\tau$-weakly-quasi-convex and $\Gamma$-weakly-smooth, and $\mathfrak{r}(\cdot)$ is an unbiased estimator for $\nabla f(\theta)$ with variance $V$. Moreover, suppose the global minimum $\theta^{*}$ belongs to $\mathcal{B}$, and the initial point $\theta_{0}$ satisfies $\|\theta_{0}-\theta^{*}\|\leq R$. Then projected gradient descent (2.2) with a proper learning rate returns $\theta_{K}$ in $K$ iterations with expected error

We defer the proof of Proposition 2.3 to Appendix A which is a simple variation of the standard convergence analysis of stochastic gradient descent (see, for example, Bottou (1998)). Finally, we note that the sum of two quasi-convex functions may no longer be quasi-convex. However, if a sequence functions is $\tau$-WQC with respect to a common point $\theta^{*}$, then their sum is also $\tau$-WQC. This follows from the linearity of gradient operation.

Suppose functions $f_{1},\dots,f_{n}$ are individually $\tau$-weakly-quasi-convex in $\mathcal{B}$ with respect to a common global minimum $\theta^{*}$ , then for non-negative $w_{1},\dots,w_{n}$ the linear combination $f=\sum_{i=1}^{n}w_{i}f_{i}$ is also $\tau$-weakly-quasi-convex with respect to $\theta^{*}$ in $\mathcal{B}$.

Population risk in frequency domain

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

The leading term $\|D-\hat{D}\|^{2}$ is convex in $\hat{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.5, and we restrict our attention to the remaining terms.

We now use basic concepts from control theory (see Heij et al. (2007); Hespanha (2009) for more background) to express the idealized risk (3.2) in Fourier domain. The transfer function of the linear system is

Note that $G(z)$ is a rational function over the complex numbers of degree $n$ and hence we can find polynomials $s(z)$ and $p(z)$ such that $G(z)=\frac{s(z)}{p(z)}\,,$ with the convention that the leading 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

The sequence $r=(r_{0},r_{1},\ldots,r_{t},\ldots)=(CB,CAB,\ldots,CA^{t}B,\ldots)$ 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 $\widehat{G}(z)=\hat{C}(zI-\hat{A})^{-1}B=\hat{s}(z)/\hat{p}(z)\,.$ The idealized risk (3.2) is only a function of the impulse response $\hat{r}$ of the learned system, and therefore it is only a function of $\widehat{G}(z)$.

Recall that $C=[c_{1},\dots,c_{n}]$ is defined in Section 1.7. For future reference, we note that by some elementary calculation (see Lemma B.1), we have

which implies that $s(z)=c_{1}+\dots+c_{n}z^{n-1}$ and $p(z)=z^{n}+a_{1}z^{n-1}+\dots+a_{n}$.

With these definitions in mind, we are ready to express idealized risk in Fourier domain.

Suppose $p_{\hat{a}}(z)$ has all its roots inside unit circle, then the idealized risk $g(\hat{a},\hat{C})$ can be written in the Fourier domain as

Proof Note that $G(e^{i\theta})$ is the Fourier transform of the sequence $\{r_{k}\}$ and so is $\widehat{G}(e^{i\theta})$ the Fourier transformThe Fourier transform exists since $\|r_{k}\|^{2}=\|\hat{C}\hat{A}^{k}\hat{B}\|^{2}\leq\|\hat{C}\|\|\hat{A}^{k}\|\|\hat{B}\|\leq c\rho(\hat{A})^{k}$ where $c$ doesn’t depend on $k$ and $\rho(\hat{A})<1$. of $\hat{r}_{k}$. Therefore by Parseval’ Theorem, we have that

Now that we have a convenient expression for risk in Fourier domain, we can prove that the idealized risk $g(\hat{A},\hat{C})$ is weakly-quasi-convex when $\hat{a}$ is not so far from the true $a$ in the sense that $p_{a}(z)$ and $\hat{p}_{a}(z)$ have an angle less than $\pi/2$ for every $z$ on the unit circle. We will use the convention that $a$ and $\hat{a}$ refer to the parameters that specify $A$ and $\hat{A},$ respectively.

For $\tau>0$ and every $\hat{C}$, the idealized risk $g(\hat{A},\hat{C})$ is $\tau$-weakly-quasi-convex over the domain

Proof We first analyze a single term $h=|\widehat{G}(z)-G(z)|^{2}$. Recall that $\widehat{G}(z)=\hat{s}(z)/\hat{p}(z)$ where $\hat{p}(z)=p_{\hat{a}}(z)=z^{n}+\hat{a}_{1}z^{n-1}+\dots+\hat{a}_{n}$. Note that $z$ is fixed and $h$ is a function of $\hat{a}$ and $\hat{C}$. Then it is straightforward to see that

Since $\hat{s}(z)$ and $\hat{p}(z)$ are linear in $\hat{C}$ and $\hat{a}$ respectively, by chain rule we have that

Plugging the formulas (3.7) and (3.8) for ${\frac{\partial h}{\partial\hat{s}(z)}}$ and ${\frac{\partial h}{\partial\hat{p}(z)}}$ into the equation above, we obtain that

Hence $h=|\widehat{G}(z)-G(z)|^{2}$ is $\tau$-weakly-quasi-convex with $\tau=2\min_{|z|=1}\Re\left\{\frac{p(z)}{\hat{p}(z)}\right\}$. This implies our claim, since by Proposition 3.2, the idealized risk $g$ is convex combination of functions of the form $|\widehat{G}(z)-G(z)|^{2}$ for $|z|=1$. Moreover, Proposition 2.5 shows that convex combination preserves weak quasi-convexity.

For future reference, we also prove that the idealized risk is $O(n^{2}/\tau_{1}^{4})$-weakly smooth.

The idealized risk $g(\hat{A},\hat{C})$ is $\Gamma$-weakly smooth with $\Gamma=O(n^{2}/\tau_{1}^{4})$.

Proof By equation (3.8) and the chain rule we get that

therefore we can bound the norm of the gradient by

Similarly, we could show that $\left\|{\frac{\partial g}{\partial\hat{a}}}\right\|^{2}\leq O(n^{2}/\tau_{1}^{2})\cdot g(\hat{A},\hat{C})$.

2 Justifying idealized risk

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

Assume that $\xi_{t}$ and $x_{t}$ are drawn i.i.d. from an arbitrary distribution with mean and variance $1.$ Then the population risk $f(\widehat{\Theta})$ can be written as,

The idealized risk is upper bound of the population risk $f(\widehat{\Theta})$ according to equation (3.1) and (3.9). We don’t have to quantify the gap between them because later in Algorithm 1, we will directly optimize the idealized risk by constructing an estimator of its gradient, and thus the optimization will guarantee a bounded idealized risk which translates to a bounded population risk. See Section 5 for details.

Proof [Proof of Lemma 3.5] Under the distributional assumptions on $\xi_{t}$ and $x_{t},$ we can calculate the objective functions above analytically. We write out $y_{t},\hat{y}_{t}$ in terms of the inputs,

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

Recall that under the controllable canonical form (1.4), $B=e_{n}$ is known and therefore $\hat{B}=B$ is no longer a variable. Then the expected objective function (3.11) simplifies to

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.

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 $\mathcal{C}$ in the complex plane,

We say a polynomial $p(z)$ is $\alpha$-acquiescent if $\{p(z)/z^{n}\colon|z|=\alpha\}\subseteq\mathcal{C}.$ A linear system with transfer function $G(z)=s(z)/p(z)$ is $\alpha$-acquiescent if the denominator $p(z)$ is.

Suppose $a\in\mathcal{B}_{\alpha}$, then the roots of polynomial $p_{a}(z)$ have magnitudes bounded by $\alpha$. Therefore the controllable canonical form $A=\textup{CC}(a)$ defined by $a$ has spectral radius $\rho(A)<\alpha$.

Proof Define holomorphic function $f(z)=z^{n}$ and $g(z)=p_{a}(z)=z^{n}+a_{1}z^{n-1}+\dots+a_{n}$. We apply the symmetric form of Rouche’s theorem Estermann (1962) on the circle $\mathcal{K}=\{z:|z|=\alpha\}$. For any point $z$ on $\mathcal{K}$, we have that $|f(z)|=\alpha^{n}$, and that $|f(z)-g(z)|=\alpha^{n}\cdot|1-p_{a}(z)/z^{n}|$. Since $a\in\mathcal{B}_{\alpha}$, we have that $p_{a}(z)/z^{n}\in\mathcal{C}$ for any $z$ with $|z|=\alpha$. Observe that for any $c\in\mathcal{C}$ we have that $|1-c|<1+|c|$, therefore we have that

Hence, using Rouche’s Theorem, we conclude that $f$ and $g$ have same number of roots inside circle $\mathcal{K}$. Note that function $f=z^{n}$ has exactly $n$ roots in $\mathcal{K}$ and therefore $g$ have all its $n$ roots inside circle $\mathcal{K}$.

The following lemma establishes the fact that $\mathcal{B}_{\alpha}$ is a monotone family of sets in $\alpha$. The proof follows from the maximum modulo principle of the harmonic functions $\Re(z^{n}p(1/z))$ and $\Im(z^{n}p(1/z))$. We defer the short proof to Section C.1. We remark that there are larger convex sets than $\mathcal{B}_{\alpha}$ 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 $\mathcal{B}_{\alpha}.$

For any $0<\alpha<\beta$, we have that $\mathcal{B}_{\alpha}\subset\mathcal{B}_{\beta}$.

Our next lemma entails that acquiescent systems have well behaved impulse responses.

Suppose $a\in\mathcal{B}_{\alpha}$ for some $\alpha\leq 1$. Then the companion matrix $A=\textup{CC}(a)$ satisfies

Moreover, it holds that for any $k\geq 0$,

Let $f_{\lambda}=\sum_{k=0}^{\infty}e^{i\lambda k}\alpha^{-k}A^{k}B$ be the Fourier transform of the series $\alpha^{-k}A^{k}B$. Then using Parseval’s Theorem, we have

where at the last step we used the fact that $(I-wA)^{-1}B=\frac{1}{p_{a}(w^{-1})}[w^{-1},w^{-2}\dots,z^{-n}]^{\top}$ (see Lemma B.1), and that $\alpha\leq 1$. Since $a\in\mathcal{B}_{\alpha}$, we have that $|q_{a}(\alpha^{-1}e^{i\lambda})|\geq\tau_{1}$, and therefore $p_{a}(\alpha e^{-i\lambda})=\alpha^{n}e^{-in\lambda}q(e^{i\lambda}/\alpha)$ has magnitude at least $\tau_{1}\alpha^{n}$. Plugging in this into equation (4.4), we conclude that

Finally we establish the bound for $\|A^{k}B\|^{2}$. By Lemma 4.3, we have $\mathcal{B}_{\alpha}\subset\mathcal{B}_{1}$ for $\alpha\leq 1$. Therefore we can pick $\alpha=1$ in equation (4.3) and it still holds. That is, we have that

This also implies that $\|A^{k}B\|^{2}\leq 2\pi n/\tau_{1}^{2}$.

Learning acquiescent systems

Next we show that we can learn acquiescent systems.

Suppose the true system $\Theta$ is $\alpha$-acquiescent and satisfies $\|C\|\leq 1$. Then with $N$ samples of length $T\geq\Omega(n+1/(1-\alpha))$, stochastic gradient descent (Algorithm 1) with projection set $\mathcal{B}_{\alpha}$ returns parameters $\widehat{\Theta}=(\hat{A},\hat{B},\hat{C},\hat{D})$ with population risk

where $O(\cdot)$-notation hides polynomial dependencies on $1/(1-\alpha)$, $1/\tau_{0},1/\tau_{1},\tau_{2}$, and $R=\|a\|$. The expectation is taken over the randomness of the algorithms and the examples.

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). 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/\sqrt{N}$ (for relatively small $N$).

Computational complexity: Step 2 in each iteration of the algorithm takes $O(Tn)$ arithmetic operations, and the projection step takes $O(n^{3.5})$ time to solve an linear programming problem. The project step is unlikely to be required in practice and may be an artifact of our analysis.

We can further balance the variance of the estimator with the number of samples by breaking each long sequence of length $T$ into $\Theta(T/n)$ short sequences of length $\Theta(n)$, and then run back-propagation (1) on these $TN/n$ shorter sequences. This leads us to the following bound which gives the right dependency in $T$ and $N$ as we expected: $TN$ should be counted as the true number of samples for the sequence-to-sequence model.

Under the assumption of Theorem 5.1, Algorithm 2 returns parameters $\widehat{\Theta}$ with population risk

where $O(\cdot)$-notation hides polynomial dependencies on $1/(1-\alpha)$, $1/\tau_{0},1/\tau_{1},\tau_{2}$, and $R=\|a\|.$

We remark the the gradient computation procedure takes time linear in $Tn$ 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^{\prime}\gg\max\{n,1/{(1-\alpha)}\}$. By the explicit formula for $f(\widehat{\Theta})$ (Lemma 3.5) and the fact that $\|CA^{k}B\|$ decays exponentially for $k\gg n$ (Lemma 4.4), we can bound the population risk on sequences of different lengths. Concretely, let $f_{T^{\prime}}(\widehat{\Theta})$ denote the population risk on sequence of length $T^{\prime}$, we have for all $T^{\prime}\gg\max\{n,1/{(1-\alpha)}\}$,

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.

We could hope to achieve linear convergence by showing that the empirical risk also satisfies the weakly-quasi-convexity. Then, we can re-use the samples and hope to use strong optimization tools (such as SVRG) to achieve the linear convergence. This is beyond the scope of this paper and left to future work.

Our proof of Theorem 5.1 simply consists of three parts: a) showing the idealized risk is weakly quasi-convex in the convex set $\mathcal{B}_{\alpha}$ (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).

Proof [Proof of Lemma 5.3] It suffices to show that for all $\hat{a},a\in\mathcal{B}_{\alpha}$, it satisfies $\hat{a}\in\mathcal{N}_{\tau}(a)$ for $\tau=\Omega(\tau_{0}\tau_{1}/\tau_{2})$. Indeed, by the monotonicity of the family of sets $\mathcal{B}_{\alpha}$ (Lemma 4.3), we have that $\hat{a},a\in\mathcal{B}_{1}$, which by definition means for every $z$ on unit circle, $p_{a}(z)/z^{n},p_{\hat{a}}(z)/z^{n}\in\mathcal{C}$. By definition of $\mathcal{C}$, for any point $w,\hat{w}\in\mathcal{C}$, the angle $\phi$ between $w$ and $\hat{w}$ is at most $\pi-\Omega(\tau_{0})$ and ratio of the magnitude is at least $\tau_{1}/\tau_{2}$, which implies that $\Re(w/\hat{w})=|w|/|\hat{w}|\cdot\cos(\phi)\geq\Omega(\tau_{0}\tau_{1}/\tau_{2})$. Therefore $\Re(p_{a}(z)/p_{\hat{a}}(z))\geq\Omega(\tau_{0}\tau_{1}/\tau_{2})$, and we conclude that $\hat{a}\in\mathcal{N}_{\tau}(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 $h_{0}$, 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 $h_{0}$ is unknown (or even to be learned). Moreover, the ability of handling unknown $h_{0}$ allows us to break a 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 $h_{0}$. We get around it by simply ignoring first $T_{1}=\Omega(T)$ outputs of the system and setting the corresponding errors to 0. Since the influence of $h_{0}$ to any outputs later than time $k\geq T_{1}\gg\max\{n,1/{(1-\alpha)}\}$ is inverse exponentially small, we could safely assume $h_{0}=0$ when the error earlier than time $T_{1}$ 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 Sutskever et al. (2014), 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.

Under the assumption of Theorem 5.1, suppose $\hat{a},a\in\mathcal{B}_{\alpha}$. Then in Algorithm 1, at each iteration, $G_{A},G_{C}$ are unbiased estimators of the gradient of the idealized risk (3.2) in the sense that:

Finally, we control the variance of the gradient estimator.

The (almost) unbiased estimator $(G_{A},G_{C})$ of the gradient of $g(\hat{A},\hat{C})$ has variance bounded by

where $\Lambda=O(\max\{n,1/{(1-\alpha)}\log 1/{(1-\alpha)}\})$.

We bound directly the variance instead. It’s tedious but 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 $\hat{A},\hat{C}$ which consists of sums of terms involving $(\hat{C}\hat{A}^{k}B-CA^{k}B)$ and $\hat{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 [Proof of Theorem 5.1] We consider $g^{\prime}(\hat{A},\hat{C},\hat{D})=(\hat{D}-D)^{2}+g(\hat{A},\hat{C})$, an extended version of the idealized risk which takes the contribution of $\hat{D}$ into account. By Lemma 5.4 we have that Algorithm 1 computes $G_{A},G_{C}$ which are almost unbiased estimators of the gradients of $g^{\prime}$ up to negligible error $\exp(-\Omega({(1-\alpha)}T))$, and by Lemma C.2 we have $G_{D}$ is an unbiased estimator of $g^{\prime}$ with respect to $\hat{D}$. Moreover by Lemma 5.5, these unbiased estimator has total variance $V=\frac{O\left(n^{5}+\sigma^{2}n^{3}\right)}{T}$ where $O(\cdot)$ hides dependency on $\tau_{1}$ and ${(1-\alpha)}$. Applying Proposition 2.3 (which only requires an unbiased estimator of the gradient of $g^{\prime}$), we obtain that after $T$ iterations, we converge to a point with $g^{\prime}(\hat{a},\hat{C},\hat{D})\leq O\left(\frac{n^{2}}{N}+\sqrt{\frac{n^{5}+\sigma^{2}n^{3}}{TN}}\right)$. Then, by Lemma 3.5 we have $f(\widehat{\Theta})\leq g^{\prime}(\hat{a},\hat{C},\hat{D})+\sigma^{2}=g^{\prime}(\hat{a},\hat{C},\hat{D})+f(\Theta)\leq O\left(\frac{n^{2}}{N}+\sqrt{\frac{n^{5}+\sigma^{2}n^{3}}{TN}}\right)+f(\Theta)$ which completes the proof.

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 invariance of a system $G(z)=\frac{s(z)}{p(z)}$. For an arbitrary polynomial $u(z)$ of leading coefficient 1, we can write $G(z)$ as

Our discussion motivates the following definition.

A polynomial $p(z)$ of degree $n$ is $\alpha$-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 $\alpha$-acquiescent.

For a transfer function $G(z)$, we define it’s $\mathcal{H}_{2}$ norm as

We assume (with loss of generality) that the true transfer function $G(z)$ has bounded $\mathcal{H}_{2}$ norm, that is, $\|G\|_{\mathcal{H}_{2}}\leq 1$. This can be achieve by a rescalingIn fact, this is a natural scaling that makes comparing error easier. Recall that the population risk is essentially $\|\hat{G}-G\|_{\mathcal{H}_{2}}$, therefore rescaling $C$ so that $\|G\|_{\mathcal{H}_{2}}=1$ implies that when error $\ll 1$ we achieve non-trivial performance. of the matrix $C$.

Suppose the true system has transfer function $G(z)=s(z)/p(z)$ with a characteristic function $p(z)$ that is $\alpha$-acquiescent by extension of degree $d$, and $\|G\|_{\mathcal{H}_{2}}\leq 1$, then projected stochastic gradient descent with $m=n+d$ states (that is, Algorithm 2 with $m$ states) returns a system $\widehat{\Theta}$ with population risk

where the $O(\cdot)$ notation hides polynomial dependencies on $\tau_{0},\tau_{1},\tau_{2},1/(1-\alpha)$.

The theorem follows directly from Corollary 5.2 (with some additional care about the scaling.

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 $\alpha$-acquiescent but is $\alpha$-acquiescent by extension with small degree $d$.

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.

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)$.

Consider $G(z)=\frac{s(z)}{p(z)}:=\frac{z^{n}-0.8^{-n}}{(z-0.1)(z^{n}-0.9^{-n})}$ and $G^{\prime}(z)=\frac{s^{\prime}(z)}{p^{\prime}(z)}:=\frac{1}{z-0.1}$. Clearly these are the minimum representations of the $G(z)$ and $G^{\prime}(z)$, which also both satisfy acquiescence. On the one hand, the characteristic polynomial $p(z)$ and $p^{\prime}(z)$ are very different. On the other hand, the transfer functions $G(z)$ and $G^{\prime}(z)$ have almost the same values on unit circle up to exponentially small error,

Moreover, the transfer functions $G(z)$ and $\hat{G}(z)$ are on the order of $\Theta(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.

2 Power of improper learning in various cases

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

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 $(\alpha+1)/2$-$(\alpha+1)/2$-acquiescent by extension of degree $3$.

Let $n$ be a large enough integer and $\alpha$ be a constant. Let $J=\{1,n-1,n\}$ and $\omega=e^{2\pi i/n}$, and then define $p(z)=z^{3}\prod_{j\in[n],j\notin J}(z-\alpha\omega^{j})$. Therefore we have that

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

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.

For a polynomial $h(z)$ of degree $n$ with roots $\lambda_{1},\dots,\lambda_{n}$ inside unit circle, define the quantity $\Gamma(\cdot)$ of the polynomial $h$ as:

Suppose $p(z)$ is a polynomial of degree $n$ with distinct roots inside circle with radius $\alpha$. Let $\Gamma=\Gamma(p)$, then $p(z)$ is $\alpha$-acquiescent by extension of degree $d$ $~{}=O(\max\{(1-\alpha)^{-1}\log(\sqrt{n}\Gamma\cdot\|p\|_{\mathcal{H}_{2}}),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 $\mathcal{C}$. We believe the following lemma should be known though for completeness we provide the proof in Section D.

Suppose $p(z)$ is a polynomial of degree $n$ and leading coefficient 1 with distinct roots inside circle with radius $\alpha$, and $\Gamma=\Gamma(p)$. Then for $d=O(\max\{(\frac{1}{1-\alpha}\log\frac{\Gamma}{(1-\alpha)\zeta},0\})$, there exists a polynomial $h(z)$ of degree $d$ and leading coefficient 1 such that for all $z$ on unit circle,

Proof [Proof of Lemma 6.4] Let $\gamma=1-\alpha$. Using Lemma 6.5 with $\zeta=0.5\|p\|_{\mathcal{H}_{\infty}}^{-1}$, we have that there exists polynomial $u$ of degree $d=O(\max\{\frac{1}{1-\alpha}\log(\Gamma\|p\|_{\mathcal{H}_{\infty}}),0\})$ such that

Therefore $p(z)u(z)/z^{n+d}\in\mathcal{C}_{\tau_{0},\tau_{1},\tau_{2}}$ for constant $\tau_{0},\tau_{1},\tau_{2}$. Finally noting that for degree $n$ polynomial we have $\|h\|_{\mathcal{H}_{\infty}}\leq\sqrt{n}\cdot\|h\|_{\mathcal{H}_{2}}$, which completes the proof.

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 $\lambda_{1},\dots,\lambda_{n}$ uniformly randomly on circle with radius $\alpha<1$, and take $\lambda_{i},\bar{\lambda}_{i}$ for $i=1,\dots,n$ as the roots of $p(z)$. That is, $p(z)=(z-\lambda_{1})(z-\bar{\lambda}_{1})\dots(z-\lambda_{n})(z-\bar{\lambda}_{n})$. We show that with good probability (over the randomness of $\lambda_{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.

Towards proving Theorem 6.6, we need the following lemma about the expected distance of two random points with radius $\rho$ and $r$ in log-space.

Proof When $r\neq\rho$, let $N$ be an integer and $\omega=e^{2i\pi/N}$. Then we have that

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.

A SISO linear system is strict-input passive if and only if for some $\tau_{0}>0$ and any $z$ on unit circle, $\Re(G(z))\geq\tau_{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 $\tau_{0},\tau_{1},\tau_{2}$, define

We say a transfer function $G(z)=s(z)/p(z)$ is $(\tau_{0},\tau_{1},\tau_{2})$-strict input passive if for any $z$ on unit circle we have $G(z)\in\mathcal{C}^{+}_{\tau_{0},\tau_{1},\tau_{2}}$. Note that for small constant $\tau_{0},\tau_{1}$ and large constant $\tau_{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 $\Gamma(s)\leq\exp(O(n))$.

Suppose $G(z)=s(z)/p(z)$ is $(\tau_{0},\tau_{1},\tau_{2})$-strict-input passive. Moreover, suppose the roots of $s(z)$ have magnitudes inside circle with radius $\alpha$ and $\Gamma=\Gamma(s)\leq\exp(O(n))$ and $\|p\|_{\mathcal{H}_{2}}\leq\exp(O(n))$. Then $p(z)$ is $\alpha$-acquiescent by extension of degree $d=O_{\tau,\alpha}(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)\in\mathcal{C}_{\tau_{0},\tau_{1},\tau_{2}}$ for every $z$ on unit circle. Then $p(z)$ is $\alpha$-acquiescent by extension of degree $d=O_{\tau,\alpha}(n)$.

The proof of Theorem 6.9 is similar in spirit to that of Lemma 6.4, and is deferred to Section D.

3 Improper learning using linear regression

In this subsection, we show that under stronger assumption than $\alpha$-acquiescent by extension, we can improperly learn a linear dynamical system with linear regression, up to some fixed bias.

A polynomial $p(z)$ of degree $n$ is extremely-acquiescent by extension of degree $d$ with bias $\varepsilon$ if there exists a polynomial $u(z)$ of degree $d$ and leading coefficient 1 such that for all $z$ on unit circle,

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}\in\mathcal{C}$. Therefore, equation (6.4) above is a much stronger requirement than acquiescence by extension.We need $(1-\delta)$-acquiescence by extension in previous subsections for small $\delta>0$, though this is merely additional technicality needed for the sample complexity. We ignore this difference between $1-\delta$-acquiescence and $1$-acquiescence and for the purpose of this subsection

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)$,

Let’s define the risk function in the frequency domain as,

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 $\hat{S}(z)$) commute with scalar $p(z)$ and $\hat{p}(z)$, and that $\hat{S}(z),\hat{p}(z)$ are linear in $\hat{a},\hat{C}$.

The risk function $g(\hat{a},\hat{C})$ defined in (7.2) is $\tau$-weakly-quasi-convex in the domain

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 Brunovsky (1970),

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

Suppose transfer function $G(z)$ of a MIMO system takes form (7.1), and has norm $\|G\|_{\mathcal{H}_{2}}\leq 1$. If the common denominator $p(z)$ is $\alpha$-acquiescent by extension of degree $d$ then projected stochastic gradient descent over the state space representation (7.3) will return $\widehat{\Theta}$ with risk

Simulations

In this section, we provide proof-of-concepts experiments on synthetic data. We will demonstrate that

plain SGD tends to blow up even with relatively small learning rate, especially on hard instances

SGD with our projection step converges with reasonably large learning rate, and with over-parameterization the final error is competitive

SGD with gradient clipping has the strongest performance in terms both of the convergence speed and the final error

Here gradient clipping refers to the technique of using a normalized gradient instead of the true gradient. Specifically, for some positive hyper parameter $B$, we follow the approximate gradient

This method is commonly applied in training recurrent neural networks Pascanu et al. (2013).

Bullet 1) suggests that stability is indeed a real concern. Bullet 2) corroborates our theoretical study. Finding 3) suggests the instability of SGD partly arises from the noise in the batches, and such noise is reduced by the gradient clipping. Our experiments suggest that the landscape of the objective function may be even nicer than what is predicted by our theoretical development. It remains possible that the objective has no non-global local minima, possibly even outside the convex set to which our algorithm projects.

We generate the true system with state dimension $d=20$ by randomly picking the conjugate pairs of roots of the characteristic polynomial inside the circle with radius $\rho=0.95$ and randomly generating the vector $C$ from standard normal distribution. The distribution of the norm of the impulse response $r$ (defined in Section 3) of such systems has a heavy-tail. When the norm of $r$ is several magnitudes larger than the median it’s difficult to learn the system. Thus we select systems with reasonable $\|r\|$ for experiments, and we observe that the difficulty of learning increases as $\|r\|$ increases. The inputs of the dynamical model are generated from standard normal distribution with length $T=500$. We note that we generate new fresh inputs and outputs at every iterations and therefore the training loss is equal to the test loss (in expectation.) We use initial learning rate 0.01 in the projected gradient descent and SGD with gradient clipping. We use batch size 100 for all experiments, and decay the learning rate at 200K and 250K iteration by a factor of 10 in all experiments.

Acknowledgments

We thank Amir Globerson, Alexandre Megretski, Pablo A. Parrilo, Yoram Singer, Peter Stoica, 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

A 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 [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 descentAlthough we used weak smoothness to get a slightly better bound. We conditioned on $\theta_{k}$, and have that

where the first inequality uses weakly-quasi-convex and the rest of lines are simply algebraic manipulations. Since $\theta_{k+1}$ is the projection of $w_{k+1}$ to $\mathcal{B}$ and $\theta^{*}$ belongs to $\mathcal{B}$, we have $\|w_{k+1}-\theta^{*}\|\geq\|\theta_{k+1}-\theta^{*}\|$. Together with (A.1), and

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

where we use the assumption that $\|\theta_{0}-\theta^{*}\|\leq R$. Suppose $K\geq\frac{4R^{2}\Gamma^{2}}{V\tau^{2}}$, then we take $\eta=\frac{R}{\sqrt{VK}}$. Therefore we have that $\tau-\eta\Gamma\geq\tau/2$ and therefore

On the other hand, if $K\leq\frac{4R^{2}\Gamma^{2}}{V\tau^{2}}$, we pick $\eta=\frac{\tau}{2\Gamma}$ and obtain that

Therefore using equation (A.3) and (A.2) we obtain that when choosing $\eta$ properly according to $K$ as above,

B Toolbox

where $p_{a}(x)=x^{n}+a_{1}x^{n-1}+\dots+a_{n}$ is the characteristic polynomial of $A$.

Proof let $v=(I-wA)^{-1}B$ then we have $(I-wA)v=B$. Note that $B=e_{n}$, and $I-wA$ is of the form

Therefore we obtain that $v_{k}=wv_{k+1}$ for $1\leq k\leq n-1$. That is, $v_{k}=v_{0}w^{-k}$ for $v_{0}=v_{1}w^{1}$. Using the fact that $((I-wA)v)_{n}=1$, we obtain that $v_{0}=p_{a}(w^{-1})^{-1}$ where $p_{a}(\cdot)$ is the polynomial $p_{a}(x)=x^{n}+a_{1}x^{n-1}+\dots+a_{n}$. Then we have that $u(I-wA)^{-1}B=\frac{u_{1}w^{-1}+\dots+u_{n}w^{-n}}{p_{a}(w^{-1})}$

Suppose $x_{1},\dots,x_{n}$ are independent variables with mean 0 and covariance matrices and $I$, $U_{1},\dots,U_{d}$ are fixed matrices, then

C Missing proofs in Sections 4 and 5

For any $0<\alpha<\beta$, we have that $\mathcal{B}_{\alpha}\subset\mathcal{B}_{\beta}$.

Proof Let $q_{a}(z)=1+a_{1}z+\dots+a_{n}z^{n}$. Note that $q(z^{-1})=p_{a}(z)/z^{n}$. Therefore we note that $B_{\alpha}=\{a:q_{a}(z)\in\mathcal{C},\forall|z|=1/\alpha\}$. Suppose $a\in\mathcal{B}_{\alpha}$, then $\Re(q_{a}(z))\geq\tau_{1}$ for any $z$ with $|z|=1/\alpha$. Since $\Re(q_{a}(z))$ is the real part of the holomorphic function $q_{a}(z)$, its a harmonic function. By maximum (minimum) principle of the harmonic functions, we have that for any $|z|\leq 1/\alpha$, $\Re(q_{a}(z))\geq\inf_{|z|=1/\alpha}\Re(q_{a}(z))\geq\tau_{1}$. In particular, it holds that for $|z|=1/\beta<1/\alpha$, $\Re(q_{a}(z))\geq\tau_{1}$. Similarly we can prove that for $z$ with $|z|=1/\beta$, $\Re(q_{a}(z))\geq(1+\tau_{0})\Im(q_{a}(z))$, and other conditions for $a$ being in $\mathcal{B}_{\beta}$.