Lyapunov Exponent of Rank One Matrices: Ergodic Formula and Inapproximability of the Optimal Distribution

TL;DR

This paper proves that finding the optimal probability distribution to minimize the Lyapunov exponent for rank-one matrices is NP-hard, contrasting with the polynomial-time solvability of the deterministic case, and explores related ergodic formulas.

Contribution

It establishes NP-hardness of the Lyapunov exponent optimization problem for rank-one matrices and provides an ergodic formula, highlighting fundamental computational complexity differences.

Findings

NP-hardness of stabilizing systems with rank-one matrices

Lyapunov exponent admits a simple quadratic form in p

Contrast with polynomial-time computability of joint spectral radius

Abstract

The Lyapunov exponent corresponding to a set of square matrices $\mathcal{A} = \{A_1, \dots, A_n \}$ and a probability distribution $p$ over $\{1, \dots, n\}$ is $\lambda(\mathcal{A},p) := \lim_{k \to \infty} \frac{1}{k} \,\mathbb{E} \log \|A_{\sigma_k} \cdots A_{\sigma_2}A_{\sigma_1}\|$, where $\sigma_i$ are i.i.d. according to $p$. This quantity is of fundamental importance to control theory since it determines the asymptotic convergence rate $e^{\lambda(\mathcal{A},p)}$ of the stochastic linear dynamical system $x_{k+1} = A_{\sigma_k} x_k$. This paper investigates the following "design problem": given $\mathcal{A}$, compute the distribution $p$ minimizing $\lambda(\mathcal{A},p)$. Our main result is that it is NP-hard to decide whether there exists a distribution $p$ for which $\lambda(\mathcal{A},p)< 0$, i.e. it is NP-hard to decide whether this dynamical system can be stabilized. This hardness result holds even in the "simple"' case where $\mathcal{A}$ contains only rank-one matrices. Somewhat surprisingly, this is in stark contrast to the Joint Spectral Radius -- the deterministic kindred of the Lyapunov exponent -- for which the analogous optimization problem for rank-one matrices is known to be exactly computable in polynomial time. To prove this hardness result, we first observe via Birkhoff's Ergodic Theorem that the Lyapunov exponent of rank-one matrices admits a simple formula and in fact is a quadratic form in $p$. Hardness of the design problem is shown through a reduction from the Independent Set problem. Along the way, simple examples are given illustrating that $p \mapsto \lambda(\mathcal{A},p)$ is neither convex nor concave in general. We conclude with extensions to continuous distributions, exchangeable processes, Markov processes, and stationary ergodic processes.

Tables

Table 1. Table 1 : Complexity overview of the analysis and design questions for the Lyapunov exponent – for collections of general matrices, as well as the special case of rank-one matrices. The various versions of the design problem are for optimizing over different classes of stochastic processes. Our new results are shown in italics.

	Analysis question	Design question (i.i.d.)	Design question (exchangeable)	Design question (Markov)
Rank one	Poly-time algorithm (Theorem 2.1)	NP-hard (Theorem 4.1)	NP-hard (Theorem 5.2)	Poly-time algorithm (Theorem 5.3)
General	Undecidable  [42, Theorem 2]	NP-hard (Theorem 4.1)	NP-hard (Theorem 5.2)	?

Table 2. Table 2 : The Lyapunov exponents λ ​ ( 𝒜 , μ sym ) = λ ​ ( 𝒜 , μ asym ) = − 1 2 ⁡ ( Ψ ​ ( d 2 ) + γ + log ⁡ 4 ) 𝜆 𝒜 subscript 𝜇 sym 𝜆 𝒜 subscript 𝜇 asym 1 2 Ψ 𝑑 2 𝛾 4 \lambda(\mathcal{A},\mu_{\text{sym}})=\lambda(\mathcal{A},\mu_{\text{asym}})=-\operatorname{\frac{1}{2}}(\Psi(\tfrac{d}{2})+\gamma+\log 4) , as derived in Example  4 , for dimensions d ∈ { 2 , … , 8 } 𝑑 2 … 8 d\in\{2,\dots,8\} . Numerical values reported to six significant digits.

d	$-\frac{1}{2}\left(\Psi \left(\frac{d}{2}\right)+\gamma +\log 4\right)$
2	$-0.693147$
3	$-1$
4	$-1.19315$
5	$-1.33333$
6	$-1.44315$
7	$-1.53333$
8	$-1.60981$

Equations

$$\lambda (\mathcal{A},p):=\underset{k\to \infty }{\lim }\frac{1}{k}\mathbb{E}\log \Vert A_{{\sigma }_k}\cdots A_{{\sigma }_2}{A}_{{\sigma }_1}\Vert ,$$

$$x_{k+1}=A_{{\sigma }_k}{x}_k.$$

$$R(\mathcal{A},p):=e^{\lambda (\mathcal{A},p)}.$$

$$\rho (\mathcal{A}):=\underset{k\to \infty }{\lim }\underset{\sigma \in \{1,\dots ,n{\}}^k}{\max }\Vert A_{{\sigma }_k}\cdots A_{{\sigma }_2}{A}_{{\sigma }_1}{\Vert }^{1/k}.$$

$$\lambda (\mathcal{A},p)=\sum _{i,j=1}^n{p}_i{p}_j\log \left|u_i^T{v}_j\right|.$$

$$\frac{1}{k}\log \Vert A_{{\sigma }_k}\cdots A_{{\sigma }_2}{A}_{{\sigma }_1}{\Vert }_2$$

$$=\frac{1}{k}\sum _{t=1}^{k-1}\log \left|u_{{\sigma }_t}^T{v}_{{\sigma }_{t+1}}\right|+\frac{1}{k}\left(\log \Vert u_{{\sigma }_k}{\Vert }_2+\log \Vert v_{{\sigma }_1}{\Vert }_2\right).$$

$$\underset{k\to \infty }{\lim }\frac{1}{k}\log \Vert A_{{\sigma }_k}\cdots A_{{\sigma }_2}{A}_{{\sigma }_1}\Vert \stackrel{\text{a.s.}}{=}\sum _{ij}{p}_i{p}_j\log |u_i^T{v}_j|.$$

$$\lambda (\mathcal{A},p)=p^T{M}_{\mathcal{A}}p,$$

$$\sqrt{-(M_{\mathcal{A}}{)}_{ij}}=\sqrt{-\log |u_i^T{u}_j|}$$

$$\underset{p\in {\Delta }_n}{\min }\lambda (\mathcal{A},p)$$

$$\underset{p\in {\Delta }_n}{\min }R(\mathcal{A},p)$$

$$\underset{p\in {\Delta }_n}{\min }{p}^T{M}_{\mathcal{A}}p,$$

$$I+A_G⩽M⩽O(\log n)\left(I+A_G\right),$$

$$\underset{p\in {\Delta }_n}{\min }{p}^T(I+A_G)p⩽\underset{p\in {\Delta }_n}{\min }{p}^{T}Mp⩽O(\log n)\underset{p\in {\Delta }_n}{\min }{p}^T(I+A_G)p.$$

$$\underset{p\in {\Delta }_n}{\min }\lambda (\mathcal{A},p)\in \left[n^{-1},O(\log n)\right].$$

$$\lambda (\mathcal{A},\hat{p})-\underset{p\in {\Delta }_n}{\min }\lambda (\mathcal{A},p)⩽\delta \left(\underset{p\in {\Delta }_n}{\max }\lambda (\mathcal{A},p)-\underset{p\in {\Delta }_n}{\min }\lambda (\mathcal{A},p)\right).$$

$$\lambda (\mathcal{A},\hat{p})⩽\underset{p\in \Delta }{\min }\lambda (\mathcal{A},p)+\delta \log \frac{1}{\gamma }$$

$$\lambda(\mathcal{A},p)=p^{T}M_{\mathcal{A}}p=\sum_{ij}p_{i}p_{j}(M_{\mathcal{A}})_{ij}\leqslant\Big{(}\sum_{ij}p_{i}p_{j}\Big{)}\max_{ij}(M_{\mathcal{A}})_{ij}=\max_{ij}(M_{\mathcal{A}})_{ij}\leqslant\max_{i}(M_{\mathcal{A}})_{ii},$$

$$\lambda(\mathcal{A},p)=p^{T}M_{\mathcal{A}}p=\sum_{ij}p_{i}p_{j}(M_{\mathcal{A}})_{ij}\geqslant\Big{(}\sum_{ij}p_{i}p_{j}\Big{)}\min_{ij}(M_{\mathcal{A}})_{ij}=\min_{ij}(M_{\mathcal{A}})_{ij}.$$

$$\underset{k\to \infty }{\lim }\frac{1}{k}\log \Vert A_{{\sigma }_k}\cdots A_{{\sigma }_1}\Vert$$

$$\frac{1}{k}\log \Vert A_{{\sigma }_k}\cdots A_{{\sigma }_1}\Vert \to \int \log \left|u^T{v}^{\prime }\right|dP(u{v}^T,u^{\prime }{v}^{\prime T})$$

$$\frac{1}{k}\log \Vert A_{{\sigma }_k}\cdots A_{{\sigma }_2}{A}_{{\sigma }_1}{\Vert }_2=\frac{1}{k}\sum _{t=1}^{k-1}\log \left|u_{{\sigma }_t}^T{v}_{{\sigma }_{t+1}}\right|+\frac{1}{k}\log \Vert u_{{\sigma }_k}{\Vert }_2.$$

$$\frac{1}{k}\left(\log {\Vert u_{{\sigma }_k}\Vert }_2+\log {\Vert v_{{\sigma }_1}\Vert }_2\right)\to 0,$$

$$\frac{1}{k}\sum _{t=1}^k\log \left|u_{{\sigma }_t}^T{v}_{{\sigma }_{t+1}}\right|\to \int \log \left|u^T{v}^{\prime }\right|dP(u{v}^T,u^{\prime }{v}^{\prime T}),$$

$$\frac{1}{k}\log \Vert A_{{\sigma }_k}\cdots A_{{\sigma }_1}\Vert \to \int \int \log |u^T{v}^{\prime }|d\mu (u{v}^T)d\mu (u^{\prime }{v}^{\prime T})$$

$$\lambda (\mathcal{A},{\mu }_{\text{sym}})=\lambda (\mathcal{A},{\mu }_{\text{asym}})={\mathbb{E}}_{u,v\sim \mathsf{Unif}(S^{d-1})}\log \left|u^{T}v\right|.$$

$$\frac{d\nu (t)}{dt}=\frac{\Gamma \left(\frac{d}{2}\right)}{\sqrt{\pi }\Gamma \left(\frac{d-1}{2}\right)}{\left(1-t^2\right)}^{(d-3)/2}\text{mathds}{1}_{t\in (-1,1)},$$

$${\mathbb{E}}_{u,v\sim \mathsf{Unif}(S^{d-1})}\log \left|u^{T}v\right|$$

$$=-\frac{\Psi \left(\frac{d}{2}\right)+\gamma +\log 4}{2},$$

$$\lambda (\mathcal{A},{\mu }_{\text{sym}})=\lambda (\mathcal{A},{\mu }_{\text{asym}})=-\frac{1}{2}(\log d+\gamma +\log 2)+O(d^{-1}).$$

Full text

This hardness result holds even in the “simple” case where $\mathcal{A}$ contains only rank-one matrices. Somewhat surprisingly, this is in stark contrast to the Joint Spectral Radius – the deterministic kindred of the Lyapunov exponent – for which the analogous optimization problem for rank-one matrices is known to be exactly computable in polynomial time.

To prove this hardness result, we first observe via Birkhoff’s Ergodic Theorem that the Lyapunov exponent of rank-one matrices admits a simple formula and in fact is a quadratic form in $p$. Hardness of the design problem is shown through a reduction from the Independent Set problem. Along the way, simple examples are given illustrating that $p\mapsto\lambda(\mathcal{A},p)$ is neither convex nor concave in general. We conclude with extensions to continuous distributions, exchangeable processes, Markov processes, and stationary ergodic processes.

00footnotetext: The authors are with the Laboratory for Information and Decision Systems (LIDS), Massachusetts Institute of Technology, Cambridge MA 02139. Work partially supported by NSF AF 1565235 and NSF Graduate Research Fellowship 1122374.

1 Introduction

Given a finite collection of square matrices $\mathcal{A}=\{A_{1},\dots,A_{n}\}$ and a probability distribution $p$ over $\{1,\dots,n\}$, the corresponding Lyapunov exponent $\lambda(\mathcal{A},p)$ is defined as:

[TABLE]

where the $\sigma_{i}$ are independently and identically distributed (i.i.d.) according to $p$. Over the past half century, this quantity $\lambda(\mathcal{A},p)$ has received significant interest due to its many connections and applications to diverse fields including probability, ergodic theory, dynamical systems, functional analysis, representation theory, computer image generation of fractals, control theory, and more; see e.g.  the surveys [10, 9, 16, 21, 43, 4] and references within.

The primary application of the Lyapunov exponent that this paper focuses on is in control theory. The fundamental connection is that $\lambda(\mathcal{A},p)$ dictates the asymptotic growth rate of the stochastic linear dynamical system

[TABLE]

Indeed, by the celebrated Furstenberg-Kesten Theorem [23], the average norm $\lim_{k\to\infty}\|A_{\sigma_{k}}\cdots A_{\sigma_{2}}A_{\sigma_{1}}\|^{1/k}$ almost surely (a.s.) exists and is equal to the Lyapunov spectral radius

[TABLE]

In particular, the dynamical system (1.2) converges asymptotically111 The Lyapunov exponent characterizes a.s. convergence of (1.2). In general, this is different from convergence in mean, i.e., $\lim_{k\to\infty}\mathbb{E}\|A_{\sigma_{k}}\cdots A_{\sigma_{1}}\|=0$. For example, consider the scalar case with $n=2$, where $A_{\sigma_{k}}$ are i.i.d., and are either $2$ or $\tfrac{1}{100}$ with probability $\operatorname{\frac{1}{2}}$. Then $A_{\sigma_{k}}\cdots A_{\sigma_{1}}$ does not converge in mean (in fact $\mathbb{E}|A_{\sigma_{k}}\cdots A_{\sigma_{1}}|$ grows exponentially in $k$) due to a few exponentially large but exponentially unlikely trajectories, but does converge to [math] with probability $1$ (since $\exp(\mathbb{E}\log|A_{\sigma_{1}}|)<1$). to zero (a.s.) if and only if $R(\mathcal{A},p)<1$, i.e. if and only if $\lambda(\mathcal{A},p)<0$. Note that in this sense, the Lyapunov spectral radius is a natural stochastic generalization of the spectral radius of a single matrix. Two natural control-theoretic questions arise about the dynamical system (1.2): analysis and design.

The analysis problem is to compute the convergence rate $\lambda(\mathcal{A},p)$ of the dynamical system given (fixed) parameters $(\mathcal{A},p)$. This problem is well-known to be difficult, and approximating $\lambda(\mathcal{A},p)$ is in general algorithmically undecidable [42, Theorem 2]. Nevertheless, the significance of this problem has motivated a substantial body of research on algorithms for obtaining upper or lower bounds on $\lambda(\mathcal{A},p)$, and on identifying special cases in which this problem becomes tractable; see e.g. [13, 14, 32, 34, 17, 37, 24, 18, 30, 6, 38].

The design problem – which is the focus of this paper – is to optimize the dynamical system so that the convergence is as fast as possible, namely: given a fixed set of matrices $\mathcal{A}$, compute a distribution $p$ minimizing $\lambda(\mathcal{A},p)$. In contrast to the analysis question, there seems to be little previous work on this design problem. (A related version of this design problem uses instead the more restrictive notion of convergence in mean, see e.g. [8].) One of this paper’s primary contributions is to explain this lack of progress on the design problem by showing that it is NP-hard to determine whether there exists a distribution $p$ such that $\lambda(\mathcal{A},p)<0$. That is, unless $\textsc{P}=\textsc{NP}$, there is no polynomial time algorithm for deciding whether the dynamical system (1.2) has a stabilizing distribution $p$. Moreover, this complexity result holds even in the “simple” case where all matrices in $\mathcal{A}$ are rank one.

1.1 Contributions and outline

This paper focuses on the case when all matrices in $\mathcal{A}$ are rank one. The advantage of restricting to this simple case is that (i) the analysis problem becomes tractable as $\lambda(\mathcal{A},p)$ admits a simple, closed-form expression; and (ii) this formula enables us to prove that the design problem is hard in a complexity-theoretic sense. We elaborate below on details and implications. Table 1 contains an overview of our main results.

Section 2 studies the analysis problem when $\mathcal{A}$ is a collection of rank-one matrices. We observe in Theorem 2.1 that in this special case, $\lambda(\mathcal{A},p)$ admits a simple closed-form expression and is in fact a quadratic form in $p$. This is an ergodic-type result since it equates $\lambda(\mathcal{A},p)$ – which is a “time average” of the product of an infinite sequence of random matrices – to a simpler “space average.” We note that this formula readily extends to ergodic stationary processes taking values in sets $\mathcal{A}$ consisting of infinitely many rank-one matrices; however, for expository reasons we defer this generalization to Section 5, and provide in Section 2 just the simple case of i.i.d. processes taking values in finite sets $\mathcal{A}$, since this admits an elementary proof and requires essentially no measure theory. We remark that this formula shows that for this rank-one case, $\lambda(\mathcal{A},p)$ is computable in polynomial time (as opposed to being undecidably hard to approximate in the general case [42]).

Sections 3 and 4 then turn to the design problem. Specifically, Section 3 investigates the convexity and concavity of the function $p\mapsto\lambda(\mathcal{A},p)$ over the probability simplex $\Delta_{n}$. By our aforementioned result that $\lambda(\mathcal{A},p)$ is a quadratic form in this rank-one case, these convexity/concavity properties are equivalent to the (conditional) definiteness of the Hessian $\tfrac{\partial^{2}}{\partial p^{2}}\lambda(\mathcal{A},p)$. We give simple, small examples illustrating that in general this matrix is neither (conditionally) positive semidefinite nor negative semidefinite, and then explain this phenomenon by making a connection to the fact that the Martin distance on the $(1,n)$ Grassmanian is not a metric.

Of course, this lack of convexity of the function $p\mapsto\lambda(\mathcal{A},p)$ does not in itself prove that the design problem $\min_{p\in\Delta_{n}}\lambda(\mathcal{A},p)$ is hard. Section 4 formally shows that it is NP-hard to even determine the sign of $\min_{p\in\Delta_{n}}\lambda(\mathcal{A},p)$. That is, unless $\textsc{P}=\textsc{NP}$, there is no polynomial-time algorithm for determining whether the dynamical system (1.2) has a stabilizing distribution $p$. This explains the lack of prior algorithmic progress on the design problem. The key idea in the proof is a reduction from the Independent Set problem, which is known to be NP-hard to approximate and is also known to be expressible as the minimization of a certain quadratic form over the probability simplex. The intuition for this reduction is, roughly speaking, that the optimal distribution for $\lambda(\mathcal{A},p)$ should concentrate mass on a set of matrices in $\mathcal{A}$ which are as mutually “orthogonal” to each other as possible, and such a set of matrices can be thought of as an independent set in a certain associated graph. (See the proof of Theorem 4.1 for precise details.)

These hardness results uncover an interesting and somewhat surprising (at least to the authors) dichotomy between the complexity of optimizing the Lyapunov exponent and its deterministic analogue, the Joint Spectral Radius (JSR). Recall that the JSR of a collection of square matrices $\mathcal{A}=\{A_{1},\dots,A_{n}\}$ is defined as:

[TABLE]

The JSR is a natural generalization of the spectral radius of a single matrix and characterizes the maximal (i.e., worst-case) asymptotic growth rate of the linear dynamical system (1.2) over all possible choices of $A_{\sigma_{k}}$. Due to its many applications (see e.g. [29] and the references within), the JSR has also attracted a large body of research. Similarly to the Lyapunov exponent, the JSR is known to in general have algorithmic undecidability and NP-hardness complexity barriers [42, 7], yet nevertheless there has been significant work on approximation algorithms for $\rho(\mathcal{A})$ and on identifying special cases in which the problem becomes more tractable [36, 39, 26, 1, 29]. In particular, it is known that (optimizing) the JSR of rank-one matrices is equivalent to a certain discrete optimization problem over an associated graph, namely the Maximum-Cycle-Mean problem, which can be solved in polynomial time [27, 2, 33, 31]. This is in stark contrast to our result that the Lyapunov exponent of rank-one matrices is NP-hard to optimize, even approximately.

In Subsection 4.2, we note that in the special case when the set $\mathcal{A}$ is “well-conditioned,” there is a polynomial-time approximation algorithm for $\min_{p\in\Delta_{n}}\lambda(\mathcal{A},p)$.

Finally, Section 5 considers extensions to more general stochastic processes $A_{\sigma_{1}}$, $A_{\sigma_{2}}$, $\dots$. The corresponding analysis and design problems are discussed in Subsection 5.1 and 5.2, respectively. Specifically, Subsection 5.1 shows that a similar formula holds for the “Lyapunov exponent” of ergodic stationary processes of rank-one matrices. This allows for the study of more exotic stochastic linear dynamical systems (1.2). Using this formula, we compute the convergence rate for setups including i.i.d. processes over continuous distributions, exchangeable processes, and Markov processes. As an example application, we explicitly calculate the Lyapunov exponent for the induced distribution over matrices $uu^{T}$, where $u$ is distributed according to the $d$-dimensional spherical measure. Subsection 5.2 then turns to the design problem for these more general setups. We show that optimizing the convergence rate of the dynamical system (1.2) is NP-hard over exchangeable processes, but can be done in polynomial time for Markov processes. The intuition for these results is as follows. First, we observe that (after some manipulation) these problems seek to optimize a linear functional over a convex set, and thus are equivalent to optimizing over the extreme points. Informally speaking, the extreme points of the sets of exchangeable processes and Markov processes, respectively, are the sets of i.i.d. process and deterministic processes. The arguments are then finished by recalling that optimizing the convergence rate over i.i.d. processes is NP-hard (discussed above), while optimizing over deterministic processes is known to be polynomial-time computable.

1.2 Notation

The set of discrete probability distributions on $n$ atoms is identified with the simplex $\Delta_{n}:=\{p\in\mathbb{R}_{\geqslant 0}^{n}:\sum_{i=1}^{n}p_{i}=1\}$. The support of a discrete probability distribution $p\in\Delta_{n}$ is denoted by $\operatorname{supp}(p):=\{i\in[n]:p_{i}>0\}$. We write $X\sim p$ to denote that $X$ is a random variable with distribution $p$. The symbols $e_{i}$, $\mathbf{1}$, and $I$ denote the $i$th standard basis vector, all-ones vector, and identity matrix, respectively, in an ambient dimension that will be clear from context. The norm $\|\cdot\|_{2}$ denotes the Euclidean norm when applied to vectors, and denotes the spectral norm when applied to matrices. All logarithms are natural logarithms, i.e., taken base $e$. We abbreviate “almost surely” by a.s., “without loss of generality” by w.l.o.g., “independent and identically distributed” by i.i.d., “positive semidefinite” by PSD, and “negative semidefinite” by NSD. Recall that a square matrix $M$ is said to be conditionally PSD (resp. conditionally NSD) if $p^{T}Mp\geqslant 0$ (resp. $\leqslant 0$) for all vectors $p$ in the orthogonal component of the subspace spanned by $\mathbf{1}$. All other specific notations are introduced in the main text.

2 Ergodic formula for Lyapunov exponent of rank one matrices

This section gives an explicit closed-form expression for the Lyapunov exponent of a collection of rank-one matrices. This is an ergodic-type result in the sense that it equates the Lyapunov exponent – which is a “time average” – to a simpler “space average.”

This result readily extends to stationary ergodic processes where the collection $\mathcal{A}$ contains infinitely many matrices (Theorem 5.1). However for simplicity of exposition, we first present the discrete case (Theorem 2.1) since it admits an elementary proof with essentially no measure theory.

Theorem 2.1**.**

For any collection of rank-one square matrices $\mathcal{A}:=\{A_{i}=u_{i}v_{i}^{T}\}_{i=1}^{n}$ and any probability distribution $p\in\Delta_{n}$,

[TABLE]

Remark 2.1**.**

The formula in Theorem 2.1 is independent of the rank-one decompositions $A_{i}=u_{i}v_{i}^{T}$. Indeed, a rescaling of $u_{i}$ and $v_{i}$ changes the cost matrix $[\log|u_{i}^{T}v_{j}|]_{ij}$ by a skew-symmetric matrix, which does not affect the underlying quadratic form $\sum_{ij}p_{i}p_{j}\log\left\lvert u_{i}^{T}v_{j}\right\rvert$.

It is important to emphasize that, in general, the Lyapunov exponent of a set of matrices cannot be computed in closed form. It is well-known that (under mild conditions) the Lyapunov exponent is equal to a simpler “space average” over a certain invariant distribution [22]; however, the complexity of such ergodic formulas lies in the invariant distribution, which is often difficult – or impossible – to characterize or manipulate analytically. The key difference for the rank-one setting studied here is that in this special case, the invariant distribution has a simple closed form. In words, this is because in this rank-one case, the contraction of the next matrix $A_{\sigma_{k+1}}$ depends on the previous matrices essentially only through the last matrix $A_{\sigma_{k}}$.

Proof of Theorem 2.1.

Note that the definition of $\lambda(\mathcal{A},p)$ is independent of the norm $\|\cdot\|$; for convenience of computation, let us take it to be the spectral norm. Let $\sigma_{1},\sigma_{2},\dots$ be i.i.d. according to $p$, and expand the (random) quantity

[TABLE]

We may assume w.l.o.g. that $A_{i}\neq 0$ for all $i\in\operatorname{supp}(p)$, since otherwise the proof is immediate as both sides of (2.1) equal $-\infty$. Thus the latter term $\frac{1}{k}(\log\|u_{\sigma_{k}}\|_{2}+\log\|v_{\sigma_{1}}\|_{2})$ in (2.2) converges to [math] as $k\to\infty$.

Now let $G_{\mathcal{A}}$ be the complete directed graph on $n$ nodes, one for each matrix $A_{i}$, with weight $\log|u_{i}^{T}v_{j}|$ on the edge from node $i$ to node $j$. Consider the $n$-state Markov chain on $G_{\mathcal{A}}$ with transition probability $p_{j}$ from node $i$ to node $j$, irrespective of $i$. The key point is that after multiplying by $\tfrac{k}{k-1}$, the sum $\frac{1}{k-1}\sum_{t=1}^{k-1}\log|u_{\sigma_{t}}^{T}v_{\sigma_{t+1}}|$ in (2.2) is precisely the average weight of edges travelled when running this Markov chain for $k-1$ steps from initial state $\sigma_{1}\sim p$. Now since this Markov chain has stationary distribution $p$, a standard ergodic result on discrete-time, finite-state Markov chains (see e.g. [25]) shows that this quantity converges a.s. to $\sum_{ij}p_{i}p_{j}\log|u_{i}^{T}v_{j}|$. Thus by (2.2) and the fact $\lim_{k\to\infty}\tfrac{k-1}{k}=1$, we conclude that

[TABLE]

By sub-multiplicativity of the spectral norm, $|\frac{1}{k}\log\lVert{A_{\sigma_{k}}\cdots A_{\sigma_{2}}A_{\sigma_{1}}}\rVert_{2}|\leqslant\max_{i\in\operatorname{supp}(p)}|\log\|A_{i}\|_{2}|$ is uniformly bounded for all $k\in\mathbb{N}$. Thus the claimed formula for $\lambda(\mathcal{A},p)$ follows from the a.s. convergence proven in (2.3) and an application of Lebesgue’s Dominated Convergence Theorem. ∎

Remark 2.2**.**

Theorem 2.1 shows that the Lyapunov exponent of rank-one matrices depends only on the products $A_{i}A_{j}$ of length two. This is in contrast to the Joint Spectral Radius of rank-one matrices, which depends on all matrix products of length up to $n$ [2, Corollary 2.1].

Remark 2.3**.**

Theorem 2.1 shows that the Lyapunov exponent of rank-one matrices is not only computable, but also computable in polynomial time. This is in sharp contrast to the Lyapunov exponent of general matrices (i.e., not necessarily rank-one), which is undecidably hard to approximate [42, Theorem 2].

Remark 2.4**.**

It is an interesting open question whether such results extend to collections of matrices with some small fixed rank $k>1$. In particular, is there a closed-form expression for the Lyapunov exponent, analogously to Theorem 2.1? And if so, how do the phenomena mentioned in Remark 2.2 and 2.3 change?

3 Convexity properties of the Lyapunov exponent as a function of the distribution

From Theorem 2.1, we know that $\lambda(\mathcal{A},p)$ is a quadratic form in $p$:

[TABLE]

where $(M_{\mathcal{A}})_{ij}:=\operatorname{\frac{1}{2}}(\log|u_{i}^{T}v_{j}|+\log|u_{j}^{T}v_{i}|)$. This resolves the analysis problem of computing the Lyapunov exponent given a collection of rank-one matrices $\mathcal{A}$ and a probability distribution $p$. In this and the following section, we now turn to the design problem: Given a collection of matrices $\mathcal{A}$, find the distribution $p\in\Delta_{n}$ minimizing $\lambda(\mathcal{A},p)$. As mentioned in the introduction, this is a fundamental problem in control theory with applications including e.g., fast convergence in the corresponding stochastic linear dynamical system (1.2).

We begin with an example illustrating that even in simple cases, the optimal distribution can depend in a non-obvious way on the geometric configuration of the matrices in $\mathcal{A}$.

Example 1**.**

We consider several cases of symmetric rank-one matrices $\mathcal{A}=\{u_{i}u_{i}^{T}\}_{i=1}^{3}$ (plotted in Figure 1) to illustrate the design problem $\min_{p\in\Delta_{n}}\lambda(\mathcal{A},p)$. This problem has a simple and intuitive geometric interpretation in terms of the linear dynamical system222This can be viewed e.g., as a randomized instance of Dykstra’s alternating-projection method for finding a point at the intersection of convex sets (namely the lines $u_{i}$ here) [20]. $x_{k+1}=A_{\sigma_{k}}x_{k}=u_{\sigma_{k}}u_{\sigma_{k}}^{T}x_{k}$, which corresponds to projecting (and scaling) the current state vector $x_{k}$ onto a randomly chosen vector $u_{\sigma_{k}}$. Specifically, $\min_{p\in\Delta_{n}}\lambda(\mathcal{A},p)$ corresponds to finding the distribution $p\in\Delta_{3}$ over lines $\{u_{i}\}_{i=1}^{3}$ such that repeatedly projecting (and scaling) onto a random vector $u_{i}$ with probability $p_{i}$ converges as quickly as possible. In all examples, $u_{1}=[1,10]^{T}$ and $u_{2}=[10,1]^{T}$ are close to orthogonal in that they have small inner product; this is why the optimal distribution is concentrated on them in Figure 1(a). In Figure 1(b), $u_{3}$ has smaller norm, and the optimal distribution is concentrated just on it. In Figure 1(c), the optimal distribution is supported on all three matrices. Finally, we note that the optimal distribution for the maximization problem $\max_{p\in\Delta_{n}}\lambda(\mathcal{A},p)$ is easily characterized: it is always concentrated on the vector(s) of largest norm (see Lemma 4.2 below).

In the rest of the section, we investigate the (lack of) convexity and concavity of the function $p\mapsto\lambda(\mathcal{A},p)$ for a fixed collection of rank-one matrices $\mathcal{A}=\{u_{i}v_{i}^{T}\}_{i=1}^{n}$. Algorithmic and complexity results for the design problem $\min_{p\in\Delta_{n}}\lambda(\mathcal{A},p)$ are then discussed in the following section.

Note that since $p\in\Delta_{n}$ is a probability distribution in our context, the function $\Delta_{n}\ni p\mapsto\lambda(\mathcal{A},p)=p^{T}M_{\mathcal{A}}p$ is convex (resp. concave) if and only if the matrix $M_{\mathcal{A}}$ is conditionally PSD (resp. conditionally NSD). We begin by observing that these conditional semidefiniteness properties of $M_{\mathcal{A}}$ are invariant under rescaling of the matrices in $\mathcal{A}$.

Lemma 3.1**.**

The conditional semidefiniteness of $M_{\mathcal{A}}$ is invariant under non-zero rescalings of matrices in $\mathcal{A}$.

Proof.

It suffices to show this invariance holds under rescaling $u_{i}$ by some $\alpha\neq 0$. This scaling changes $M_{\mathcal{A}}$ to $M_{\mathcal{A}}+\Delta$ where $\Delta:=\operatorname{\frac{1}{2}}\log|\alpha|(e_{i}\mathbf{1}^{T}+\mathbf{1}e_{i}^{T})$. Now for any vector $x$ orthogonal to $\mathbf{1}$, $x^{T}\Delta x=0$, thus in particular $x^{T}(M_{\mathcal{A}}+\Delta)x=x^{T}M_{\mathcal{A}}x$. Therefore $M_{\mathcal{A}}+\Delta$ is conditionally PSD (resp. NSD) if and only if $M_{\mathcal{A}}$ is. ∎

As seen next, $M_{\mathcal{A}}$ is neither conditionally PSD nor conditionally NSD in general. This is true even if all matrices in $\mathcal{A}$ are symmetric, i.e. when each $u_{i}=v_{i}$, so the rest of this section henceforth restricts to this symmetric case $\mathcal{A}=\{u_{i}u_{i}^{T}\}_{i=1}^{n}$ for simplicity of exposition.

Example 2**.**

Let $\mathcal{A}=\{e_{1}e_{1}^{T},e_{2}e_{2}^{T}\}\subset\mathbb{R}^{2\times 2}$, where $e_{1}:=[1,0]^{T}$ and $e_{2}:=[0,1]^{T}$. Then $M_{\mathcal{A}}=\begin{bmatrix}0&-\infty\\ -\infty&0\end{bmatrix}$ is not conditionally NSD since $\begin{bmatrix}1&-1\end{bmatrix}M_{\mathcal{A}}\begin{bmatrix}1\\ -1\end{bmatrix}=\infty\not\leqslant 0$.

Example 3**.**

Let $\mathcal{A}=\{e_{1}e_{1}^{T},e_{2}e_{2}^{T},(e_{1}+e_{2})(e_{1}+e_{2})^{T}\}\subset\mathbb{R}^{2\times 2}$, where $e_{1}:=[1,0]^{T}$ and $e_{2}:=[0,1]^{T}$. Then $M_{\mathcal{A}}=\begin{bmatrix}0&-\infty&0\\ -\infty&0&0\\ 0&0&\log 2\end{bmatrix}$ is not conditionally PSD since $\begin{bmatrix}1&1&-2\end{bmatrix}M_{\mathcal{A}}\begin{bmatrix}1\\ 1\\ -2\end{bmatrix}=-\infty\not\geqslant 0$.

Note that Examples 2 and 3 are easily modified so that $M_{\mathcal{A}}$ only has finite entries: add $\varepsilon\mathbf{1}$ to each of the vectors $u_{i}$, for some small $\varepsilon>0$.

We conclude the section by connecting this lack of convexity of $p\mapsto\lambda(\mathcal{A},p)$ to the Martin distance on the space of lines (i.e., the $(1,n)$ Grassmanian). For simplicity, we continue to restrict to the symmetric case $A_{i}=u_{i}u_{i}^{T}$ and assume that the vectors $u_{i}$ are normalized to have unit Euclidean norm (this does not change convexity by Lemma 3.1). Since $-M_{\mathcal{A}}$ has diagonal entries of zero (because of the normalization $\|u_{i}\|_{2}=1$) and off-diagonal entries that are strictly positive (by Cauchy-Schwarz), a celebrated result of Schoenberg ensures that $-M_{\mathcal{A}}$ is conditionally NSD if and only if there is an embedding $z_{1},\dots,z_{n}\in\ell_{2}$ such that the Euclidean distance $\|z_{i}-z_{j}\|_{2}=\sqrt{-(M_{\mathcal{A}})_{ij}}$ for each $i,j\in\{1,\dots,n\}$ [40, Theorem 1]. In particular, this implies that if $M_{\mathcal{A}}$ were conditionally PSD, then $\{\sqrt{-(M_{\mathcal{A}})_{ij}}\}_{i,j=1}^{n}$ satisfy the triangle inequality. However,

[TABLE]

is precisely the Martin distance between the lines in directions $u_{i}$ and $u_{j}$, and the Martin distance is well-known to not satisfy the triangle inequality (i.e. it is not a metric333E.g., the construction from Example 3: $u_{1}=e_{1}$, $u_{2}=e_{2}$, and $u_{3}=e_{1}+e_{2}$. Then $\sqrt{-\log|u_{1}^{T}u_{3}|}+\sqrt{-\log|u_{2}^{T}u_{3}|}=0+0<\infty=\sqrt{-\log|u_{1}^{T}u_{2}|}$, violating the triangle inequality.) [15, Section 12]. Therefore $M_{\mathcal{A}}$ cannot always be conditionally PSD.

4 Optimizing the distribution

We now investigate the design problem from an algorithmic viewpoint. First we show in Subsection 4.1 that this problem is NP-hard in general; then in Subsection 4.2, we briefly note that in a certain “well-conditioned” case, there is a polynomial-time approximation algorithm.

4.1 Inapproximability

The main result here (Theorem 4.1) is that even the stabilization decision problem for Lyapunov exponents is NP-hard. That is, unless $\textsc{P}=\textsc{NP}$, there is no polynomial-time algorithm for deciding whether there exists a distribution $p\in\Delta_{n}$ satisfying $\lambda(\mathcal{A},p)<0$. Moreover, these hardness results hold even in the “simple” case where $\mathcal{A}$ contains only rank-one symmetric matrices.

Theorem 4.1**.**

Given a collection $\mathcal{A}$ of $n$ symmetric rank-one matrices, it is NP-hard to decide whether

[TABLE]

is negative.

A few remarks. First, note that Theorem 4.1 can be equivalently restated in terms of the Lyapunov spectral radius $R(\mathcal{A},p)=e^{\lambda(\mathcal{A},p)}$ (see (1.3)). That is, it is NP-hard to decide whether

[TABLE]

is less than $1$ for a collection $\mathcal{A}$ of $n$ symmetric rank-one matrices. Recall from the introduction that $R(\mathcal{A},p)<1$ is equivalent to the (a.s.) asymptotic convergence of the stochastic linear dynamical system (1.2).

Next, we emphasize that this hardness result is somewhat surprising (at least to the authors) in the sense that the analogous optimization problem for the Joint Spectral Radius – which is the deterministic kindred of the Lyapunov exponent (see the introduction for details) – can be computed exactly in polynomial time for rank-one matrices [2, Remark 2.1].

We now turn to proving Theorem 4.1. Note that by Theorem 2.1, it suffices to show the NP-hardness of deciding the sign of the following quadratic optimization problem over the simplex:

[TABLE]

where $(M_{\mathcal{A}})_{ij}=\log\left\lvert u_{i}^{T}u_{j}\right\rvert$. We observed in the previous Section 3 that, in general, this optimization problem (4.3) is not convex since $M_{\mathcal{A}}$ is not conditionally PSD. Although this lack of convexity is a first step towards suggesting the hardness of the minimization problem (4.3), it of course does not constitute a hardness proof in itself.

We will prove Theorem 4.1 via a reduction from approximating the Independent Set problem. It is helpful to first recall two classical results, namely the Motzkin-Straus formulation for the Independent Set problem as a quadratic optimization over the simplex, and the hardness of approximation of the Independent Set problem. Below, $\alpha(G)$ denotes the independence number of a graph $G$, i.e. the size of the largest subset of nodes that has no edges between them. $A_{G}$ denotes the adjacency matrix of $G$, i.e. the entry $(A_{G})_{ij}$ is $1$ if $(i,j)$ is an edge of $G$ and is [math] otherwise.

Theorem 4.2**.**

[35*, Corollary 1]**

For any undirected graph $G$, $\frac{1}{\alpha(G)}=\min_{p\in\Delta_{n}}p^{T}(I+A_{G})p$.*

Theorem 4.3**.**

[44, Theorem 1.1]444The complexity result stated in [44] – as well in much of the related literature – is for the problem of computing the size of the largest clique in a graph. However, this problem is well-known to be equivalent to the independent set problem, since a set of nodes is a clique in a graph if and only if it is an independent set in the complement graph. For all $\varepsilon>0$, it is NP-hard to approximate the independence number of a graph on $n$ nodes to within a multiplicative factor of $n^{1-\varepsilon}$.

The reduction for proving Theorem 4.1 is as follows. Given an undirected graph $G$ on $n$ nodes and with edge set $E$, let $B:=3nI+\exp[I+A_{G}]$, where the exponential is applied entrywise. Note that $B$ is symmetric and diagonally dominant since for each $i\in[n]$, $B_{ii}=3n+e\geqslant n+(e-1)\cdot\text{degree}(i)-1=\sum_{j\neq i}B_{ij}$. Thus $B$ is PSD and in particular admits a factorization $B=U^{T}U$ for some $U\in\mathbb{R}^{d\times n}$, $d\leqslant n$. Denote the columns of $U$ by $u_{1},\dots,u_{n}\in\mathbb{R}^{d}$, and let $M\in\mathbb{R}^{n\times n}$ be the matrix with entries $M_{ij}:=\log|u_{i}^{T}u_{j}|=\log|(U^{T}U)_{ij}|=\log|B_{ij}|=\log B_{ij}$. Now the key observation is that $\min_{p\in\Delta_{n}}p^{T}Mp$ is an $O(\log n)$ multiplicative approximation to $\tfrac{1}{\alpha(G)}$:

Lemma 4.1**.**

$\frac{1}{\alpha(G)}\leqslant\min_{p\in\Delta_{n}}p^{T}Mp\leqslant O(\log n)\frac{1}{\alpha(G)}$.

Proof.

Since $M=(I+A_{G})+(\log(3n+e)-1)I$, we have

[TABLE]

where the matrix inequalities are in the entrywise sense. Thus

[TABLE]

The proof is then complete by applying Theorem 4.2. ∎

The last ingredient we will need for the proof of Theorem 4.1 is the simple but helpful observation that $\lambda(\mathcal{A},p)$ additively shifts if all matrices in $\mathcal{A}$ are scaled by the same constant. This follows immediately from the definition of the Lyapunov exponent – or equivalently, by homogeneity of the Lyapunov spectral radius.

Observation 4.1**.**

Let $\mathcal{A}$ be some collection of square matrices, and let $c\mathcal{A}:=\{cA:A\in\mathcal{A}\}$ for some positive scalar $c$. Then $\lambda(c\mathcal{A},p)=\lambda(\mathcal{A},p)+\log c$ for all $p\in\Delta_{n}$.

Proof of Theorem 4.1.

Assume there exists a polynomial-time algorithm for this decision problem. Now consider an undirected graph $G$ on $n$ nodes, and note that the aforementioned matrix $M$ and vectors $\{u_{i}\}_{i=1}^{n}$ are computable in polynomial time in $n$. Letting $\mathcal{A}:=\{u_{i}u_{i}^{T}\}_{i=1}^{n}$, we have by Theorem 2.1 that $\lambda(\mathcal{A},p)=p^{T}Mp$. Thus by Lemma 4.1 and the trivial bound $\alpha(G)\in[1,n]$, it holds that

[TABLE]

By Observation 4.1, we may – via binary search – refine the confidence interval in (4.6) to size $n^{-1}$ using $O(\log(\tfrac{O(\log n)-n^{-1}}{n^{-1}}))=O(\log n)$ queries to the assumed decision algorithm on $c\mathcal{A}=\{cA:A\in\mathcal{A}\}$, for appropriate scalings $c$. This yields a polynomial-time555To be precise, this requires checking that the bit-complexity of the inputs to the decision algorithm is also polynomial in $n$; however, this is not an issue. Specifically, in light of Observation 4.1 and the original confidence interval (4.6), each scaling $c$ used in the binary search satisfies $\log c\in[n^{-1},O(\log n)]$. Moreover, since we refine the search interval only to size $O(n^{-1})$, it suffices to truncate each $\log c$ to precision $O(n^{-1})$. Together, this implies that the bit complexity of the scaled matrices $c\mathcal{A}$ is amplified from that of $\mathcal{A}$ by at most a polynomial factor in $n$. approximation to $\min_{p\in\Delta_{n}}\lambda(\mathcal{A},p)$ to additive error $n^{-1}$, and therefore by Lemma 4.1 yields a polynomial-time approximation $\alpha(G)$ within a multiplicative factor of $O(\log n)$. But this is NP-hard by Theorem 4. ∎

We remark that although the above reduction was for NP-hardness of deciding the sign of $\min_{p\in\Delta_{n}}\lambda(\mathcal{A},p)$, it also demonstrates an NP-hardness barrier for approximating the value of $\min_{p\in\Delta_{n}}\lambda(\mathcal{A},p)$. This is recorded briefly below. Note that we show hardness for approximating $\min_{p\in\Delta_{n}}\lambda(\mathcal{A},p)$ to additive error (rather than multiplicative) since this implies hardness for multiplicatively approximating the optimal Lyapunov spectral radius $R(\mathcal{A},p)=e^{\lambda(\mathcal{A},p)}$, which is the relevant rate of the dynamical system (1.2).

Theorem 4.4**.**

For all $\varepsilon>0$, it is NP-hard to approximate $\min_{p\in\Delta_{n}}\lambda(\mathcal{A},p)$ to an additive error of $n^{-\varepsilon}$ for a collection $\mathcal{A}$ of $n$ symmetric rank-one matrices.

Proof.

We begin in the same way as in the proof of Theorem 4.1: given an undirected graph $G$ on $n$ nodes, construct the matrix $M$ and vectors $\{u_{i}\}_{i=1}^{n}$, let $\mathcal{A}:=\{u_{i}u_{i}^{T}\}_{i=1}^{n}$, and note by Theorem 2.1 that $\lambda(\mathcal{A},p)=p^{T}Mp$. Now if $\min_{p\in\Delta_{n}}\lambda(\mathcal{A},p)$ can be approximated to an additive error of $n^{-\varepsilon}$ in polynomial time, then by Lemma 4.1 this yields a polynomial-time algorithm that approximates $\alpha(G)$ to a multiplicative factor of $n^{1-\varepsilon}O(\log n)$. But this is NP-hard by Theorem 4. ∎

4.2 PTAS for “well-conditioned” sets of rank-one matrices

Here, we note that there is a polynomial-time approximation scheme (PTAS) for approximating the design problem $\min_{p\in\Delta_{n}}\lambda(\mathcal{A},p)$ in the case where the set of rank-one matrices $\mathcal{A}$ is “well-conditioned.” (Note this does not conflict with the hardness results shown above, as explained below.)

The main idea is that since $\lambda(\mathcal{A},p)$ is a quadratic form in $p$ in this rank-one case by Theorem 2.1, we may use the well-known PTAS for minimizing polynomials over the simplex [12, Theorem 1.7]. Applying their PTAS to this design problem yields, for any fixed constant $\delta>0$, a polynomial-time algorithm that, given any collection $\mathcal{A}$ of rank-one matrices, outputs a distribution $\hat{p}\in\Delta_{n}$ satisfying

[TABLE]

Note that here “polynomial time” means time that is polynomial in the input size (namely $n$ and the bit complexity of $\mathcal{A}$), but not in $\delta$ (see the discussion in [12, Section 1.4]).

Intuitively, this approximation guarantee (4.7) is good whenever the quantity $\max_{p\in\Delta_{n}}\lambda(\mathcal{A},p)-\min_{p\in\Delta_{n}}\lambda(\mathcal{A},p)$ governing the error bound is not very large. One concrete setup where this provably occurs is when the set $\mathcal{A}=\{u_{i}u_{i}^{T}\}_{i=1}^{n}$ is “well-conditioned” in the sense that all pairwise angles $\measuredangle(u_{i},u_{j})$ are uniformly bounded away from $\tfrac{\pi}{2}$ – or in words, all pairs of $u_{i}$ and $u_{j}$ are “uniformly away” from being orthogonal. This is made precise in the following result. For simplicity of exposition, we assume that the matrices in $\mathcal{A}$ are symmetric (i.e. of the form $u_{i}u_{i}^{T}$) and normalized ($\|u_{i}\|_{2}=1$), but these assumptions can be easily removed.

Theorem 4.5**.**

Fix any $\delta,\gamma>0$. There is an algorithm that, given any $\mathcal{A}=\{u_{i}u_{i}^{T}\}_{i=1}^{n}$ such that each $\|u_{i}\|_{2}=1$ and $|\cos\measuredangle(u_{i},u_{j})|\geqslant\gamma$, outputs $\hat{p}\in\Delta_{n}$ satisfying

[TABLE]

in time that is polynomial in $n$ and the input size of $\mathcal{A}$ (but not necessarily in $\delta$ or $\gamma$).

We make use of the following simple bounds.

Lemma 4.2**.**

Let $\mathcal{A}=\{u_{i}u_{i}^{T}\}_{i=1}^{n}$. Then:

• •

$\max_{p\in\Delta_{n}}\lambda(\mathcal{A},p)=2\max_{i}\log\|u_{i}\|_{2}$.

• •

$\min_{p\in\Delta_{n}}\lambda(\mathcal{A},p)\geqslant\min_{i,j}\log|u_{i}^{T}u_{j}|$.

Proof.

By Theorem 2.1, $\lambda(\mathcal{A},p)=p^{T}M_{\mathcal{A}}p$ where $(M_{\mathcal{A}})_{ij}=\log|u_{i}^{T}u_{j}|$. For the maximization problem, note that $2\max_{i}\log\|u_{i}\|$ is achieved by $p=e_{i}$. Moreover, this value is optimal since for any $p\in\Delta_{n}$,

[TABLE]

where the final inequality is by the definition of $M_{\mathcal{A}}$ and Cauchy-Schwarz. The bound on the minimization problem follows similarly: for any $p\in\Delta_{n}$,

[TABLE]

∎

Proof of Theorem 4.5.

By Lemma 4.2, $\max_{p\in\Delta_{n}}\lambda(\mathcal{A},p)-\min_{p\in\Delta_{n}}\lambda(\mathcal{A},p)\leqslant\max_{i}\log\left(\|u_{i}\|_{2}^{2}\right)-\min_{ij}\log\left(|u_{i}^{T}u_{j}|\right)\leqslant\log\frac{1}{\gamma}$. The proof is complete by (4.7). ∎

Note that the PTAS in Theorem 4.5 is not at odds with the NP-hardness results proved in Theorems 4.1 and 4.4. This is because the runtime is not polynomial in the accuracy parameter $\delta$, and to obtain the approximation bounds discussed in the hardness results, $\delta$ must be taken to be scaling with $n$ – resulting in a runtime that is not polynomial in $n$.

5 Extension to ergodic stationary processes

So far in this paper, we have considered the analysis and design problems for the Lyapunov exponent in the standard setup where the random matrices are i.i.d. from a discrete distribution. In this section, we discuss extensions to the more general setup of ergodic stationary processes. As corollaries, we discuss setups including i.i.d. processes from continuous distributions, exchangeable processes, and Markov processes. This allows for the study of the stability of more exotic stochastic linear dynamical systems (1.2). The corresponding analysis and design problems are discussed in Subsections 5.1 and 5.2, respectively.

5.1 Ergodic formula for rank-one matrices

Here we show that, with only a minor modification, the rank-one Lyapunov exponent formula in Theorem 2.1 for i.i.d. processes extends to ergodic stationary processes. The proof is similar, the primary difference being the use of Birkhoff’s Ergodic Theorem.

Below, let $\mathcal{A}=\{uv^{T}:u,v\in\mathbb{R}^{d}\}$, and let it be equipped with the Borel $\sigma$-algebra (arising equivalently from any norm on $\mathbb{R}^{d\times d}$). We consider processes $A_{\sigma_{1}},A_{\sigma_{2}},\dots$ taking values in $\mathcal{A}$. Recall that this process is said to be stationary if $(A_{\sigma_{1}},\dots,A_{\sigma_{k}})$ has the same distribution as $(A_{\sigma_{\tau+1}},\dots,A_{\sigma_{\tau+k}})$ for every $k,\tau\in\mathbb{N}$. Recall also that this process is said to be ergodic if every shift-invariant event has probability [math] or $1$. We refer the reader to e.g., [19, Chapter 7] for further background on such processes.

Theorem 5.1**.**

*Let $A_{\sigma_{1}},A_{\sigma_{2}},\dots$ be an ergodic stationary process taking values in $\mathcal{A}$, and let $P$ denote the marginal distribution of $(A_{\sigma_{1}},A_{\sigma_{2}})$. If the integrability conditions $\mathbb{E}|\log\|A_{\sigma_{1}}\||<\infty$ and $\int\left\lvert\log\left\lvert u^{T}v^{\prime}\right\rvert\right\rvert dP(uv^{T},u^{\prime}v^{\prime T})<\infty$ are satisfied, then the (random) quantity *

[TABLE]

exists a.s., and moreover converges to

[TABLE]

a.s. and in $L_{1}$.

Note that the formula (5.1) is independent of the decompositions $u^{T}v$ of the rank-one matrices. This follows by a simple calculation that is identical to the one in Remark 2.1.

Proof.

Note that the theorem statement is independent of the norm $\|\cdot\|$ by equivalence of finite-dimensional norms; for convenience of computation, let us take it to be the spectral norm. Consider rank-one decompositions $A_{\sigma_{k}}=u_{\sigma_{k}}v_{\sigma_{k}}^{T}$ where w.l.o.g. we take $\|v_{\sigma_{k}}\|_{2}=1$. Just as in (2.2), expand the (random) quantity

[TABLE]

To deal with the final term $\frac{1}{k}\log\|u_{\sigma_{k}}\|_{2}$, observe that by stationarity and the first integrability assumption, we have $\mathbb{E}|\log\lVert{u_{\sigma_{k}}}\rVert_{2}|=\mathbb{E}|\log\lVert{u_{\sigma_{k}}}\rVert_{2}+\log\lVert{v_{\sigma_{k}}}\rVert_{2}|=\mathbb{E}|\log\lVert{A_{\sigma_{k}}}\rVert|=\mathbb{E}|\log\lVert{A_{\sigma_{1}}}\rVert|<\infty$. Thus the random variable $\log\|u_{\sigma_{k}}\|_{2}$ is bounded in $L_{1}$ (and thus also bounded a.s.), and so

[TABLE]

where the convergence is a.s. and in $L_{1}$. Now by ergodicity, stationarity, and the second integrability assumption, we may apply Birkhoff’s Ergodic Theorem (see e.g. [19, Theorem 7.2.1]) to obtain

[TABLE]

where the convergence is a.s. and in $L_{1}$. Combining (5.2), (5.3), and (5.4) completes the proof. ∎

We remark that the ergodic formula (5.1) in Theorem 5.1 depends on the process only through its second-order statistics. This occurs since the matrices are all rank-one. We now discuss several examples of ergodic stationary processes for which Theorem 5.1 applies.

5.1.1 I.i.d. processes

A first, simple example of ergodic stationary processes is an i.i.d. process. Specifically, consider the setup where $A_{\sigma_{1}},A_{\sigma_{2}},\dots$ are i.i.d. random variables with (the same) distribution $\mu$ over $\mathcal{A}$. Then Theorem 5.1 shows that, assuming some mild integrability conditions, we have the (a.s. and $L_{1}$) convergence:

[TABLE]

Note that this holds for continuous distributions $\mu$, i.e. distributions with infinite support. In the special case that $\mu$ has finite support, we recover Theorem 2.1.

Example 4**.**

Let $S^{d-1}=\{z\in\mathbb{R}^{d}:\|z\|_{2}=1\}$ denote the $d-1$ sphere, for $d>1$, and let $\mathsf{Unif}(S^{d-1})$ denote the uniform distribution over it (i.e. the $d-1$ dimensional spherical measure). In this example, we explicitly compute the Lyapunov exponent over “normalized” rank-one matrices $\mathcal{A}:=\{uv^{T}:u,v\in S^{d-1}\}$ for the following two probability measures:

• •

$\mu_{\text{sym}}$ is the induced distribution of $uu^{T}\in\mathcal{A}$, where $u\sim\mathsf{Unif}(S^{d-1})$.

• •

$\mu_{\text{asym}}$ is the induced distribution of $uv^{T}\in\mathcal{A}$, where $u,v\sim\mathsf{Unif}(S^{d-1})$ and are independent.

By Theorem 5.1, the Lyapunov exponent is the same for both $\mu_{\text{sym}}$ and $\mu_{\text{asym}}$, and is equal to:

[TABLE]

To compute this, we first compute the distribution $\nu$ of $u^{T}v$. By rotational invariance of the spherical measure $\mathsf{Unif}(S^{d-1})$, $\nu$ is equal to the distribution of $z_{1}$ where $z\sim\mathsf{Unif}(S^{d-1})$. This distribution is well-known to have density [41]

[TABLE]

where $\Gamma(\cdot)$ denotes the Gamma function. Using this, the integral in the right hand side of (5.6) – and thus the Lyapunov exponents for both $\mu_{\text{sym}}$ and $\mu_{\text{asym}}$ – can be computed to be:

[TABLE]

where $\Psi(\cdot)$ denotes the Digamma function and $\gamma\approx 0.5772$ denotes the Euler-Mascheroni constant.

The value of (5.8) is reported in Table 2 for small values $d\in\{2,\dots,8\}$. Note that the value for $d=3$ is equal to $-1$; this simple value arises because of the well-known fact that the distribution $\nu$ of the first coordinate of a random vector on the $2$-sphere, is the uniform distribution over $[-1,1]$ (this can be seen directly from (5.7)). Indeed, using this one can easily derive the special case $\mathbb{E}_{u,v\sim\mathsf{Unif}(S^{2})}\log|u^{T}v|=\operatorname{\frac{1}{2}}\int_{-1}^{1}\log|z|dz=\int_{0}^{1}\log zdz=(z\log z-z)_{0}^{1}=-1$.

The asymptotic behavior of (5.8) as the dimension $d\to\infty$ is

[TABLE]

due to the well-known inequality $\log z-z^{-1}\leqslant\Psi(z)\leqslant\log z-(2z)^{-1}$ that holds for all $z>0$ [3]. Note that the leading term $\log(\tfrac{1}{\sqrt{d}})$ makes sense intuitively since if $u,v\sim\mathsf{Unif}(S^{d-1})$ are independent, then $|u^{T}v|$ concentrates tightly around $\tfrac{1}{\sqrt{d}}$.666Here is a simple, standard way to derive this informally. The starting point is that if $\tilde{u}\sim\mathcal{N}(\mathbf{0},I)$, the normal distribution centered at the origin and with identity covariance, then $\tilde{u}/\|\tilde{u}\|_{2}\sim\mathsf{Unif}(S^{d-1})$ by rotational invariance. Thus $u^{T}v$ is equal in distribution to $\tilde{u}^{T}\tilde{v}/(\|\tilde{u}\|_{2}\|\tilde{v}\|_{2})$ where $u$ and $v$ are i.i.d. $\mathcal{N}(\mathbf{0},I)$. Now the numerator $\tilde{u}^{T}\tilde{v}$ is equal in distribution to $\tilde{v}_{1}\sim\mathcal{N}(0,1)$ by rotational invariance, and the denominator $\|\tilde{u}\|_{2}\|\tilde{v}\|_{2}$ concentrates tightly around $\sqrt{d}$ since $\|\tilde{u}\|_{2}^{2}=\sum_{i=1}^{d}\tilde{u}_{i}^{2}$ is $\Theta(\sqrt{d})$ by the Law of Large Numbers. Thus $|u^{T}v|$ is approximately equal to $|\tilde{v}_{1}|/\sqrt{d}$, where $\tilde{v}_{1}\sim\mathcal{N}(0,1)$, and this concentrates tightly around $\Theta(1/\sqrt{d})$. The dropped terms in this discussion are given precisely in the formula (5.8) and asymptotically by (5.9).

5.1.2 Exchangeable processes

Although exchangeable processes are not ergodic stationary processes777A simple example of an exchangeable process that is not ergodic is the process $X_{1},X_{2},\dots$ where all random variables take the same value $X$, and $X$ is a Bernoulli random variable., we can understand them in terms of such processes using De Finetti’s Theorem [11]. Specifically, consider the setup of an exchangeable process $A_{\sigma_{1}},A_{\sigma_{2}},\dots$ taking values in $\mathcal{A}$. Recall that (infinite) exchangeability means that for each $k\in\mathbb{N}$ and each permutation $\pi$ of $\{1,\dots,k\}$, $(A_{\sigma_{1}},\dots,A_{\sigma_{k}})$ has the same distribution as $(A_{\pi(\sigma_{1})},\dots,A_{\pi(\sigma_{k})})$. Now De Finetti’s Theorem – or, more precisely,  [28]’s generalization thereof to complete, separable metric spaces – states that the law of this process is equal to a mixture $\nu$ over i.i.d. product distributions over $\mathcal{A}$. It therefore follows from Theorem 5.1 that under the same mild integrability conditions, we have the (a.s. and $L_{1}$) convergence of

[TABLE]

5.1.3 Markov processes

Another example of an ergodic stationary process is a Markov process over a countable state space (for simplicity) with an irreducible transition kernel (for ergodicity) and with an initial distribution equal to the stationary distribution (for stationarity), see e.g. [19, Example 7.1.7]. That is, consider the setup where the non-initial matrices $A_{\sigma_{2}},A_{\sigma_{3}},\dots$ are drawn from an irreducible Markov kernel $Q$ over a countable subset of the rank-one matrices $\mathcal{A}$, and the initial matrix $A_{\sigma_{1}}$ is drawn from the stationary distribution $\pi_{Q}$ of the chain. Then Theorem 5.1 guarantees that, under mild integrability asumptions, we have the (a.s. and $L_{1}$) convergence

[TABLE]

Note that this recovers Theorem 2.1 since the original Lyapunov exponent setup is the special case where the columns of $Q$ are constant, i.e. $Q=\mathbf{1}p^{T}$, and the stationary distribution is $\pi_{Q}=p$. We also remark that (5.11) still holds when the assumption that the initial matrix $A_{\sigma_{1}}$ has stationary distribution $\pi_{Q}$ is relaxed.888This can be proved using a nearly identical argument as in the proof of Theorem 2.1.

5.2 Optimizing the Lyapunov exponent

We now turn to the design problem of optimizing the convergence rate of the stochastic linear dynamical system (1.2) over different types of stochastic processes $A_{\sigma_{1}},A_{\sigma_{2}},\dots$. We already showed in Subsection 4.1 that this design problem is NP-hard for i.i.d. processes. Below, we show that the design problem is also NP-hard for exchangeable processes (Subsection 5.2.1) yet admits a polynomial-time algorithm for Markov processes (Subsection 5.2.2).

5.2.1 Exchangeable processes

Here, we show that the design problem is NP-hard over exchangeable processes:

Theorem 5.2**.**

Given a collection $\mathcal{A}$ of $n$ symmetric rank-one matrices, it is NP-hard to decide whether there exists an exchangeable process $A_{\sigma_{1}},A_{\sigma_{2}},\dots$ taking values in $\mathcal{A}$, for which the stochastic linear dynamical system (1.2) converges (a.s.) for every initial state.

To give intuition for this result, recall the following three facts: (i) by De Finetti’s Theorem, exchangeable processes are mixtures of i.i.d. distributions; (ii) the Lyapunov exponent is a linear functional in this mixing measure (see (5.10)); and (iii) every linear program with a compact feasible set admits an optimal extreme point. Combining these three facts implies that, among exchangeable processes, there is an optimal i.i.d. process. This explains the hardness in Theorem 5.2, since the design problem over i.i.d. processes is NP-hard by Theorem 4.1. We now formalize this intuition:

Proof of Theorem 5.2.

By (5.10), there exists such a stabilizing exchangeable process if and only if there exists such a stabilizing i.i.d. process. But deciding whether there exists such a stabilizing i.i.d. process is NP-hard by Theorem 4.1. ∎

5.2.2 Markov processes

Here we show that the design problem is significantly easier when optimizing over irreducible Markov processes: there is a polynomial-time algorithm that not only can detect whether the associated dynamical system (1.2) is stabilizable, but also can compute an irreducible Markov process for which the convergence is fastest. This is in stark contrast to the i.i.d. and exchangeable setups, for which even the stabilization decision problem is NP-hard (see Theorems 4.1 and 5.2).

It will be helpful to introduce some notation. We say a square matrix is irreducible if for every simultaneous permutation of its rows and columns, the resulting matrix is a block-diagonal matrix with at most one diagonal block, not counting zero diagonal blocks999We note that the standard definition of an irreducible matrix does count zero diagonal blocks; this slight tweak allows us to more easily consider optimizing over subsets of $\mathcal{A}$ without having to introduce additional notation. A similar remark goes for our definition of irreducible Markov transition matrices.. We say $Q\in\mathbb{R}_{\geqslant 0}^{n\times n}$ is an irreducible Markov transition matrix on $n$ states if $Q$ is irreducible and $Q\mathbf{1}=\mathbf{1}$; let $\mathcal{Q}_{n}$ denote the set of all such matrices. Each $Q\in\mathcal{Q}_{n}$ has a unique stationary distribution which we will denote by $\pi_{Q}$.

For simplicity, we assume below an exact arithmetic model. Inexact arithmetic is easily handled by computing logarithms to arbitrary precision.

Theorem 5.3**.**

There is a polynomial-time algorithm that, given any collection of $n$ rank-one matrices $\mathcal{A}=\{u_{i}v_{i}^{T}\}_{i=1}^{n}$, outputs a minimizer of

[TABLE]

Note that the objective function in (5.12) is precisely the convergence rate of the corresponding dynamical system (1.2) by the discussion in Subsection 5.1.3. It therefore follows immediately that stability of the dynamical system can be detected in polynomial time: simply run the algorithm in Theorem 5.3 and check whether the optimization problem (5.12) has strictly negative value.

Corollary 5.1**.**

There is a polynomial-time algorithm that, given a collection $\mathcal{A}$ of $n$ rank-one matrices, decides whether there exists an irreducible Markov process $A_{\sigma_{1}},A_{\sigma_{2}},\dots$ taking values in $\mathcal{A}$ such that the stochastic linear dynamical system (1.2) converges (a.s.) for every initial state.

We now turn to the proof of Theorem 5.3. We proceed in two steps, informally sketched as follows. First, we observe that the optimization problem (5.12) can be equivalently recast as a linear program over circulations (Lemma 5.1). Then, by a basic result about decomposing circulations into cycles, we observe that this linear program is equivalent to the Minimum-Cycle-Mean problem in a certain associated graph (Lemma 5.2). The proof is then completed by recalling a well-known polynomial-time algorithm for the Minimum-Cycle-Mean problem.

We now formalize step $1$. Below, we say $F\in\mathbb{R}_{\geqslant 0}^{n\times n}$ is a (normalized) circulation of size $n$ if $F\mathbf{1}=F^{T}\mathbf{1}$ and $\sum_{ij}F_{ij}=1$; let $\mathcal{F}_{n}$ denote the set of all such matrices.

Lemma 5.1**.**

$\min_{Q\in\mathcal{Q}_{n}}\sum_{ij}(\pi_{Q})_{i}Q_{ij}\log|u_{i}^{T}v_{j}|=\min_{F\in\mathcal{F}_{n}}\sum_{ij}F_{ij}\log|u_{i}^{T}v_{j}|$.

Proof.

For any $Q\in\mathcal{Q}_{n}$, the matrix $F$ with entries $F_{ij}:=(\pi_{Q})_{i}Q_{ij}$ is easily checked to be in $\mathcal{F}_{n}$ since $F\mathbf{1}=F^{T}\mathbf{1}=\pi_{Q}$ by definition of $\pi_{Q}$ being the stationary distribution for $Q$. This proves the inequality “$\geqslant$”. For the other direction “$\leqslant$”, consider an optimal $F\in\mathcal{F}_{n}$ for the optimization problem on the right hand side. We may assume w.l.o.g. that $F$ is irreducible since a circulation can always be decomposed as a convex combination of irreducible circulations (this follows e.g. from [5, Exercise 7.14]), and each of these irreducible circulations must also be optimal by linearity of the cost function $F\mapsto\sum_{ij}F_{ij}\log|u_{i}^{T}v_{j}|$. Let $\pi:=F\mathbf{1}=F^{T}\mathbf{1}$, let $j^{*}$ be any fixed index in $\{1,\dots,n\}$ such that $\pi_{j^{*}}>0$ (such an index exists since $F\in\mathcal{F}_{n}$), and let $Q$ be the matrix with entries

[TABLE]

It is straightforward to check that $Q\in\mathcal{Q}_{n}$, that it has stationary distribution $\pi_{Q}=\pi$, and that $(\pi_{Q})_{i}Q_{ij}=F_{ij}$ for all $i,j$. This completes the proof. ∎

We now formalize step $2$. Below, $G_{\mathcal{A}}$ denotes the complete directed graph on $n$ nodes, one for each matrix $A_{i}$, with weight $\log|u_{i}^{T}v_{j}|$ on the edge from node $i$ to node $j$. The mean weight of a cycle in a $G_{\mathcal{A}}$ is the sum of the weights of the edges in the cycle, divided by the total number of edges. The minimum-cycle-mean of $G_{\mathcal{A}}$ is the minimum mean weight of a cycle, where the minimum is over all cycles in $G_{\mathcal{A}}$.

Lemma 5.2**.**

$\min_{F\in\mathcal{F}_{n}}\sum_{ij}F_{ij}\log|u_{i}^{T}v_{j}|$* is equal to the minimum-cycle-mean of $G_{\mathcal{A}}$. *

Proof.

Recall that a circulation (flow) can be decomposed as the convex combination of cycle (flows) see e.g. [5, Exercise 7.14]. Thus by linearity of the objective function $F\mapsto\sum_{ij}F_{ij}\log|u_{i}^{T}v_{j}|$, there exists an optimal circulation $F\in\mathcal{F}_{n}$ and a cycle $C$ such that $F_{ij}$ is equal to $1/|C|$ if $(i,j)$ is an edge in $C$, and [math] otherwise. The direction “$\geqslant$” then follows since $\sum_{ij}F_{ij}\log|u_{i}^{T}v_{j}|$ is equal to the mean weight of $C$. The other direction “$\leqslant$” follows by constructing $F\in\mathcal{F}_{n}$ from an optimal cycle $C$ in the same way. ∎

We are now ready to conclude the proof of Theorem 5.3.

Proof of Theorem 5.3.

By Lemmas 5.1 and 5.2, $\min_{Q\in\mathcal{Q}_{n}}\sum_{ij}(\pi_{Q})_{i}Q_{ij}\log|u_{i}^{T}v_{j}|$ is equal to the minimum-cycle-mean of $G_{\mathcal{A}}$. An optimal such cycle (i.e. one with minimum mean weight) can be computed in $O(n^{3})$ time [31]. Let $F$ denote the $n\times n$ matrix with entries $F_{ij}$ equal to $1/|C|$ if the cycle contains the edge $(i,j)$ and equal to [math] otherwise. Clearly $F\in\mathcal{F}_{n}$ and $\sum_{ij}F_{ij}\log|u_{i}^{T}v_{j}|$ is equal to the mean weight of $C$. Now construct $Q\in\mathcal{Q}_{n}$ from $F$ by (5.13) so that $\sum_{ij}F_{ij}\log|u_{i}^{T}v_{j}|=\sum_{ij}(\pi_{Q})_{i}Q_{ij}\log|u_{i}^{T}v_{j}|$. Thus $Q$ is an optimal solution for (5.12). Since $F$ and $Q$ are clearly computable from $C$ in polynomial time, this completes the proof. ∎

Remark 5.1**.**

It is known that the Joint Spectral Radius of a collection of rank-one matrices can be computed from the minimum-cycle-mean of $G_{\mathcal{A}}$ [2, Theorem 2.2]. We have shown here that the analogous optimization problem over irreducible Markov processes is also equivalent to the minimum-cycle-mean of $G_{\mathcal{A}}$.