Iterative building of Barabanov norms and computation of the joint spectral radius for matrix sets

The problem of construction of Barabanov norms for analysis of properties of the joint (generalized) spectral radius of matrix sets has been discussed in a number of publications. The method of Barabanov norms was the key instrument in disproving the Lagarias-Wang Finiteness Conjecture. The related constructions were essentially based on the study of the geometrical properties of the unit balls of some specific Barabanov norms. In this context the situation when one fails to find among current publications any detailed analysis of the geometrical properties of the unit balls of Barabanov norms looks a bit paradoxical. Partially this is explained by the fact that Barabanov norms are defined nonconstructively, by an implicit procedure. So, even in simplest cases it is very difficult to visualize the shape of their unit balls. The present work may be treated as the first step to make up this deficiency. In the paper two iteration procedure are considered that allow to build numerically Barabanov norms for the irreducible matrix sets and simultaneously to compute the joint spectral radius of these sets.

where A , for a matrix A, is the matrix norm generated by the vector norm · in R m , that is A = sup x =1 Ax . Then the limit ρ(A ) := lim sup n→∞ (ρ n (A )) 1/n (1) does not depend on the choice of the norm · and is called the joint spectral radius of the matrix set A [33]. When r = 1 this definition coincides with the famous Gelfand formula for the spectral radius of a matrix [13] as in this case A = {A} is a singleton matrix set andρ n (A ) = A n .
For matrix sets A consisting of a finite amount of matrices, as is our case, the quantitiesρ(A ) andρ(A ) coincide with each other [5] and their common value is denoted as This last formula may serve as a basis for a posteriori estimating the accuracy of computation of ρ(A ). The first algorithms of a kind in the context of control theory problems have been suggested in [6], for linear inclusions in [2], and for problems of wavelet theory in [8][9][10]. Later the computational efficiency of these algorithms was essentially improved in [14,26]. Unfortunately, the common feature of all such algorithms is that they do not provide any bounds for the number of computational steps required to get desired accuracy of approximation of ρ(A ). Recently, in [19] explicit computable estimates for the rate of convergence of the quantities A n 1/n to ρ(A) end their extension to the case of the joint spectral radius were obtained. Probably, these results will help to make more constructive the problem of evaluating of ρ(A ) by the generalized Gelfand formula (1) .
Some works suggest different formulas to compute ρ(A ). So, in [7] it is shown that ρ(A ) = lim sup n→∞ max Ai j ∈A |tr(A in · · · A i2 A i1 )| 1/n , where, as usual, tr(·) denotes the trace of a matrix. Given a norm · in R m , denote Then the spectral radius of the matrix set A can be defined by the equality where infimum is taken over all norms in R d [12,33], and therefore ρ(A ) ≤ A for any norm · in R m . For irreducible matrix sets, 1 the infimum in (2) is attained, and for such matrix sets there are norms · in R m , called extremal norms, for which A ≤ ρ(A ).
In analysis of the joint spectral radius ideas suggested by N.E. Barabanov [2][3][4] play an important role. These ideas have got further development in a variety of publications among which we would like to distinguish [35]. Barabanov). Let the matrix set A = {A 1 , . . . , A r } be irreducible. Then the quantity ρ is the joint (generalized) spectral radius of the set A iff there is a norm · in R m such that Throughout the paper, a norm satisfying (4) will be called a Barabanov norm corresponding to the matrix set A . Note that Barabanov norms are not unique.
Similarly, as is shown in [31,Thm 3.3] and [32], the value of ρ equals to ρ(A ) if and only if for some central-symmetric convex body 2 S the following equality holds where conv(·) denotes the convex hull of a set and ρS := {ρx : x ∈ S}. As is noted in [31], the relation (5) was proved by A.N. Dranishnikov and S.V. Konyagin, so it is natural to call the central-symmetric set S the Dranishnikov-Konyagin-Protasov set. The set S can be treated as the unit ball of some norm · in R d (recently this norm is usually called the Protasov norm). Note that Barabanov and Protasov norms are the extremal norms, that is they satisfy the inequality (3). In [28,29,36] it is shown that Barabanov and Protasov norms are dual to each other. Remark that formulas (3), (4) and (5) define the joint or generalized spectral radius for a matrix set in an apparently computationally nonconstructive manner. In spite of that, namely such formulas underlie quite a number of theoretical constructions (see, e.g., [1,18,21,27,35,36]) and algorithms [30] for computation of ρ(A ).
Different approaches for constructing Barabanov norms to analyze properties of the joint (generalized) spectral radius are discussed, e.g., in [15,17] and [34,Sect. 6.6]. In [18,21] the method of Barabanov norms was the key instrument in disproving the Lagarias-Wang Finiteness Conjecture. The related constructions were essentially based on the study of the geometrical properties of the unit balls of some specific Barabanov norms. In [16,17] the method of extremal polytope norms was the key tool in investigation of the finiteness properties of pairs of 2 × 2 matrices.
In this context the situation when one fails to find among current publications any detailed analysis of the geometrical properties of the unit balls of Barabanov norms looks a bit paradoxical. Partially this is explained by the fact that Barabanov norms are defined nonconstructively, by an implicit procedure. So, even in simplest cases it is very difficult to visualize the shape of their unit balls. The present work may be treated as the first step to make up this deficiency.
In the paper, an iteration procedure is considered that allows to build numerically Barabanov norms for the irreducible matrix sets and simultaneously to compute the joint spectral radius of these sets. A similar iteration procedure is also discussed in [22].
The structure of the paper is as follows. In Introduction we give basic definitions and present the motivation of the work. In Section 2 the main iteration procedure is introduced and Main Theorem stating convergence of this procedure is formulated. The iteration procedure under consideration is called the max-relaxation procedure since in it the next approximation to the Barabanov norm is constructed as the maximum of the current approximation and some auxiliary norm. In Section 3 proof Main Theorem is given. In Section 4, to build simplest examples, the max-relaxation scheme is adapted for computations with 2 × 2 matrices. Results of numerical tests are illustrated by two examples. At last, in concluding Section 5 we discuss some shortcomings of the proposed approach and formulate further unresolved problems.
Throughout the paper, a continuous function γ(t, s) defined for t, s > 0 and satisfying min{t, s} < γ(t, s) < max{t, s} for t = s, will be called an averaging function. Examples for averaging functions are: Given some averaging function γ(·, ·), a norm · 0 in R m , and a vector e ∈ R m such that e 0 = 1, construct recursively the norms · n and · • n , n = 1, 2, . . ., in accordance with the following rules: MR1: if the norm · n has been already defined compute the quantities MR2: define the norms · n+1 and · • n+1 : Remark that the number of operations needed to perform one step of algorithm MR1-MR2 is of order rm 2 ν(ε), where ν(ε) is the number of operations needed to compute, for an arbitrary vector x ∈ R m , the value of the norm x with a relative accuracy ε. In general, the value ν(ε) is of order ε −m . So, the total number of operations needed to perform n steps of algorithm MR1-MR2 has the same rate of growth as nrm 2 ε −m .
Remark also, that the procedure (6) of calculation of ρ ± n resembles the technique of iterative approximation of the joint spectral radius suggested in [14].
Main Theorem. For any irreducible matrix set A , nonzero vector e ∈ R m , initial norm · 0 , and any averaging function γ(t, s), the sequences {ρ ± n } constructed by the iteration procedure MR1, MR2 converge to ρ(A ), and the sequence of norms · n converges uniformly on each bounded set to some Barabanov norm · * of the matrix set A . Moreover, the sequence {ρ − n } is nondecreasing, the sequence {ρ + n } is nonincreasing, and ρ − n ≤ ρ(A ) ≤ ρ + n for all n = 1, 2, . . . , which provides an a posteriori estimate for the computational error of ρ(A ).

Proof of Main Theorem.
Let us suppose that we managed to prove the following assertions: A1: the sequences {ρ + n } and {ρ − n } are convergent; A2: the limits of the sequences {ρ + n } and {ρ − n } coincide: A3: the norms · • n converge pointwise to a limit · * . Then the function · * will be a semi-norm in R m . Moreover, by (8) each norm · • n meets the normalization condition e • n = 1, and hence e * = lim which implies x * ≡ 0. Note also that due to (8) the norms · • n differ from · n only by numerical factors. Therefore, the quantities ρ ± n can be defined as Then, passing to the limit in (9), one can conclude that the semi-norm x * satisfies the Barabanov condition But as shown in [21,Thm. 3], any semi-norm x * ≡ 0 satisfying the Barabanov condition for an irreducible matrix set is a Barabanov norm. Thus, under assumptions A1, A2 and A3, the iteration procedure (6)-(8) allows to build a Barabanov norm and to find the joint spectral radius of the matrix set A .
So, to complete the proof of Main Theorem we need to justify assertions A1, A2 and A3 which will be done in Sections 3.1-3.6. In Section 3.1 we establish convergence of the sequence of norms { · • n } to some limit which allows to prove in Lemma 3.2 that Assertion A3 is a corollary of Assertions A1 and A2. Section 3.2 demonstrates that the quantities {ρ − n } form the family of lower bounds for the joint spectral radius ρ of the matrix set A , while the quantities {ρ + n } form the family of upper bounds for ρ. In Section 3.3 we prove that the sequences {ρ ± n } are bounded and monotone which implies the existence of the limits ρ − = lim n→∞ ρ − n and ρ + = lim n→∞ ρ + n . At last, in Sections 3.4 and 3.5 we prove that ρ − = ρ + which allows to justify in Section 3.6 the validity of Assertions A1 and A2 and thus to finalize the proof of Main Theorem.

3.1.
Convergence of the sequence of norms { · • n }. Given a pair of norms · ′ and · ′′ in R m define the quantities Since all norms in R m are equivalent to each other, the quantities e − ( · ′ , · ′′ ) and e + ( · ′ , · ′′ ) are correctly defined and Therefore the quantity which is called the eccentricity of the norm · ′ with respect to the norm · ′′ (see, e.g., [36]), is also correctly defined.
Lemma 3.1. Let · * be a Barabanov norm for the matrix set A . Then and the sequence of the numbers ecc( · n , · * ) is nonincreasing.
Proof. Note first that by (8) each norm · • n differs from the corresponding norm · n only by a numerical factor. From this, by the definition (10), (11) of the eccentricity of one norm with respect to another, the equality (12) follows.
Denote by ρ the joint spectral radius of the matrix set A . Then, by definitions of the function e + (·) and of the Barabanov norm · * , from (6), (7) we obtain: Similarly, by definitions of the function e − (·) and of the Barabanov norm · * , from (6), (7) we obtain: By dividing termwise the inequality (13) on (14) we get Hence, the sequence {ecc( · n , · * )} is nonincreasing. Denote by N loc (R m ) the topological space of norms in R m with the topology of uniform convergence on bounded subsets of R m .
Therefore the norms · n , n ≥ 1, are equicontinuous and uniformly bounded on each bounded subset of R m . Moreover, their values are also uniformly separated from zero on each bounded subset of R m separated from zero. From here by the Arzela-Ascoli theorem the statement of the corollary follows.
Corollary 2. If at least one of subsequences of norms from { · • n } converges in N loc (R m ) to some Barabanov norm then the whole sequence { · • n } also converges in N loc (R m ) to the same Barabanov norm.
Proof. Let { · • n k } be a subsequence of { · • n } which converges in N loc (R m ) to some Barabanov norm · * . Then by definition of the eccentricity of one norm with respect to another ecc( · • n k , · * ) → 1 as k → ∞. Here by Lemma 3.1 the eccentricities ecc( · • n , · * ) are nonincreasing in n, and then the following stronger relation holds Note now that by the definition (10), (11) of the eccentricity of one norm with respect to another from which by (15) it follows that the sequence of norms { · • n } converges in space N loc (R m ) to the norm · * . Lemma 3.2. Assertion A3 is a corollary of Assertions A1 and A2.
Proof. By Corollary 1 the sequence of norms { · • n } has a subsequence { · • n k } that converges in space N loc (R m ) to some norm · * . Then, passing to the limit in (9) as n = n k → ∞, we get by Assertions A1 and A2: which means that · * is a Barabanov norm for the matrix set A . This and Corollary 2 then imply that the sequence { · • n } converges in space N loc (R m ) to the Barabanov norm · * . Assertion A3 is proved.
So, in view of Lemma 3.2 to prove that the iteration procedure (6)-(8) is convergent it suffices to verify only that Assertions A1 and A2 hold.

3.2.
Relations between ρ ± n and ρ. The following lemma provides a way to estimate the spectral radius of a matrix set. Lemma 3.3. Let α, β be numbers such that in some norm · the inequalities hold. Then α ≤ ρ ≤ β, where ρ is the joint spectral radius of the matrix set A .
Proof. Let · * be some Barabanov norm for the matrix set A . Since all norms in R m are equivalent, there are constants σ − > 0 and σ + < ∞ such that Consider, for each k = 1, 2, . . ., the functions Then, as is easy to see, Similarly consider, for each k = 1, 2, . . ., the functions For these functions, by definition of Barabanov norms the following identity hold which is stronger than (17). Now, note that (16) and the definition of the functions ∆ k (x) and ∆ So, Lemma 3.3 and the definition (6) of ρ ± n imply that the quantities {ρ − n } form the family of lower bounds for the joint spectral radius ρ of the matrix set A , while the quantities {ρ + n } form the family of upper bounds for ρ. This allows to estimate a posteriori errors of computation of the joint spectral radius with the help of the iteration procedure (6)-(8).

Convergence of the sequences {ρ
Here by the definition (6) of the quantities ρ ± n the right-hand part of the chain of equalities can be estimated as follows: Therefore, by definition of the norm x n+1 , and then, ρ − n ≤ ρ − n+1 ≤ ρ + n+1 ≤ ρ + n . So, the following lemma holds.
To prove that ρ − = ρ + , below it will be supposed the contrary, which will lead us to a contradiction.

Transition to a new sequence of norms.
To simplify further constructions we will switch over to a new sequence of norms for which the quantities ρ ± n will be independent of n.
By Corollary 1 the sequence of the norms · • n is compact in space N loc (R m ). Hence, there is a subsequence of indices {n k } such that the the norms · • n k = · n k / e n k converge to some norm · • 0 satisfying the normalization condition e • 0 = 1. Then, passing to the limit in (9), by Lemma 3.4 we obtain: Now by induction the following statement can be easily proved.
Proof. The statement of the lemma is obvious for x = 0. So, suppose that x ∈ ω n , x = 0. In this case (21) and the inequalities ρ − ≤ ρ + imply From here by the definition (20) of the norm · • n+1 we get the required equality: The lemma is proved.
Proof. Let x ∈ ω n+1 . If x = 0 then clearly x ∈ ω n , so suppose that x = 0. By definitions of the set ω n+1 and of the norm · • n the following equalities take place: Then from (22) it follows that max max But by the conditions of the lemma ρ − < ρ + . Then γ = γ(ρ − , ρ + ) > ρ − , and the right-hand part of the above equalities is strictly less than max i A i x • n . A contradiction, since the left-hand part of the same equalities is no less than max i A i x • n . The above contradiction is caused by the assumption (23), and therefore it is proved that the condition x = 0 ∈ ω n+1 implies the strict inequality In this case from (22) it follows that Let us show that the equality (24) implies Indeed, supposing the contrary, by definition of the quantity ρ − , there should be valid the strict inequality max i A i x • n > ρ − x • n . Then the left-hand part of the equality (24) should be strictly greater than ρ − x • n , that is greater than the righthand part of the same equality, which is impossible. This last contradiction shows that the equality (25) holds as soon as x = 0 ∈ ω n+1 , which means by (19) that x ∈ ω n . Corollary 3. If ρ − < ρ + then ω := ∩ n≥0 ω n = 0 and Proof. By Lemma 3.7 {ω n } is the family of embedded closed conic 3 sets. Then the intersection ω of these sets is also a closed conic set such that ω = {0}. By definition of the set ω, if x ∈ ω then x ∈ ω n for every integer n ≥ 0. Hence, by Lemma 3.6 x • n+1 = x • n , from which the equalities (26) follow.
3.6. Completion of the proof of Assertion A2. By Corollary 3 there is a nonzero vector g on which all the norms · • n take the same values: Then, due to uniform boundedness of the eccentricities of the norms · • n with respect to some Barabanov norm · * (this fact can be proved by verbatim repetition of the analogous proof for the norms · n ), the norms · • n form a family which is uniformly bounded and equicontinuous with respect to the Barabanov norm · * :

Hence, by the Arzela-Ascoli theorem the family of norms
From the definition (20) of the norms · • n it follows also that Then the norms · • n are monotone increasing in n and bounded (with respect to the Barabanov norm · * ) and therefore they pointwise converge to some norm · • . Moreover, since the family of norms { · • n } is equicontinuous with respect to the Barabanov norm · * , the norms · • n converge to the norm · • in space N loc (R m ). Now, passing to the limit in the relations which follow from (20), we obtain On the other hand, passing to the limit in the first relation of (19), we obtain Relations (27) and (28) imply the inequality ρ + ≤ γ which contradicts the assumption ρ − < ρ + because by definition of the function γ(·, ·) the condition ρ − < ρ + implies the inequality γ = γ(ρ − , ρ + ) < ρ + .
The obtained contradiction completes the proof of the equality ρ − = ρ + as well as of the convergence of the iteration procedure (6)-(8).
Then the norm x of the vector x with polar coordinates (r, ϕ) is determined by the equality and the unit sphere in the norm · is determined as the geometrical locus of the vectors x polar coordinates of which satisfy the relations see Fig. 1.
Now, let R n (ϕ) be the function defining in the polar coordinates the graph of the unit sphere x n = 1 of the norm · n determined by the iteration procedure (6)- (8). Rewrite the relations (6)- (8) in terms of the functions R n (ϕ). To do it we should express the quantities A i x n , i = 0, 2, . . . , r, in terms of the functions R n (ϕ).
By (29) A i x n = r(A i x)R n (ϕ(A i x)).

Example 1. Consider the family
The functions Φ i (ϕ), H i (ϕ), R n (ϕ), R * n (ϕ) were chosen to be piecewise linear with 3000 nodes uniformly distributed over the interval [−π, π]. It was needed 13 iterations of algorithm MR1-MR2 with the averaging function γ(t, s) = t+s 2 implemented in MATLAB to compute the joint spectral radius ρ(A ) with the absolute accuracy 10 −3 . The computed value of the joint spectral radius is ρ(A ) = 1.617. The computed unit sphere of the Barabanov norm · * after the 13th iteration of algorithm MR1-MR2 is shown on Fig. 2 on the left.
As is seen from Fig. 2, in Example 1 the sets A 1 x = ρ and A 2 x = ρ have exactly 4 intersection points. This was theoretically proved in [18,21] for the case when one of the matrices A 1 , A 2 is lower triangle and the other is upper triangle, and their entries are nonnegative. In [18,21] this fact was one of key points in disproving the Finiteness Conjecture. We do not know whether this fact is valid in a general case or not, but numerical tests based on algorithm MR1-MR2 with several dozens pairs of matrices A 1 , A 2 testify for this fact. Here the functions Φ i (ϕ), H i (ϕ), R n (ϕ), R * n (ϕ) were also chosen to be piecewise linear with 3000 nodes uniformly distributed over the interval [−π, π]. It was needed 31 iterations of algorithm MR1-MR2 with the averaging function γ(t, s) = t+s 2 implemented in MATLAB to compute the joint spectral radius ρ(A ) with the absolute accuracy 10 −3 . The computed value of the joint spectral radius is ρ(A ) = 1.347. The computed unit sphere of the Barabanov norm · * after the 31d iteration of algorithm MR1-MR2 is shown on Fig. 2 on the right. The MATLAB code for this Example can be found in [20]. 5. Concluding remarks. The max-relaxation algorithm suggested in the paper allows to calculate the joint spectral radius of a finite matrix family (of arbitrary matrix size and arbitrary amount of matrices in the set) with any required accuracy and to evaluate a posteriori the computational error. Still, this algorithm gives rise to a set of open problems.
Problem. While the quantities {ρ ± n } provide a posteriori bounds for the accuracy of approximation of ρ(A ) the question about the accuracy of approximation of the Barabanov norm · * by the norms · • n is open. It seems, the difficulty in resolving this problem is caused by the fact that in general the Barabanov norms for a matrix family are determined ambiguously. Namely to overcome this difficulty we preferred to consider relaxation algorithm instead of direct one. Moreover, as can be shown, both theoretically and numerically, if to set x n+1 = γ −1 n max i A i x n in (7) then the obtained direct computational analog of algorithm MR1-MR2 may turn out to be non-convergent.
Problem. The question about the rate of convergence of the sequences {ρ + n } and {ρ − n } to the joint spectral radius is also open.
Remark also that in this paper mainly the algorithm for building of Barabanov norms rather than its computational details was studied. The numerical aspects of implementation of this algorithm require additional analysis.
Problem. An estimation of the computational cost of the max-relaxation algorithm is required. The dependance of the algorithm on parameters r, m and the choice of the averaging function γ(t, s) is acute, too.