A note on the eigenvectors of perturbed matrices with applications to linear positive systems

A result is presented describing the eigenvectors of a perturbed matrix, for a class of structured perturbations. One motivation for doing so is that positive eigenvectors of nonnegative, irreducible matrices are known to induce norms — acting much like Lyapunov functions — for linear positive systems, which may help estimate or control transient dynamics. The results apply to both discreteand continuous-time linear positive systems. The theory is illustrated with several examples.


Introduction
Perturbation theory describes the effects of an often unknown perturbation ∆ on the perturbed matrix A + ∆. It has applications in a broad range of disciplines including systems and control theory, engineering, numerical analysis and numerical linear algebra. In recognition of the importance (and intricacies) of perturbation theory, there are several monographs dedicated to its study including [47], [31] and [55]. Perturbation theory arises naturally in a context of linear dynamical systems -ubiquitous objects in the mathematical sciences. In discrete-time and finite-dimensions these take the form where A is an n × n matrix for some n ∈ N. Perturbations to a linear system replaces (1.1) by Equations of the form (1.2) are relevant when the original A in (1.1) is not reliably known, so that ∆ captures parametric or structural uncertainty in the dynamics. Alternatively, ∆ in (1.2) may denote an introduced or forced change to the dynamics. For instance, the choice of state-feedback u(t) := Kx(t) in the controlled linear system x(t + 1) = Ax(t) + Bu(t) , x(0) = x 0 , t ∈ N 0 , gives rise to (1.2) with ∆ = BK. For non-zero x 0 , the asymptotic behaviour of (1.1) and (1.2) is determined by the spectral radii r(A) and r(A + ∆), respectively. In the context of linear systems, a natural question is to ask how r(A + ∆) is determined, particularly when A in (1.2) is known, but ∆ is not. An answer to that question was given by Hinrichsen & Pritchard in [24,25] where the concept of the stability radius was developed in a control theoretic setting. The stability radius provides a notion of maximal, local robustness to perturbations in that when r(A) < 1, the stability radius is the largest β > 0 with the property that r(A + ∆) < 1 , ∀ ∆ such that ∆ < β .
In particular, if β is finite, then there exists a destabilising perturbation ∆ 0 with ∆ 0 = β and r(A + ∆ 0 ) = 1. We note that the stability radius is dependent on the norm that the set of perturbations is equipped with. We comment as well that the stability radius was developed independently in other contexts and is, in fact, an instance of the more general Wald's maximin model [54], see [44] for a helpful and interesting discussion. In numerous applications, the perturbations are structured, so that ∆ = BP C , B ∈ R n×m , C ∈ R p×n for some m, p ∈ N. (1. 3) The matrices B and C in (1.3) are known and describe the perturbation structure. The unknown m × p matrix P denotes the perturbation magnitude. In this set-up, the complex stability radius (when it is finite) is equal to the smallest P 0 over all P ∈ C m×p , where P 0 destabilises A + BP 0 C. When P is constrained to be real, that is P ∈ R m×p , then the resulting (minimal) norm of a destabilising perturbation is equal to the real stability radius (again, assuming that such a real perturbation exists). It is well-known that the complex and real stability radii need not coincide [22], see also [21]. Furthermore, whilst the complex stability radius is readily computable with a formula appearing in [25], computing the real stability radius is much more complicated in general, see [23,Section 5.3]. Suffice it to say, the stability radius is a ubiquitous and well-studied tool in robust control theory, a discipline predicated on tolerating uncertainty in controlled dynamical systems, such as in the sense of the ∆ that appears in (1.2). The stability radius is complemented with other concepts such as that of µ-values, introduced by Doyle in [8], see more recently, [23, or, for example, [30].
By using the Matrix Inversion Lemma, also known as the Sherman-Morrison-Woodbury formula [29, p. 19], it is possible to derive a relationship between the spectrum of the perturbed matrix A + BP C and the perturbation P in terms of the model data A, B and C. In the present note we exploit this relationship to provide formulae for the left and right eigenvectors of A + BP C. Despite their simplicity, and the wealth of existing knowledge on perturbation theory, we cannot find these results elsewhere in the literature. They demonstrate that determining eigenvectors of the n × n matrix A + BP C reduces to determining the eigenvectors of a min{m, p} dimensional matrix (the rank of the perturbation BP C), which in applications may be much smaller than n. The formulae are of interest in examples where eigenvectors are the desired object, such as the Google PageRank [37]. A second application, which we proceed to introduce and motivate, is that perturbed positive eigenvectors associated with a perturbed positive matrix specifying the perturbed linear system (1.2)-(1.3) induce norms that, in turn, provide estimates of the transient dynamics of (1.2)-(1.3).
Transient dynamics refers to short-term behaviour of a dynamical system, particularly dynamics which deviate away from those at steady-state. Even seemingly simple linear models (1.1) can exhibit exotic transient dynamics that are often neglected and overlooked and which may have serious implications in a myriad of physical contexts. Unlike asymptotic behaviour, transient dynamics of linear systems are not well captured by an eigenvalue (and thus by the stability radius either), as highlighted in the seminal work [53] related to the transition from laminar to turbulent flow in hydrodynamic stability theory. Instead, transition to turbulence occurs when the pseudospectrum crosses the stability/instability threshold. We refer the reader to [51] or [52] for more background and history on the development of the pseudospectrum. The term pseudospectrum is perhaps more commonly used by numerical analysts and the term spectral value sets, introduced in [19,20], is used by those in the control theory community for an equivalent concept.
and in this case t → V (x(t)) defined by C n ∋ x → V (x) := x, P x , where ·, · denotes the usual inner-product on C n , decreases along solutions of (1.1).
Dynamical systems that leave a positive cone invariant are called positive dynamical systems, or simply positive systems. There are different conventions in the academic literature regarding the use of the words "positive" and "nonnegative" in this context. The nonnegative orthant in R n , denoted R n + , with the usual partial ordering of componentwise inequality of vectors in R n is perhaps the most widely used positive cone in applications, for instance. Here invariance of the positive cone captures the essential feature that state-variables of positive systems, typically modelling abundances or concentrations, must be nonnegative. The study of positive systems is motivated by numerous models arising in a diverse range of fields from biology, chemistry, ecology and economics to genetics, medicine and engineering [17, p. xv There has subsequently been much attention devoted to generalisations of the Perron-Frobenius Theorem to nonlinear functions, including, for example [14,42,35]. Owing to their importance in applications there are several textbooks dedicated to the study of positive systems, such as [4,33,9], and to perturbations of positive semigroups [2]. In a context of perturbation theory, when A in (1.1) and B and C in (1.3) are all nonnegative, then the real and complex stability radii are known to be equal, and readily computable, see [27,26] (and [45] for continuous-time systems).
A stability radius approach to perturbation theory for linear positive dynamical systems and their transients has been considered by Hinrichsen & Plischke in [18]. There the authors note, as is well-known, that right and left eigenvectors (ensured by the celebrated Perron-Frobenius Theorem) associated with stable linear positive systems induce weighted infinity-and one-norms, respectively, which consequently act as Lyapunov functions. Using both these, and other candidate so-called Lyapunov vectors, estimates for transient dynamics of continuous time linear systems (1.4) are derived. The current work is similar in outlook to [18], although we comment that there the authors restrict attention to (exponentially) stable systems, which is not required presently. The Lyapunov functions induced by a left positive eigenvector of a nonnegative or Metzler (see below) matrix are used extensively elsewhere in the analysis and control of positive systems. Applications include, for example, robust control [38,39,5] and the control of switched-positive systems [12]. Moreover, these Lyapunov functions are examples of so-called max-separable and sum-separable Lyapunov functions [40] in the non-linear control literature see, for example, [7] and the references therein.
The results we derive apply to linear positive systems in continuous-time as well, specified bẏ Here the matrix M ∈ R n×n is Metzler (also known as essentially positive or quasi-positive), meaning that every off-diagonal entry is nonnegative and z 0 ∈ R n + . Recall that Metzler matrices characterise the matrices for which the matrix exponential e M t for t ≥ 0 -and hence the solution of (1.4) -is componentwise nonnegative.
Our main result of this note is Theorem 2.1, stated and proven in Section 2, which provides formulae for eigenvectors of matrices subject to structured perturbations. The formulae for the eigenvectors do not require that the matrices are nonnegative, and we present them for full generality. Section 3 applies these results to Lyapunov functions for perturbed linear systems. We illustrate the theory with examples in Section 4 and Section 5 contains a brief discussion.
Notation: Most notation we use is standard or is defined as it is introduced. The symbols N, R and C denote the sets of positive integers, real numbers and complex numbers, respectively. The symbol N 0 denotes the set of nonnegative integers. For n ∈ N, we let n := {1, 2, . . . , n}.
We recall that a square matrix A ∈ R n×n with entries a ij is said to be reducible if there exist non-empty, disjoint subsets J 1 , J 2 ⊆ n such that J 1 ∪ J 2 = n and a ij = 0 for all (i, j) ∈ J 1 × J 2 . If A is not reducible then it is said to be irreducible. If A is additionally nonnegative, then A is irreducible if, and only if, for each i, j ∈ n there exists k ∈ N such that the (i, j)-th entry of A k is positive. We write that a vector or matrix is strictly positive if every component is positive. We let σ(A) and r(A) denote the spectrum of A and spectral radius of A, respectively.

Eigenvectors of matrices with structured perturbations
defined for all z ∈ C that are not poles of G. For any P ∈ R m×p and λ ∈ C (a) λ ∈ σ(A) satisfies λ ∈ σ(A + BP C) if, and only if, 1 ∈ σ(G(λ)P ). The geometric multiplicity of λ as an eigenvalue of A + BP C is equal to that of one as an eigenvalue of G(λ)P or P G(λ). Remark 2.2. Converses to parts (b) and (c) of the above theorem hold, which we describe here for completeness and are in fact used to prove the assertion in (a) regarding geometric multiplicities. Namely, if θ T is a left eigenvector of A + BP C corresponding to the eigenvalue λ ∈ σ(A), then θ T BP and θ T B are left eigenvectors of G(λ)P and P G(λ) corresponding to the eigenvalue one, respectively. Similarly, if ω is a right eigenvector of A + BP C corresponding to the eigenvalue λ ∈ σ(A), then Cω and P Cω are right eigenvectors of G(λ)P and P G(λ) corresponding to the eigenvalue one, respectively.
Proof of Theorem 2.1: The first part of assertion (a) is known and follows from [22, Proposition 2.3] (alternatively compare with, for example, [28,Theorem 4.3]). The second part shall be derived in the proof of (b), which we proceed to next. We let g.m. M (µ) denote the geometric multiplicity of the eigenvalue µ of the square matrix M and recall that 1 ∈ σ(G(λ)P ) if, and only if, 1 ∈ σ(P G(λ)).
If ζ T is a left eigenvector of P G(λ) corresponding to the eigenvalue one, that is as required. A similar calculation to that in (2.3) demonstrates that if ζ T j are linearly independent left eigenvectors of P G(λ) corresponding to the eigenvalue one, then ζ T j P C(λI − A) −1 are linearly independent, and thus g.
Now suppose that g.m. A+BP C (λ) = r and let θ T 1 , . . . , θ T r denote r linearly independent left eigenvectors of A + BP C, corresponding to the eigenvalue λ, so that Since λ ∈ σ(A), it follows from (2.7) that both θ T j B = 0 and θ T j BP = 0. Rearranging (2.7) we obtain which establishes that θ T j BP and θ T j B are left eigenvectors of G(λ)P and P G(λ), respectively. If there exists constants ω 1 , . . . , ω r ∈ C such that where we have used (2.8), and since the θ T j are linearly independent. We conclude that the r vectors θ T j BP are linearly independent, as are the r vectors θ T j B, whence Combining (2.4), (2.6) and (2.10) we see that and hence the claimed one-to-one correspondence between eigenvectors holds. The proof of statements (a) and (b) is now complete.
The proof of (c) is analogous to that of (b) and is therefore omitted. Parts (d) and (e) now follow from (b) and (c), respectively, noting that as A + BP C is nonnegative and irreducible, the Perron-Frobenius Theorem implies that λ = r(A + BP C) is a simple eigenvalue. Hence, the strictly positive left and right eigenvectors of A + BP C corresponding to λ are uniquely given (up to multiplication by a positive scalar) by ξ T C(λI − A) −1 and (λI − A) −1 BP η, for ξ T and η denoting left and right eigenvectors of G(λ)P , respectively, corresponding to the simple eigenvalue one (necessarily simple by part (a)).
The formulae derived in Theorem 2.1 may be of use in applications where eigenvectors are the desired object, as we seek to illustrate in Example 4.2. A second application of the theorem, that we pursue in Section 3, is that (perturbed) strictly positive eigenvectors induce norms which may be used to provide estimates of the transient dynamics of (perturbed) linear positive systems. Before that we provide some remarks on the above theorem. Remark 2.3. (i) The relationship in Theorem 2.1 (a) between eigenvalues λ of A + BP C, and one being an eigenvalue of G(λ)P , holds for all n eigenvalues of A + BP C, not just its spectral radius. When A+BP C is nonnegative, determining the difference between λ ∈ σ(A+BP C) and λ = r(A+BP C) is not always immediate. We draw attention to [36,Theorem 2.2], however, for rank-one perturbations which provides sufficient conditions for when λ ∈ σ(A + bpc T ) does imply that λ = r(A + bpc T ).
(ii) If, in addition to the assumptions of Theorem 2.1 parts (d) and (e), A is assumed irreducible and 0 ≤ BP C = 0 then, from [3, p. 27] it follows that and thus λ = r(A + BP C) implies that λ ∈ σ(A).

(iii) If the assumption that A + BP C is irreducible in Theorem 2.1 for parts (d) and (e) is relaxed, and
A + BP C is instead there only assumed to be nonnegative, then, although r(A + BP C) is still a (nonnegative) eigenvalue of A + BP C, the associated nonnegative eigenvector (there may be others) need not be strictly positive in general. Strict positivity of eigenvectors is essential for the norms described in Section 3, although we comment on the reducible case in Section 3.3.
Determining the eigenvectors of A + BP C is an a priori operation on n × n matrices. Theorem 2.1 demonstrates that it is, in fact, an min{m, p} dimensional problem -the rank of the perturbation BP C -which may be much smaller than n. To demonstrate this property, we present as a corollary the case where B = b ∈ R n ; the case where C = c T ∈ R 1×n is analogous. The condition 1 = q T G(λ) in Corollary 2.4 (a) is a scalar equation in at most p + 1 unknowns -the p entries of q T and λ. Therefore, computing the eigenvectors of G(λ)q T or q T G(λ) is not required in this case.

Perturbations of norms and Lyapunov functions for linear positive systems
Here we recall how eigenvectors of the nonnnegative, irreducible matrix A induce norms that help describe the dynamics of the linear positive dynamical system (1.1). We first consider the case where the matrices involved are irreducible, in discrete-and continuous-time in Sections 3.1 and 3.2, respectively, before commenting on the reducible case in Section 3.3.

Discrete-time
Given irreducible A ∈ R n×n + , let v T and w denote strictly positive left and right eigenvectors of A corresponding to the eigenvalue r = r(A) of A, which exist and are unique up to multiplication by a positive scalar by the Perron-Frobenius Theorem. Let V, W : R n + → R + denote the norms on the positive cone R n + given by where w i denotes the (necessarily positive) i-th component of w, for i ∈ n. Note that V is linear, whilst W is not. Both V and W are only determined up to multiplication by a positive scalar, which may be fixed by fixing the norm of v T and w. In words, V and W are weighted one-and infinity-norms on R n + , respectively, and satisfy the norm equivalences: As mentioned in the introduction, the functions V and W appear elsewhere in the positive systems literature, with different terminology. For instance, the function V is called a linear copositive Lyapunov function in [5] and is an example of a sum-separable Lyapunov function from [40]. The function W is an example of a max-separable Lyapunov function, also from [40].
Along solutions x of the discrete-time linear positive system (1.1), V satisfies the equalities Similarly, arguing now for W , we estimate for t ∈ N 0 4) and comment that the inequality in (3.4) is an equality if x(t) = w, for t ∈ N 0 . The inequality (3.4) admits the estimate W (x(t)) ≤ r t W (x 0 ) , t ∈ N 0 .  1). However, the following observations hold for any r > 0. For large t, the right hand sides of (3.3) and (3.5) are dominated by r t , determining the asymptotic dynamics of the linear system (1.1) (unless x 0 = 0). For small t, the terms V (x 0 ) and W (x 0 ) play a larger role in determining V (x(t)) and W (x(t)), respectively, and may capture the transients of x(t), in either a weighted one-or infinity-norm, respectively. For fixed c 1 , c 2 > 0 the equations V (x 0 ) = c 1 and W (z 0 ) = c 2 describe n-hyperplanes and n-hyperrectangles, respectively, intersected with the positive cone R n + in the unknowns x 0 , z 0 ∈ R n + .
Combining the calculations of this section with Theorem 2.1 yields the obvious corollary.
Corollary 3.2. Given A ∈ R n×n + , B ∈ R n×m and C ∈ R p×n , assume that P ∈ R m×p and r ≥ 0 are such that (a) A + BP C ∈ R n×n + is irreducible; let ξ T , ζ T , ν and η be as in Theorem 2.1. Then V 1 , V 2 , W 1 , W 2 : R n + → R + defined by are norms on R n + (once the eigenvectors from Theorem 2.1 are chosen to have positive components). Further, the solutions x of the perturbed, discrete-time linear positive system When r(A) = 1 and A is primitive then the function V in (3.1) provides estimates of the one-norm of the asymptotic solution x of (1.1). Particularly, when A is primitive it is well-known that for any r(A) the solution x of (1.1) satisfies where v T and w are positive left and right eigenvectors of A corresponding to r(A), respectively. Note that the right hand side of (3.8) is independent of the scalings of v T and w chosen. Thus, if v T , w are (uniquely) chosen so that w 1 = 1 and v T w = 1 , (3.9) then (3.8) may be rewritten as lim where V is as in (3.1). Therefore, in the situation that r(A) = 1 it follows that that is, the solution x(t) converges in one-norm to V (x 0 ), and is asymptotically parallel to w. Clearly, lim t→∞ x(t) 1 depends linearly on x 0 1 and hence over all initial states with one-norm one we have and min for some k, ℓ ∈ n. Moreover, the above maximum and minimum are respectively attained at x 0 = e k and x 0 = e ℓ . By replacing A and V above by A + BP C and V 1 or V 2 from (3.6), respectively, and appealing to Theorem 2.1 and Corollary 3.2, the same comments apply to the solution of the perturbed linear system (3.7). The equation (

Continuous-time
Here we present parallel results to those in Sections 2 and 3.1 only now for continuous-time linear positive systems (1.4), that is,ż (t) = M z(t) , z(0) = z 0 , t ∈ R + .
In (1.4), z 0 ∈ R n + for n ∈ N and M ∈ R n×n is Metzler, meaning that M ij ≥ 0 for all i, j ∈ n such that i = j. Noting that the unique solution of (1.4) is given by where R + ∋ t → e M t denotes the usual matrix exponential, it is well-known that the family of semigroups (e M t ) t≥0 is componentwise nonnegative for all t ∈ R + if, and only if, M is Metzler, see, for example [43, Section 3.1] or [41, Theorem 3] for a proof.
In continuous-time the asymptotic dynamics of (1.4) (certainly for z 0 = 0) are determined by the spectral abscissa of M , which is defined as The following result is an analogue of the Perron-Frobenius Theorem for irreducible Metzler matrices, and is well-known. (1) a ∈ σ(M ) and a = r(µI + M ) − µ; (2) if λ ∈ σ(M ) and λ = a, then Re λ < a.
Furthermore, under the additional assumption that M is irreducible, the following statements hold.  is differentiable and moreover The differential equation (3.12) has solution Arguing now for t → W (w(t)) (which need not be classically differentiable), we use the series expansion of the matrix exponential and (3.11) to see that (3.14) Therefore, for t ∈ R + W (z(t)) = max where we have used (3.14). The inequality in (3.15) is an equality if z 0 = w. Combining the calculations of this section with Theorem 3.5 yields the obvious corollary.
Corollary 3.7. Given M ∈ R n×n + , B ∈ R n×m and C ∈ R p×n , assume that P ∈ R m×p and a ∈ R are such that (a) M + BP C ∈ R n×n + is Metzler and irreducible; let ξ T , ζ T , ν and η be as in Theorem 3.5

The reducible case
The material in Section 3 thus far has assumed irreducibility to ensure that V and W in (3.1) (as well as V i and W i in Corollaries 3.2 and 3.7) are well-defined norms. The irreducibility assumption may be dropped at the potential loss of accuracy and increased conservatism, as we proceed to describe.
Given A ∈ R n×n + , it follows that A ε := A + εQ is strictly positive and hence irreducible, for all ε > 0 and strictly positive Q ∈ R n×n + . The spectral radius of A is continuous with respect to the entries of A and so with r ε := r(A + εQ) and r = r(A), then 0 ≤ r ε − r → 0 as ε → 0. Letting v T ε and w ε denote strictly positive left and right eigenvectors of A ε corresponding to r ε , respectively, we may define V ε and W ε as in (3.1). Straightforward adjustments to the calculations (3.2) and (3.4) demonstrate that along solutions x of (1.1) V ε (x(t + 1)) ≤ r ε V ε (x(t)) and W ε (x(t + 1)) ≤ r ε W ε (x(t)) t ∈ N 0 . (3.17) If r(A) < 1, then ε > 0 and strictly positive Q ∈ R n×n + may be chosen so that r ε < 1, and so (3.17) demonstrates that V ε and W ε are Lyapunov functions for (1.1) -in other words, V ε and W ε still capture the same qualitative behaviour of solutions of (1.1). Corollary 3.2 naturally extends to the reducible case by replacing A with A ε . When r ≥ 1, however, the inequalities in (3.17) (compare with the equality in (3.2)) are less informative and may be conservative. Further, recall that our motivating application for Corollary 3.2 is to see how structured perturbations affect transient dynamics in linear positive systems, that is, how the transients of (3.7) compare to those of (1.1). The solution that we propose uses the norms V i and W i in (3.6). The introduction of ε > 0 and strictly positive Q ∈ R n×n + obfuscates the relative contributions of εQ and the perturbation BP C to (V ε ) i and (W ε ) i .
A different perspective is that irreducibility of A ensures that every (non-zero) solution of the linear positive dynamical system (1.1) experiences the same asymptotic rate of growth or decline (in the norm V ), captured by (3.2). This qualitative property of independence of asymptotic rates of growth on initial conditions -reminiscent of the independence of initial distributions on the limiting distribution of ergodic Markov chains, see [13,Theorem 11',p.95] -need not hold for reducible matrices. Consider the simple example of a reducible A given by where A 1 and A 3 are both irreducible. Since A leaves a proper cone invariant if, for example, r(A 1 ) < 1 and r(A 3 ) > 1, then clearly, the solutions of (1.1) from the initial states x 0 0 and 0 x 0 shall exhibit different asymptotic rates of growth. Although adding εQ to A would make the material of Section 3 applicable, it is arguably qualitatively more appropriate to instead consider the eigenvectors and norms induced by the irreducible components A 1 and A 3 .
The above comments are also applicable in the continuous-time case, by noting that A ε := A + εQ is Metzler if A is, ε > 0 and Q ∈ R n×n + is strictly positive.

Examples
Example 4.1. Block-wise matrix inversion -the Matrix Inversion Lemma -states that for A ∈ C n×n , b, c ∈ C and d ∈ C also known as the Abcd lemma, provided that A −1 exists and the Schur complement of A, G(0) = d − c T A −1 bc = 0. The expression (4.1) is useful because it determines M −1 from lower dimensional quantities and is particularly effective when computing A −1 is elementary (or A −1 is already known). We seek formulae for the eigenvectors of M in terms of A, b, c T and d.
By decomposing M into the sum of a matrix and a rank-one perturbation where G is as in (2.1). For such λ ∈ σ(M ) the corresponding left and right eigenvectors of M are given by by (4.2), and respectively. We note that the same formulae as (4.2)-(4.4) (the latter two up to multiplication by a positive scalar) are obtained when M is instead decomposed as The above arguments do not require that M is nonnegative. If M is nonnegative then by [36,Theorem 2.1], it follows that M may only have one positive eigenvalue greater than that of r(A), which therefore must necessarily equal r(M ).  1). Examples include simple random walks with absorbing or reflecting boundaries or simple queuing models, as well as applications in group theory in [1]. Amongst numerous suitable monographs on stochastic processes we refer the reader to, for example, [10,Ch. 16] or [34,Ch. 1]. Let X t denote the state of a discrete-time (time-homogeneous or stationary) Markov chain, for t ∈ N 0 , taking values in {1, 2, . . . , n} with associated transition probabilities Defining A ∈ R n×n + with (i, j) th entry a ij it follows that A is a (left) stochastic matrix (also termed a probability matrix, transition matrix, substitution matrix or Markov matrix) since n i=1 a ij = 1 , ∀ j ∈ {1, 2, . . . , n}. (4.6) In words, (4.6) states that every column sum of A is one, and thus as r(A) ∈ σ(A) However, evidently (4.6) also yields that for any x 0 ∈ R n + with x 0 1 = 1 imply that, as is well-known, that w is the limiting, stationary distribution of (X t ) t∈N0 .
Here we use Theorem 2.1 to give an example of how the stationary distribution changes analytically with perturbations to the transition matrix. Consider a Markov chain with n = 3 states, transition matrix T and stationary distribution w 0 : Any perturbation to T must preserve the properties (4.5) and (4.6) and we consider the two-parameter perturbation   1 2 where we note that r(A) < 1. By construction each column sum of the matrix in (4.8) is equal to one, and thus we expect that 1 = r(A+BP C) (with B = I) for any p, q as in (4.8). For consistency we demonstrate that condition Theorem 2.1 (a) that 1 ∈ σ(A + BP C) if, and only if, 1 ∈ σ(C(I − A) −1 P ) = σ(G(1)P ), always holds. Indeed, an elementary calculation shows that We now seek the right eigenvector of T corresponding to r(T ) = 1. A right eigenvector η of G(1)P corresponding to the eigenvalue one is given by and hence, by Theorem 2.1 (c), w = (I − A) −1 BP η = (I − A) −1 P η is a simple right eigenvector of T corresponding to r(T ) = 1. Once normalised so that w 1 = 1 we find that (4.9) Correctly, we see from (4.9) that when p = q = 0 we recover w(0, 0) = w 0 in (4.7). As the number of parameters is low, the components of w have been graphed over (p, q)-parameter space in Fig. 1.
We conclude this example by noting that as the state-dimension n = 3 is low, the right eigenvalue in (4.9) could be derived directly from T in (4.8), either by hand or through symbolic computing. We have chosen n = 3 so as to demonstrate the concepts involved without obscuring them with calculations -the real value of Theorem 2.1 and Theorem 3.5 is that they apply when n is very large, and computing w in (4.9) directly is computationally expensive.
Example 4.3. Matrix population projection models are linear systems of the form (1.1) and are a tool for modelling stage-structured populations. They are used broadly from conservation and harvesting to evolutionary theory and we refer the reader to the monograph of Caswell [6] for further background. A reasonable assumption for meaningful ecological models (see [48]) is that A is primitive, so that (3.8) holds, that is, Here w and v T are positive right and left eigenvectors of A, respectively, corresponding to r(A) and are typically called the stable-stage structure and the reproductive vector, respectively. The latter receives its name as v T contains the reproductive values (as in [11]) of each stage-class [15]. The product v T x 0 that appears in the numerator of (3.8) contains the contributions to the asymptotic population from the initial population distribution x 0 . The nonnegative constant that appears on the right hand side of (3.8) is defined as the population inertia [32] of A from x 0 . Noting that the population inertia from x 0 = w is one; the population inertia of A from arbitrary x 0 is a long term multiplicative ratio of the size of the population projected from y 0 compared to that projected from stable stage-structure w. Combining (3.8) with Theorem 2.1 for the perturbed eigenvectors demonstrates that is equal to the population inertia of A + BP C from x 0 , where ξ T , ζ T , η and ν are as in Theorem 2.1. The expressions in (4.10) are the expected generalisations of the rank-one, one-parameter perturbations presented in [49]. Example 4.4. Our final example is based on [36, Section 3.1] which considers the effects of three parameters on the spectral radius of matrix population projection model from [50] for the invasive weed Cirsium vulgare. The linear model has four stage-classes denoting the seed bank, small, medium and large weeds, respectively, and is given by: where the three parameters of interest are the germination rate g, summer survival s of small plants and mortality caused by floral herbivory h. As survival/mortality rates the parameters satisfy 0 ≤ s, g, h ≤ 1 and A is easily observed to be primitive whenever 0 < s, g, h < 1. The nominal values g nom = 0.2142, s nom = 0.516 and h nom = 0.942 from [50] give rise to r(A) = 1.58 > 1 where A := A(s nom , g nom , h nom ).
To write A as a structured perturbation of A depending on the parameters s, g and h we write  51. In [36] the authors seek to describe the required changes to g, s or h that asymptotically stabilise the weed population, that is, give rise to r(A) = 1. The r(A) = 1 surface is plotted in (g, s, h)-parameter space in Fig. 2(a) and has been found from Theorem 2.1 (a) by varying s and g in the interval [0.1, 0.9] and determining h by solving det(I − P G(r(A))) = det(I − P G(1)) = 0 .
Here G(1) = C(I − A) −1 B = (I − A) −1 B, as C = I in (4.11). The surface in Fig. 2(a) is a reproduction of [36, Fig. 3.1]. Any (g, s, h)-triple lying on the surface results in r(A) = 1. However, the analysis in [36] does not consider the impact of the any potential management strategy (that is, perturbation) on the resulting dynamics -particularly x(t) 1 as t → ∞.
Supposing that (s * , g * , h * ) are such that r(A) = 1 then, as in Remark 3.3 (iii), it follows that as in (3.9). In (4.12), ζ T and ν are any right and left eigenvectors of P G(1) corresponding to the eigenvalue one. In addition to the r(A) = 1 surface in Fig. 2(a), contours of max V 2 (z) have been overlaid (as well as plotted in the h = 0.8 plane). Combined, it is possible to inspect how both r(A) and max V 2 (z) vary over the parameter space. In Fig. 2(b), the surface max V 2 (z) has been plotted. In the present applied context of managing an invasive weed, perturbations that result in r(A) ≤ 1 are desirable. We note, however, that although small values of s, g require smaller values of h to lead to r(A) = 1, they also result in both a larger transient and asymptotic population abundance. Finally, 20 projections of (1.1) with (s, g, h) = (0.79, 0.58, 0.9954) are plotted in Fig. 2(c) from 18 random initial conditions. Additionally, the two lines initial conditions to x 0 1 = e 1 and x 0 4 = e 4 are plotted so that

Conclusion
We have considered structured perturbations of matrices and especially the co-dependencies of eigenvalues and eigenvectors. The perturbed eigenvalues capture the asymptotic behaviour of the perturbed time-invariant linear systems (1.2) or (3.16), and here we have followed existing stability radius arguments. Meanwhile, the perturbed eigenvectors capture other, example specific, features of the perturbed model. When the perturbed matrix has a fixed spectrum, for example as a consequence of "poleplacement" (also known as pole-shifting, see [46,Ch. 5]) or because of other constraints such as a need to preserve column stochasticty in matrix models of Markov chains, we find that the perturbed left and right eigenvectors themselves are then highly constrained by the structure of the perturbation. When the dimension of the matrix model is high but the rank of the perturbation is low, then the dimension of the problem of finding eigenvectors is significantly reduced. Similarly, when the number of parameters defining the perturbation is low (typically no greater than three) then the co-dependencies of eigenvalues and (entries of the) eigenvectors may be displayed graphically. We believe that our results are useful in a variety of contexts: they can be used to investigate the stable stage structure in ecological or population models, equivalently the stationary distribution in perturbed Markov Chains; to study the dependence of population inertia on vital rates, and; to consider how linear Lyapunov functions, as determined by left and right eigenvectors, respond to perturbations. These applications are considered through four simple examples and the connection to estimating or controlling transient dynamics via induced norms (or Lyapunov functions) is complementary to the results of [18].