arXiv:1102.4180v1  [cs.RO]  21 Feb 2011
apport 

de recherche
ISSN 0249-6399
ISRN INRIA/RR--7544--FR+ENG
Robotics
INSTITUT NATIONAL DE RECHERCHE EN INFORMATIQUE ET EN AUTOMATIQUE
Characterizing and approximating eigenvalue sets of
symmetric interval matrices
Milan Hladík — David Daney — Elias Tsigaridas
N° 7544
Février 2011
Centre de recherche INRIA Sophia Antipolis – Méditerranée
2004, route des Lucioles, BP 93, 06902 Sophia Antipolis Cedex
Téléphone : +33 4 92 38 77 77 — Télécopie : +33 4 92 38 77 65
Characterizing and approximating eigenvalue sets
of symmetric interval matrices
Milan Hladík∗† , David Daney† , Elias Tsigaridas‡
Theme : Robotics
Perception, Cognition, Interaction
Équipe-Projet Coprin
Rapport de recherche n° 7544 — Février 2011 — 19 pages
Abstract:
We consider the eigenvalue problem for the case where the input
matrix is symmetric and its entries perturb in some given intervals. We present
a characterization of some of the exact boundary points, which allows us to
introduce an inner approximation algorithm, that in many case estimates exact
bounds. To our knowledge, this is the first algorithm that is able to guaran-
tee exactness. We illustrate our approach by several examples and numerical
experiments.
Key-words:
Interval matrix, symmetric matrix, interval analysis, eigenvalue,
eigenvalue bounds
∗
Charles
University,
Faculty
of
Mathematics
and
Physics,
Department
of
Ap-
plied
Mathematics,
Malostranské
nám.
25,
11800,
Prague,
Czech
Republic,
e-mail:
milan.hladik@matfyz.cz
† INRIA Sophia-Antipolis Méditerranée, 2004 route des Lucioles, BP 93, 06902 Sophia-
Antipolis Cedex, France, e-mail: FirstName.LastName(AT)inria.fr
‡ Computer Science Department, Aarhus University, Denmark. e-mail: elias@cs.au.dk
Caractérisation et approximation des ensembles
de valeurs propres des matrices symétriques
intervalles
Résumé : We consider the eigenvalue problem for the case where the input
matrix is symmetric and its entries perturb in some given intervals. We present
a characterization of some of the exact boundary points, which allows us to
introduce an inner approximation algorithm, that in many case estimates exact
bounds. To our knowledge, this is the first algorithm that is able to guaran-
tee exactness. We illustrate our approach by several examples and numerical
experiments.
Mots-clés :
Matrice d’intervalles, matrice symétrique, analyse par intervalles,
valeur propre, borne sur les valeurs propres
Characterizing and approximating eigenvalue sets of symmetric interval matrices3
1
Introduction
Computing eigenvalues of a matrix is a basic linear algebraic task used through-
out in mathematics, physics and computer science. Real life makes this prob-
lem more complicated by imposing uncertainties and measurement errors on the
matrix entries. We suppose we are given some compact intervals in which the
matrix entries can perturb. The set of all possible real eigenvalues forms a com-
pact set, and the question that we deal with in this paper is how to characterize
and compute it.
The interval eigenvalue problem has its own history. The first results are
probably due to Deif [10] and Rohn & Deif [34]: bounds for real and imaginary
parts for complex eigenvalues were studied by Deif [10], while Rohn & Deif
[34] considered real eigenvalues. Their theorems are applicable only under an
assumption on sign pattern invariancy of eigenvectors, which is not easy to
verify (cf. [8]). A boundary point characterization of eigenvalue set was given
by Rohn [30], and it was used by Hladík et al. [17] to develop a branch & prune
algorithm producing an arbitrarily tight approximation of the eigenvalue set.
Another approximate method was given by Qiu et al. [28]. The related topic of
finding verified intervals of eigenvalues for real matrices was studied in e.g. [4].
In this paper we focus on the case of the symmetric eigenvalue problem.
Symmetric matrices naturally appear in many practical problems, but symmet-
ric interval matrices are hard to deal with. This is so, mainly due to the so-called
dependencies, that is, correlations between the matrix components. If we “for-
get” these dependencies and solve the problem by reducing it to the previous
case, then the results will be greatly overestimated, in general (but not the ex-
tremal points, see Theorem 2). From now on we consider only the symmetric
case.
Due to the dependencies just mentioned, the theoretical background for the
eigenvalue problem of symmetric interval matrices is not wide enough and there
are few practical methods. The known results are that by Deif [10] and Hertz
[15]. Deif [10] gives an exact description of the eigenvalue set together with re-
strictive its assumptions. Hertz [15] (cf. [32]) proposed a formula for computing
two extremal points of the eigenvalue set—the largest and the smallest one. As
the problem itself is very hard, it is not surprising conjectures on the problem
[29] turned out to be wrong [35].
In the recent years, several approximation algorithms have been developed.
An evolution strategy method by Yuan et al. [35] yields an inner approximation
of the eigenvalue set. By means of matrix perturbation theory, Qiu et al. [27]
proposed an algorithm for approximate bounds, and Leng & He [23] for an
outer estimation. An outer estimation was also considered by Kolev [22], but
for general case with nonlinear dependencies. Some initial bounds that are easy
and fast to compute were discussed by Hladík et al. [18]. An iterative algorithm
for outer estimation was given by Beaumont [6].
This problem has a lot of applications in the field of mechanics and engi-
neering. Let us mention for instance automobile suspension systems [28], mass
structures [27], vibrating systems [11], principal component analysis [12], and
robotics [7]. Another applications arise from the engineering area concerning
singular values and condition numbers. Using the well-known Jordan–Wielandt
transformation [13, 20, 25] we can simply reduce a singular value calculation to
a symmetric eigenvalue one.
RR n° 7544
Characterizing and approximating eigenvalue sets of symmetric interval matrices4
The rest of the paper is structured as follows: In Sec. 2 we introduce the
notation that we use throughout the paper.
In Sec. 3 we present our main
theoretical result, and in Sec. 4 we develop our algorithms for the problem.
Finally, in Sec. 5 we demonstrate our approach by a number of examples and
numerical experiments, and we conclude in Sec. 6.
2
Basic definitions and theoretical background
Let us introduce some notions first. An interval matrix is defined as
A := [A, A] = {A ∈Rm×n; A ≤A ≤A},
where A, A ∈Rm×n, A ≤A, are given matrices. By
Ac := 1
2(A + A),
A∆:= 1
2(A −A)
we denote the midpoint and the radius of A, respectively.
By an interval linear system of equations Ax = b we mean a family of
systems Ax = b, such that A ∈A, b ∈b. In a similar way we introduce interval
linear systems of inequalities and mixed systems of equations and inequalities.
A vector x is a solution of Ax = b if it is a solution of Ax = b for some A ∈A
and b ∈b. We assume that the reader is familiar with the basics of interval
arithmetic; for further details we refer to e.g. [3, 14, 26].
Let F be a family of n × n matrices. We denote the eigenvalue set of the
family F by
Λ(F) := {λ ∈R; Ax = λx, x ̸= 0, A ∈F}.
A symmetric interval matrix as defined as
AS := {A ∈A | A = AT }.
It is usually a proper subset of A. Considering the eigenvalue set Λ(A), it,
generally, represents an overestimation of Λ(AS). That is why we focus directly
on eigenvalue set of the symmetric counterpart, even though we must take into
account the dependencies between the elements, in the definition of AS.
A real symmetric matrix A ∈Rn×n has always n real eigenvalues, let us sort
them in a non-increasing order
λ1(A) ≥λ2(A) ≥· · · ≥λn(A).
We extend this notation for symmetric interval matrices
λi(AS) :=

λi(A) | A ∈AS	
.
These sets represent n compact intervals λi(AS) = [λi(AS), λi(AS)], i =
1, . . . , n; cf. [18]. The intervals can be disjoint, can overlap, or some of them,
can be identical. However, what it can not happen is one interval to be a proper
subset of another interval. The union of these intervals produces Λ(AS).
Throughout the paper we use the following notation:
RR n° 7544
Characterizing and approximating eigenvalue sets of symmetric interval matrices5
λi(A)
the ith eigenvalue of a symmetric matrix A (in a non-
increasing order)
σi(A)
the ith singular value of a matrix A (in a non-increasing
order)
vi(A)
the eigenvector associated to the ith eigenvalue of a sym-
metric matrix A
ρ(A)
the spectral radius of a matrix A
∂S
the boundary of a set S
conv S
the convex hull of a set S
diag(y)
the diagonal matrix with entries y1, . . . , yn
sgn(x)
the
sign
vector
of
a
vector
x,
i.e.,
sgn(x)
=
(sgn(x1), . . . , sgn(xn))T
∥x∥2
the Euclidean vector norm, i.e., ∥x∥2 =
√
xT x
∥x∥∞
the Chebyshev (maximum) vector norm, i.e., ∥x∥∞
=
max {|x|i; i = 1, . . . , n}
x ≤y, A ≤B
vector and matrix relations are understood component-wise
3
Main theorem
The following theorem is the main theoretical result of the the present paper.
We remind the reader that the principal m × m submatrix of a given n × n
matrix is any submatrix obtained by eliminating any n −m rows and the same
n −m columns.
Theorem 1. Let λ ∈∂Λ(AS). There is a principal submatrix e
AS ∈Rk×k of
AS such that:
• If λ = λj(AS) for some j ∈{1, . . . , n}, then
λ ∈

λi
 eAc + diag(z) eA∆diag(z)


; z ∈{±1}k, i = 1, . . . , k
	
.
(1)
• If λ = λj(AS) for some j ∈{1, . . . , n}, then
λ ∈

λi
 eAc −diag(z) eA∆diag(z)


; z ∈{±1}k, i = 1, . . . , k
	
.
(2)
Proof. Let λ ∈∂Λ(AS). Then either λ = λj(AS) or λ = λj(AS), for some
j ∈{1, . . . , n}. We assume the former case. The latter could be proved similarly.
The eigenvalue λ is achieved for a matrix A ∈A.
It is without loss of
generality to assume that the corresponding eigenvector x, ∥x∥2 = 1, is of the
form x = (0T , yT )T , where y ∈Rk and yi ̸= 0, for all 1 ≤i ≤k, and for some
k ∈{1, . . . , n}. The symmetric interval matrix AS can be written as
AS =
BS
C
CT
DS

,
where BS ⊂R(n−k)×(n−k), C ⊂R(n−k)×k, DS ⊂Rk×k.
From the basic equality Ax = λx it follows that
Cy = 0 for some C ∈C,
(3)
RR n° 7544
Characterizing and approximating eigenvalue sets of symmetric interval matrices6
and
Dy = λy for some D ∈DS.
(4)
We focus on the latter relation; it says that λ is an eigenvalue of D. We will
show that DS is the required principal submatrix e
AS and D could be written
as in (1).
From (4) we have that λ = yT Dy, and hence the partial derivatives are
∂λ
∂dij
= yiyj ̸= 0,
i, j = 1, . . . , k.
This relation strongly influences the structure of D. If yiyj > 0, then dij = dij.
This is so, because otherwise by increasing dij we also increase the value of λ,
which contradicts our assumption that λ lies on the upper boundary of Λ(AS).
Likewise, yiyj < 0 implies dij = dij. This allows us to write D in the following
more compact form
D = Dc + diag(z) D∆diag(z),
(5)
where z = sgn(y) ∈{±1}k. Therefore, λ belongs to a set as the one presented
in the right-hand side of (1), which completes the proof.
Note that not every λj(AS) or λj(AS) is a boundary point of Λ(AS). The-
orem 1 is also true for such λj(AS) or λj(AS) that are non-boundary, but
make no multiple eigenvalue (since the corresponding eigenvector is uniquely
determined). However, truthfulness of Theorem 1 for all λj(AS) and λj(AS),
j = 1, . . . , n, is still an open question. Moreover, full characterization of all
λj(AS) and λj(AS), j = 1, . . . , n, is lacking, too.
As we have already mentioned, in general, the eigenvalue set of an interval
matrix is larger than the eigenvalue set of its symmetric counterpart. This is
true even if both the midpoint and radius matrices are symmetric (see Exam-
ple 1). The following theorem says that overestimation caused by the additional
matrices is somehow limited by the intermediate area.
Theorem 2. Let Ac, A∆∈Rn×n be symmetric matrices. Then
conv Λ(AS) = conv Λ(A).
Proof. The inclusion conv Λ(AS) ⊆conv Λ(A) follows from the definition of the
convex hull.
Let A ∈A be arbitrary, λ one of its real eigenvalues, and x the corresponding
eigenvector, where ∥x∥2 = 1. Let B := 1
2(A + AT ) ∈AS, then the following
holds:
λ = xT Ax ≤max
∥y∥2=1 yT Ay = max
∥y∥2=1 yT By = λ1(B) ∈conv Λ(AS).
Similarly,
λ = xT Ax ≥min
∥y∥2=1 yT Ay = min
∥y∥2=1 yT By = λn(B) ∈conv Λ(AS).
Therefore λ ∈conv Λ(AS), and so conv Λ(A) ⊆conv Λ(AS), which completes
the proof.
RR n° 7544
Characterizing and approximating eigenvalue sets of symmetric interval matrices7
4
Inner approximation algorithms
Theorem 1 naturally yields an algorithm to compute a very sharp inner approx-
imation of Λ(AS), which could also be exact in some cases. We will present the
algorithm in the sequel (Section 4.3). First, we define some notions and propose
two simple but useful methods for inner approximations.
Any subset of S is called an inner approximation. Similarly, any set that con-
tains S is called an outer approximation. In our case, an inner approximation of
the eigenvalue set λi(AS), is denoted by µi(AS) = [µi(AS), µi(AS)] ⊆λi(AS),
and an outer approximation is denoted by ωi(AS) = [ωi(AS), ωi(AS)] ⊇λi(AS),
where 1 ≤i ≤n.
From a practical point of view, an outer approximation is usually more use-
ful. However, an inner approximation is also important in some applications.
For example, it could be used to measure quality (sharpness) of an outer approx-
imation, or it could be used to prove (Hurwitz or Schur) unstability of certain
interval matrices, cf. [31].
4.1
Local improvement
The first algorithm that we present is based on local improvement search tech-
nique. A similar method, but for another class of symmetric interval matrices
was proposed by Rohn [31]. The basic idea of the algorithm is to start with an
eigenvalue, λi(Ac), and the corresponding eigenvector, vi(Ac), of the midpoint
matrix, Ac, and then move to an extremal matrix in AS according to the sign
pattern of the eigenvector. The procedure is repeated until no improvement is
possible.
Algorithm 1 outputs the upper boundaries µi(AS) of the inner approxima-
tion [µi(AS), µi(AS)], where 1 ≤i ≤n. The lower boundaries, µi(AS), can
be obtained similarly. The validity of the procedure follows from the fact that
every considered matrix, A, belongs to AS.
Algorithm 1 (Local improvement for µi(AS), i = 1, . . . , n)
1: for i = 1, . . . , n do
2:
µi(AS) = −∞;
3:
A := Ac;
4:
while λi(A) > µi(AS) do
5:
D := diag(sgn(vi(A)));
6:
A := Ac + DA∆D;
7:
µi(AS) := λi(A);
8:
end while
9: end for
10: return µi(AS), i = 1, . . . , n.
The algorithm terminates after, at most, 2n iterations. However, usually in
practice the number of iterations is much smaller, which makes the algorithm
attractive for applications. Our numerical experiments (Section 5) indicate that
the number of iterations is rarely greater than two, even for matrices of dimen-
sion 20. Moreover, the resulting inner approximation is quite sharp, depending
on the width of intervals in AS. This is not surprising as whenever the input
RR n° 7544
Characterizing and approximating eigenvalue sets of symmetric interval matrices8
intervals are narrow enough, the algorithm produces, sometimes even after the
first iteration, exact bounds; see [10]. We refer the reader to Section 5 for a
more detailed presentation of the experiments.
4.2
Vertex enumeration
The second method that we present is based on vertex enumeration. It consists
of inspecting all matrices
Az := Ac + diag(z) A∆diag(z),
z ∈{±1}n,
z1 = 1,
and continuously improving an inner approximation µi(AS), whenever λi(Az) >
µi(AS), where 1 ≤i ≤n. The lower bounds, µi(AS), could be obtained in a
similar way using the matrices Ac −diag(z) A∆diag(z), where z ∈{±1}n, and
z1 = 1. The condition z1 = 1 follows from the fact that diag(z) A∆diag(z) =
diag(−z) A∆diag(−z), which gives us the freedom to fix one component of
z. The number of steps that the algorithm performs is 2n−1. Therefore, this
method is suitable only for matrices of moderate dimensions.
The main advantages of the vertex enumeration approach are the following.
First, it provides us with sharper inner approximation of the eigenvalue sets
than the local improvement. Second, two of the computed bounds are exact;
by Hertz [15] (cf. [32]) and Hertz [16] we have that µ1(AS) = λ1(AS) and
µn(AS) = λn(AS). The other bounds calculated by the vertex enumeration,
even though it was conjectured that there were exact [29], it turned out that
they were not exact, in general [35] The assertion by Hertz [16, Theorem 1]
that µ1(AS) = λ1(AS) and µn(AS) = λn(AS) is wrong, too; see Example 3.
Nevertheless, Theorem 1 and its proof indicate a sufficient condition: If no
eigenvector corresponding to an eigenvalue of AS has a zero component, then
the vertex enumeration yields exactly the eigenvalue sets λi(AS), i = 1, . . . , n.
The efficient implementation of this approach is quite challenging. In or-
der to overcome in practice the exponential complexity of the algorithm, we
implemented a branch & bound algorithm, which is in accordance with the sug-
gestions of Rohn [31]. However, the adopted bounds are not that tight, and the
actual running times are usually worse than the direct vertex enumeration. That
is why we do not consider further this variant. The direct vertex enumeration
scheme for computing the upper bounds, µi(AS), is presented in Algorithm 2.
4.3
Submatrix vertex enumeration
In this section we present an algorithm that is based on Theorem 1, and it usually
produces very tight inner approximations, even exact ones in some cases. The
basic idea the algorithm is to enumerate all the vertices of all the principal
submatrices of AS. The number of steps performed with this approach is
2n−1 + n2n−2 +
n
2

2n−2 + · · · + n20 = 1
2 (3n −1) .
To overcome the obstacle of the exponential number of iterations, at least in
practice, we notice that not all eigenvalues of the principal submatrices of the
matrices in AS belong to some of the eigenvalue sets λi(AS), where 1 ≤i ≤n.
For this we will introduce a condition for checking such an inclusion.
RR n° 7544
Characterizing and approximating eigenvalue sets of symmetric interval matrices9
Algorithm 2 (Vertex enumeration for µi(AS), i = 1, . . . , n)
1: for i = 1, . . . , n do
2:
µi(AS) = λi(Ac);
3: end for
4: for all z ∈{±1}n, z1 = 1, do
5:
A := Ac + diag(z) A∆diag(z);
6:
for i = 1, . . . , n do
7:
if λi(A) > µi(AS) then
8:
µi(AS) := λi(A);
9:
end if
10:
end for
11: end for
12: return µi(AS), i = 1, . . . , n.
Assume that we are given an inner approximation µi(AS) and an outer ap-
proximation ωi(AS) of the eigenvalue sets λi(AS); that is µi(AS) ⊆λi(AS) ⊆
ωi(AS), where 1 ≤i ≤n. As we will see in the sequel, the quality of the output
of our methods depends naturally on the sharpness of the outer approximation
used.
Let DS ⊂Rk×k be a principal submatrix of AS and, without loss of gener-
ality, assume that it is situated in the right-bottom corner, i.e.,
AS =
BS
C
CT
DS

,
where BS ⊂R(n−k)×(n−k) and C ⊂R(n−k)×k. Let λ be an eigenvalue of some
vertex matrix D ∈DS which is of the form (5), and let y be the corresponding
eigenvector. If the eigenvector is not unique then λ is a multiple eigenvalue and
therefore it is a simple eigenvalue of some principal submatrix of DS; in this
case we restrict our consideration to this submatrix.
We want to determine whether λ is equal to λp(AS) ∈Λ(AS) for some fixed
p ∈{1, . . ., n}, or if this is not possible, to improve the upper bound µp(AS);
the lower bound can be handled accordingly. In view of (3), it must be
0 ∈Cy ,
so that λ to be an eigenvalue of some matrix in AS. Now, we are sure that
λ ∈Λ(AS) and it remains to determine whether λ also belongs to λp(AS).
If λ ≤µp(AS), then it is useless to further considering λ, as it will not extend
the inner approximation of the pth eigenvalue set. If p = 1 or λ < ωp−1(AS),
then λ must belong to λp(AS), and we can improve the inner bound µp(AS) :=
λ. In this case the algorithm terminates early, and that is the reason we need
ωi(AS) to be as tight as possible, where 1 ≤i ≤n.
If p > 1 and λ ≥ωp−1(AS), we proceed as follows. We pick an arbitrary
C ∈C, such that Cy = 0; we refer to, e.g. [33] for details on the selection
process. Next, we select an arbitrary B ∈BS and let
A :=

B
C
CT
D

.
(6)
RR n° 7544
Characterizing and approximating eigenvalue sets of symmetric interval matrices10
We compute the eigenvalues of A, and if µp(AS) < λp(A), then we set µp(AS) :=
λp(A), otherwise we do nothing.
However, it can happen that λ = λi(AS), and we do not identify it, and
hence we do not enlarge the inner estimation µp(AS). Nevertheless, if we apply
the method for all p = 1, . . . , n and all principal submatrices of AS, then we
touch all the boundary points of Λ(AS). If λ ∈∂Λ(AS), then λ is covered by
the resulting inner approximation. In the case when λ is an upper boundary
point, we consider the maximal i ∈{1, . . . , n} such that λ = λi(AS) and then
the ith eigenvalue of the matrix (6) must be equal to λ. Similarly test are valid
for a lower boundary point.
Now we have all the ingredients at hand for the direct version of the sub-
matrix vertex enumeration approach that is presented in Algorithm 3, which
improves the upper bound µp(AS) of an inner approximation, where the index
p is still fixed. Let us also mention that in step 4 of Algorithm 3, the decom-
position of AS according to the index set J means that DS is a restriction of
AS to the rows and the columns indexed by J, BS is a restriction of AS to the
rows and the columns indexed by {1, . . . , n} \ J, and C is a restriction of AS
to the rows indexed by {1, . . . , n} \ J and the columns indexed by J.
Algorithm 3 (Direct submatrix vertex enumeration for µp(AS))
1: compute outer approximation ωi(AS), i = 1, . . . , n;
2: call Algorithm 1 to get inner approximation µi(AS), i = 1, . . . , n;
3: for all J ⊆{1, . . . , n}, J ̸= ∅, do
4:
decompose AS =

BS
C
CT DS

according to J;
5:
for all z ∈{±1}|J|, z1 = 1, do
6:
D := Dc + diag(z) D∆diag(z);
7:
for i = 1, . . . , |J| do
8:
λ := λi(D);
9:
y := vi(D);
10:
if λ > µp(AS) and λ ≤ωp(AS) and 0 ∈Cy then
11:
if p = 1 or λ < ωp−1(AS) then
12:
µp(AS) := λ;
13:
else
14:
find C ∈C such that Cy = 0;
15:
A :=

Bc C
CT D

;
16:
if λp(A) > µp(AS) then
17:
µp(AS) := λp(A);
18:
end if
19:
end if
20:
end if
21:
end for
22:
end for
23: end for
24: return µp(AS).
RR n° 7544
Characterizing and approximating eigenvalue sets of symmetric interval matrices11
4.3.1
Branch & bound improvement
In order to tackle the exponential worst case complexity of Algorithm 3, we
propose the following modification. Instead of inspecting all non-empty subsets
of {1, . . . , n} in step 3, we exploit a branch & bound method, which may skip
some useless subsets.
Let a non-empty J ⊆{1, . . . , n} be given.
The new,
possible improved, eigenvalue λ must lie in the interval λ := [µp(AS), ωp(AS)].
If this is the case, then the interval matrix AS −λI must be irregular, i.e., it
contains a singular matrix. Moreover, the interval system
(AS −λI)x = 0,
∥x∥∞= 1 ,
has a solution x, where xi = 0 for all i ̸∈J. We decompose AS −λI according to
J, and, without loss of generality, we may assume that J = {n −|J| + 1, . . . , n},
then
AS −λI =
BS −λI
C
CT
DS −λI

.
The interval system becomes
Cy = 0, (DS −λI)y = 0, ∥y∥∞= 1,
(7)
where we considered x = (0T , yT )T . This is a very useful necessary condition.
If (7) has no solution, then we cannot improve the current inner approximation.
We can also prune the whole branch with J as a root; that is, we will inspect no
index sets J′ ⊆J. The strength of this condition follows from the fact that the
system (7) is overconstrained, it has more equations than variables. Therefore,
with high probability, that it has no solution, even for larger J.
Let us make two comments about the interval system (7). First, this system
has a lot of dependencies. They are caused from the multiple occurrences of λ,
and by the symmetry of DS. If no solver for interval systems that can handle
dependencies is available, then we can solve (7) as an ordinary interval system,
“forgetting” the dependencies. The necessary condition will be weaker, but still
valid. This is what we did in our implementation.
The second comment addresses the expression ∥y∥∞= 1. We have chosen
the maximum norm to pertain linearity of the interval system. The expression
could be rewritten as −1 ≤y ≤1 (for checking solvability of (7) we can use
either normalization ∥y∥∞= 1 or ∥y∥∞≤1). Another possibility is to write
−1 ≤y ≤1, yi = 1 for some i ∈{1, . . . , |J|}.
This indicates that we can split the problem into solving |J| interval systems
Cy = 0, (DS −λI)y = 0, −1 ≤y ≤1, yi = 1 ,
where i runs, sequentially, through all the values {1, . . . , |J|}; cf. the ILS method
proposed in [17]. The advantage of this approach is that the overconstrained
interval systems have (one) more equation than the original overconstrained
system, and hence the resulting necessary condition could be stronger. Our
numerical results discussed in Section 5 concern this variant. As a solver for
interval systems we utilize the convex approximation approach by Beaumont
[5]; it is sufficiently fast and produces narrow enough approximations of the
solution set.
RR n° 7544
Characterizing and approximating eigenvalue sets of symmetric interval matrices12
4.3.2
How to conclude for exact bounds?
Let us summarize properties of the submatrix vertex enumeration method. On
the one hand the worst case complexity of the algorithm is rather prohibitive,
O(3n), but on the other, we obtain better inner approximations, and sometimes
we get exact bounds of the eigenvalue sets. Theorem 1 and the discussion in the
previous section allow us to recognize exact bounds. For any i ∈{2, . . ., n}, we
have that if λi(AS) < λi−1(AS), then µi(AS) = λi(AS); a similar inequality
holds for the lower bound. This is a rather theoretical recipe because we may
not know a priori whether the assumption is satisfied. However, we can propose
a sufficient condition: If ωi(AS) < ωi−1(AS), then the assumption is obviously
true, and we conclude µi(AS) = λi(AS); otherwise we cannot conclude.
This sufficient condition is another reason why we need a sharp outer ap-
proximation. The sharper it is, the more times we are able to conclude that the
exact bound is achieved.
Exploiting the condition we can also decrease running time of submatrix
vertex enumeration.
We call Algorithm 3 only for p ∈{1, . . . , n} such that
p = 1 or ωp(AS) < ωp−1(AS). The resulting inner approximation may be a
bit less tight, but the number of exact boundary points of Λ(AS) that we can
identify remains the same.
Notice that there is enough open space for developing better conditions. For
instance, we do not know whether µi(AS) < µi−1(AS) (computed by submatrix
vertex enumeration) can serve also as a sufficient condition for the purpose of
determining exact bounds.
5
Numerical experiments
In this section we present some examples and numerical results illustrating prop-
erties of the proposed algorithms.
The experiments we performed on a PC
Intel(R) Core 2, CPU 3 GHz, 2 GB RAM, and the source code was written in
C++. We use GLPK v.4.23 [24] for solving linear programs, CLAPACK v.3.1.1
for its linear algebraic routines, and PROFIL/BIAS v.2.0.4 [21] for interval
arithmetic and basic operations. Notice, however, that routines of GLPK and
CLAPACK[1] do not produce verified solutions; for real-life problems this may not
be acceptable.
Example 1. Consider the following symmetric interval matrix
AS =


1
2
[1, 5]
2
1
1
[1, 5]
1
1


S
.
Local improvement (Algorithm 1) yields an inner approximation
µ1(AS) = [3.7321, 6.7843],
µ2(AS) = [0.0888, 0.3230],
µ3(AS) = [−4.1072, −1.0000].
The same result is obtained by the vertex enumeration (Algorithm 2). Therefore,
µ1(AS) = λ1(AS) and µ3(AS) = λ3(AS).
An outer approximation that is
RR n° 7544
Characterizing and approximating eigenvalue sets of symmetric interval matrices13
needed by the submatrix vertex enumeration (Algorithm 3) is computed using
the methods of Hladík et al. [18, 19]. It is
ω1(AS) = [3.5230, 6.7843],
ω2(AS) = [0.0000, 1.0519],
ω3(AS) = [−4.1214, −0.2019].
Now, the submatrix vertex enumeration algorithm yields the inner approxima-
tion
µ′
1(AS) = [3.7321, 6.7843],
µ′
2(AS) = [0.0000, 0.3230],
µ′
3(AS) = [−4.1072, −1.0000].
Since the outer approximation intervals do not overlap, we can conclude that
this approximation is exact, that is, λi(AS) = µ′
i(AS), i = 1, 2, 3.
This example shows two important aspects of the interval eigenvalue prob-
lem. First, it demonstrates that the vertex enumeration does not produce exact
bounds in general. Second, the symmetric eigenvalue set can be a proper subset
of the unsymmetric one, i.e., Λ(AS) ⫋Λ(A). This could be easily seen in the
matrix


1
2
1
2
1
1
5
1
1

.
It has three real eigenvalues 4.6458, −0.6458 and −1.0000, but the second one
does not belong to Λ(AS). Indeed, using the method by Hladík et al. [17] we
obtain
Λ(A) = [3.7321, 6.7843] ∪[−0.6458, 0.3230] ∪[−4.1072, −1.0000].
Example 2. Consider the example given by Qiu et al. [27] (see also [18, 35]):
AS =




[2975, 3025]
[−2015, −1985]
0
0
[−2015, −1985]
[4965, 5035]
[−3020, −2980]
0
0
[−3020, −2980]
[6955, 7045]
[−4025, −3975]
0
0
[−4025, −3975]
[8945, 9055]




S
.
The local improvement (Algorithm 1) yields an inner approximation
µ1(AS) = [12560.8377, 12720.2273], µ2(AS) = [7002.2828, 7126.8283],
µ3(AS) = [3337.0785, 3443.3127],
µ4(AS) = [842.9251, 967.1082].
The vertex enumeration (Algorithm 2) produces the same result. Hence we can
state that µ1(AS) and µ4(AS) are optimal.
To call the last method, submatrix vertex enumeration (Algorithm 3) we
need an outer approximation. We use the following by [18]
ω1(AS) = [12560.6296, 12720.2273], ω2(AS) = [6990.7616, 7138.1800],
ω3(AS) = [3320.2863, 3459.4322], ω4(AS) = [837.0637, 973.1993].
RR n° 7544
Characterizing and approximating eigenvalue sets of symmetric interval matrices14
Now, submatrix vertex enumeration yields the same inner approximation as the
previous methods. However, now we have more information. Since the outer
approximation interval are mutually disjoint, the obtained results are the best
possible. Therefore, µi(AS) = λi(AS), where i = 1, . . . , 4.
Example 3. Herein, we present two examples for approximating the singular
values of an interval matrix.
Let A ∈Rm×n and q := min{m, n}.
By the
Jordan–Wielandt theorem [13, 20, 25] the singular values, σ1(A) ≥· · · ≥σq(A),
of A are identical with the q largest eigenvalues of the symmetric matrix

0
AT
A
0

.
Thus, if we consider the singular value sets σ1(A), . . . , σq(A) of some interval
matrix A ∈Rm×n, we can identify them as the q largest eigenvalue sets of the
symmetric interval matrix
M :=

0
AT
A
0
S
.
(1) Consider the following interval matrix from [9].
A =


[2, 3]
[1, 1]
[0, 2]
[0, 1]
[0, 1]
[2, 3]


Both the local improvement and the vertex enumeration result the same inner
approximation, i.e.
µ1(M) = [2.5616, 4.5431],
µ2(M) = [1.2120, 2.8541].
Thus, σ1(A) = 4.5431. Additionally, consider the following outer approximation
from [18].
ω1(M) = [2.0489, 4.5431],
ω2(M) = [0.4239, 3.1817].
Using Algorithm 3, we obtain
µ′
1(M) = [2.5616, 4.5431],
µ′
2(M) = [1.0000, 2.8541].
Now we can claim that σ2(A) = 1, since ω2(M) > 0. Unfortunately, we cannot
conclude about the exact values of the remaining quantities, since the two outer
approximation intervals overlap. We know only that σ1(A) ∈[2.0489, 2.5616]
and σ2(A) ∈[2.8541, 3.1817].
(2) The second example comes from Ahn & Chen [2]. Let A be the following
interval matrix
A =


[0.75, 2.25]
[−0.015, −0.005]
[1.7, 5.1]
[3.55, 10.65]
[−5.1, −1.7]
[−1.95, −0.65]
[1.05, 3.15]
[0.005, 0.015]
[−10.5, −3.5]

.
Both local improvement and vertex enumeration yield the same result, i.e.
µ1(M) = [4.6611, 13.9371],
µ2(M) = [2.2140, 11.5077],
µ3(M) = [0.1296, 2.9117].
RR n° 7544
Characterizing and approximating eigenvalue sets of symmetric interval matrices15
Hence, σ1(A) = 13.9371.
As an outer approximation we use the following
intervals, using [18].
ω1(M) = [4.3308, 14.0115],
ω2(M) = [1.9305, 11.6111],
ω3(M) = [0.0000, 5.1000].
Running the submatrix vertex enumeration, we get the inner approximation
µ′
1(M) = [4.5548, 13.9371],
µ′
2(M) = [2.2140, 11.5077],
µ′
3(M) = [0.1296, 2.9517].
We cannot conclude that σ3(A) = µ3(A) = 0.1296, because ω3(M) has a
nonempty intersection with the fourth largest eigenvalue set, which is equal
to zero. Also the other singular value sets remain uncertain, but within the
computed inner and outer approximation.
Notice that µ′
1(M) < µ1(M), whence µ′
1(M) < λ1(M) = σ1(A) disproving
the Hertz’s Theorem 1 from [16] that the lower and upper limits of λ1(M) and
λn(M) are computable by the vertex enumeration method. It is true only for
λ1(M) and λn(M).
Example 4. In this example we present some randomly generated examples of
large dimensions. The entries of the midpoint matrix, Ac, are taken randomly in
[−20, 20] using the uniform distribution. The entries of the radius matrix A∆are
taken randomly, using the uniform distribution in [0, R], where R is a positive
real number. We applied our algorithm on the interval matrix M := AT A,
because it has a convenient distribution of eigenvalue set—some are overlapping
and some are not. Sharpness of results is measured using the quantity
1 −eT µ∆(M S)
eT ω∆(M S),
where e = (1, . . . , 1)T . This quantity lies always within the interval [0, 1]. The
closer to zero it is, the tighter the approximation. In addition, if it is zero, then
we achieved exact bounds. The initial outer approximation, ωi(M S), where
1 ≤i ≤n, was computed using the method due of Hladík et al.
[18], and
filtered by the method proposed by Hladík et al. in [19]. Finally, it was refined
according to the comment in Section 4.3.2. For the submatrix vertex enumer-
ation algorithm we implemented the branch & bound improvement, which is
described in Section 4.3.1 and 4.3.2.
The results are displayed in Table 1.
Example 5. In this example we present some numerical results on approxi-
mating singular value sets as introduced in Example 3. The input consists of an
interval (rectangular) matrix A ⊆Rm×n which is selected randomly as in the
previous example.
Table 2 presents our experiments. The time in the table corresponds to the
computation of the approximation of only the q largest eigenvalue sets of the
Jordan–Wielandt matrix.
RR n° 7544
Characterizing and approximating eigenvalue sets of symmetric interval matrices16
n
R
Algorithm 1
Algorithm 2
Algorithm 3
sharpness
time
sharpness
time
sharpness
time
5
0.001
0.05817
0.00 s
0.05041
0.00 s
0.00000
0.04 s
5
0.01
0.07020
0.00 s
0.05163
0.00 s
0.00000
0.03 s
5
0.1
0.26273
0.00 s
0.23389
0.00 s
0.17332
0.04 s
5
1
0.25112
0.00 s
0.23644
0.00 s
0.20884
0.01 s
10
0.001
0.08077
0.00 s
0.07412
0.09 s
0.00000
1.15 s
10
0.01
0.13011
0.01 s
0.11982
0.08 s
0.04269
1.29 s
10
0.1
0.27378
0.01 s
0.25213
0.09 s
0.12756
3.17 s
10
1
0.56360
0.01 s
0.52330
0.09 s
0.52256
2.58 s
15
0.001
0.07991
0.02 s
0.07557
7.3 s
0.00000
16.47 s
15
0.01
0.21317
0.02 s
0.19625
6.5 s
0.11341
2 m 29 s
15
0.1
0.36410
0.02 s
0.34898
7.0 s
0.34869
4 m 58 s
15
1
0.76036
0.02 s
0.73182
7.2 s
0.73182
7.5 s
20
0.001
0.09399
0.06 s
0.09080
7 m 21 s
0.00000
13 m 46 s
20
0.01
0.24293
0.06 s
0.22976
7 m 6 s
0.12574
1 h 14 m 55 s
20
0.1
0.24293
0.06 s
0.22976
7 m 14 s
0.12574
1 h 15 m 41 s
20
1
0.82044
0.06 s
0.79967
7 m 33 s
0.79967
7 m 39 s
25
0.001
0.14173
0.13 s
0.13397
6 h 53 m 0 s
0.02871
9 h 32 m 54 s
Table 1: Eigenvalues of random interval symmetric matrices AT A of dimension
n × n.
m
n
R
Algorithm 1
Algorithm 2
Algorithm 3
sharpness
time
sharpness
time
sharpness
time
5
5
0.01
0.08945
0.00 s
0.07716
0.10 s
0.00000
0.53 s
5
5
0.1
0.09876
0.01 s
0.09270
0.08 s
0.00000
0.73 s
5
5
1
0.43560
0.01 s
0.31419
0.10 s
0.26795
4.34 s
5
10
0.01
0.11320
0.02 s
0.10337
5.79 s
0.00000
7.91 s
5
10
0.1
0.13032
0.02 s
0.12321
5.98 s
0.00000
8.40 s
5
10
1
0.35359
0.02 s
0.33176
5.52 s
0.22848
21.53 s
5
15
0.01
0.10603
0.05 s
0.09424
5 m 31 s
0.00000
5 m 36 s
5
15
0.1
0.17303
0.04 s
0.16758
5 m 33 s
0.00000
7 m 58 s
5
15
1
0.46064
0.05 s
0.39708
5 m 32 s
0.31847
15 m 47 s
10
10
0.01
0.10211
0.06 s
0.09652
8 m 3 s
0.00000
8 m 19 s
10
10
0.1
0.13712
0.07 s
0.13387
8 m 10 s
0.00000
14 m 12 s
10
10
1
0.39807
0.07 s
0.35580
7 m 52 s
0.30279
26 h 48 m 38 s
10
15
0.01
0.09561
0.12 s
0.09116
5 h 51 m 53 s
0.00000
5 h 54 m 56 s
Table 2: Singular values of random interval matrices of dimension m × n.
RR n° 7544
Characterizing and approximating eigenvalue sets of symmetric interval matrices17
6
Conclusion and future directions
We proposed a new solution theorem for the symmetric interval eigenvalue prob-
lem, which describes some of the boundary points of the eigenvalue set. Unfor-
tunately the complete characterisation is still a challenging open problem.
We developed an inner approximation algorithm (submatrix vertex enumer-
ation), which in the case where the eigenvalue sets are disjoint, and the interme-
diate gaps are wide enough, output exact results. To our knowledge, even under
this assumption, this is the first algorithm that can guarantee exact bounds.
Based on our numerical experiments suggest that the local search algorithm
is superior to the submatrix vertex enumeration algorithm when the input ma-
trices are not of very small dimension.
Acknowledgment.
ET is partially supported by an individual postdoctoral
grant from the Danish Agency for Science, Technology and Innovation.
References
[1] CLAPACK
–
Linear
Algebra
PACKage
written
for
C.
http://www.netlib.org/clapack/.
[2] H.-S. Ahn and Y. Q. Chen. Exact maximum singular value calculation of
an interval matrix. IEEE Trans. Autom. Control, 52(3):510–514, 2007.
[3] G. Alefeld and J. Herzberger. Introduction to interval computations. Aca-
demic Press, London, 1983.
[4] G. Alefeld and G. Mayer. Interval analysis: Theory and applications. J.
Comput. Appl. Math., 121(1-2):421–464, 2000.
[5] O. Beaumont.
Solving interval linear systems with linear programming
techniques. Linear Algebra and Its Applications, 281(1-3):293–309, 1998.
[6] O. Beaumont. An algorithm for symmetric interval eigenvalue problem.
Technical Report IRISA-PI-00-1314, Institut de recherche en informatique
et systèmes aléatoires, Rennes, France, 2000.
[7] D. Chablat, P. Wenger, F. Majou, and J. Merlet.
An interval analysis
based study for the design and the comparison of 3-dof parallel kinematic
machines. Int. J. Robot. Res., 23(6):615–624, 2004.
[8] A. Deif and J. Rohn.
On the invariance of the sign pattern of matrix
eigenvectors under perturbation. Linear Algebra Appl., 196:63–70, 1994.
[9] A. S. Deif. Singular values of an interval matrix. Linear Algebra Appl.,
151:125–133, 1991.
[10] A. S. Deif.
The interval eigenvalue problem.
Z. Angew. Math. Mech.,
71(1):61–64, 1991.
[11] A. D. Dimarogonas. Interval analysis of vibrating systems. J. Sound Vib.,
183(4):739–749, 1995.
RR n° 7544
Characterizing and approximating eigenvalue sets of symmetric interval matrices18
[12] F. Gioia and C. N. Lauro. Principal component analysis on interval data.
Comput. Stat., 21(2):343–363, 2006.
[13] G. H. Golub and C. F. Van Loan. Matrix computations. 3rd ed. Johns
Hopkins University Press, 1996.
[14] E. Hansen and G. W. Walster. Global optimization using interval analysis.
2nd ed., revised and expanded. Marcel Dekker, New York, 2004.
[15] D. Hertz. The extreme eigenvalues and stability of real symmetric interval
matrices. IEEE Trans. Autom. Control, 37(4):532–535, 1992.
[16] D. Hertz. Interval analysis: Eigenvalue bounds of interval matrices. In
C. A. Floudas and P. M. Pardalos, editors, Encyclopedia of optimization.,
pages 1689–1696. Springer, New York, 2009.
[17] M. Hladík, D. Daney, and E. Tsigaridas. An algorithm for the real in-
terval eigenvalue problem.
Research Report RR-6680, INRIA, France,
http://hal.inria.fr/inria-00329714/en/, October 2008. sumbitted to
J. Comput. Appl. Math.
[18] M. Hladík, D. Daney, and E. Tsigaridas. Bounds on eigenvalues and sin-
gular values of interval matrices. Research Report inria-00370603:1, IN-
RIA, France, http://hal.inria.fr/inria-00370603/en/, March 2009.
sumbitted to SIAM J. Matrix Anal. Appl.
[19] M. Hladík, D. Daney, and E. Tsigaridas. A filtering method for the in-
terval eigenvalue problem.
Research Report RR-7057, INRIA, France,
http://hal.inria.fr/inria-00422966/en/, October 2009. sumbitted to
Computing.
[20] R. A. Horn and C. R. Johnson. Matrix analysis. Cambridge University
Press, Cambridge, 1985.
[21] O. Knüppel, D. Husung, and C. Keil. PROFIL/BIAS – a C++ class library.
http://www.ti3.tu-harburg.de/Software/PROFILEnglisch.html.
[22] L. V. Kolev. Outer interval solution of the eigenvalue problem under general
form parametric dependencies. Reliab. Comput., 12(2):121–140, 2006.
[23] H. Leng and Z. He.
Computing eigenvalue bounds of structures with
uncertain-but-non-random parameters by a method based on perturbation
theory. Commun. Numer. Methods Eng., 23(11):973–982, 2007.
[24] A.
Makhorin.
GLPK
–
GNU
Linear
Programming
Kit.
http://www.gnu.org/software/glpk/.
[25] C. D. Meyer. Matrix analysis and applied linear algebra. SIAM, Philadel-
phia, 2000.
[26] A. Neumaier. Interval methods for systems of equations. Cambridge Uni-
versity Press, Cambridge, 1990.
[27] Z. Qiu, S. Chen, and I. Elishakoff. Bounds of eigenvalues for structures with
an interval description of uncertain-but-non-random parameters.
Chaos
Solitons Fractals, 7(3):425–434, 1996.
RR n° 7544
Characterizing and approximating eigenvalue sets of symmetric interval matrices19
[28] Z. Qiu, P. C. Müller, and A. Frommer. An approximate method for the
standard interval eigenvalue problem of real non-symmetric interval matri-
ces. Commun. Numer. Methods Eng., 17(4):239–251, 2001.
[29] Z. Qiu and X. Wang. Solution theorems for the standard eigenvalue prob-
lem of structures with uncertain-but-bounded parameters. J. Sound Vib.,
282(1-2):381–399, 2005.
[30] J. Rohn.
Interval matrices: Singularity and real eigenvalues.
SIAM J.
Matrix Anal. Appl., 14(1):82–91, 1993.
[31] J. Rohn. An algorithm for checking stability of symmetric interval matrices.
IEEE Trans. Autom. Control, 41(1):133–136, 1996.
[32] J.
Rohn.
A
handbook
of
results
on
interval
linear
problems.
http://www.cs.cas.cz/rohn/handbook, 2005.
[33] J. Rohn. Solvability of systems of interval linear equations and inequalities.
In M. Fiedler, J. Nedoma, J. Ramík, J. Rohn, and K. Zimmermann, editors,
Linear optimization problems with inexact data., chapter 2, pages 35–77.
Springer, New York, 2006.
[34] J. Rohn and A. Deif. On the range of eigenvalues of an interval matrix.
Comput., 47(3-4):373–377, 1992.
[35] Q. Yuan, Z. He, and H. Leng. An evolution strategy method for computing
eigenvalue bounds of interval matrices. Appl. Math. Comput., 196(1):257–
265, 2008.
Contents
1
Introduction
3
2
Basic definitions and theoretical background
4
3
Main theorem
5
4
Inner approximation algorithms
7
4.1
Local improvement . . . . . . . . . . . . . . . . . . . . . . . . . .
7
4.2
Vertex enumeration . . . . . . . . . . . . . . . . . . . . . . . . . .
8
4.3
Submatrix vertex enumeration
. . . . . . . . . . . . . . . . . . .
8
4.3.1
Branch & bound improvement
. . . . . . . . . . . . . . .
11
4.3.2
How to conclude for exact bounds? . . . . . . . . . . . . .
12
5
Numerical experiments
12
6
Conclusion and future directions
17
RR n° 7544
Centre de recherche INRIA Sophia Antipolis – Méditerranée
2004, route des Lucioles - BP 93 - 06902 Sophia Antipolis Cedex (France)
Centre de recherche INRIA Bordeaux – Sud Ouest : Domaine Universitaire - 351, cours de la Libération - 33405 Talence Cedex
Centre de recherche INRIA Grenoble – Rhône-Alpes : 655, avenue de l’Europe - 38334 Montbonnot Saint-Ismier
Centre de recherche INRIA Lille – Nord Europe : Parc Scientifique de la Haute Borne - 40, avenue Halley - 59650 Villeneuve d’Ascq
Centre de recherche INRIA Nancy – Grand Est : LORIA, Technopôle de Nancy-Brabois - Campus scientifique
615, rue du Jardin Botanique - BP 101 - 54602 Villers-lès-Nancy Cedex
Centre de recherche INRIA Paris – Rocquencourt : Domaine de Voluceau - Rocquencourt - BP 105 - 78153 Le Chesnay Cedex
Centre de recherche INRIA Rennes – Bretagne Atlantique : IRISA, Campus universitaire de Beaulieu - 35042 Rennes Cedex
Centre de recherche INRIA Saclay – Île-de-France : Parc Orsay Université - ZAC des Vignes : 4, rue Jacques Monod - 91893 Orsay Cedex
Éditeur
INRIA - Domaine de Voluceau - Rocquencourt, BP 105 - 78153 Le Chesnay Cedex (France)
http://www.inria.fr
ISSN 0249-6399
apport 

     technique