An interval-matrix branch-and-bound algorithm for bounding eigenvalues

We present and explore the behaviour of a branch-and-bound algorithm for calculating valid bounds on the kth largest eigenvalue of a symmetric interval matrix. Branching on the interval elements of the matrix takes place in conjunction with the application of Rohn's method (an interval extension of Weyl's theorem) in order to obtain valid outer bounds on the eigenvalues. Inner bounds are obtained with the use of two local search methods. The algorithm has the theoretical property that it provides bounds to any arbitrary precision (assuming infinite precision arithmetic) within finite time. In contrast with existing methods, bounds for each individual eigenvalue can be obtained even if its range overlaps with the ranges of other eigenvalues. Performance analysis is carried out through nine examples. In the first example, a comparison of the efficiency of the two local search methods is reported using 4000 randomly generated matrices. The eigenvalue bounding algorithm is then applied to five randomly generated matrices with overlapping eigenvalue ranges. Valid and sharp bounds are indeed identified given a sufficient number of iterations. Furthermore, most of the range reduction takes place in the first few steps of the algorithm so that significant benefits can be derived without full convergence. Finally, in the last three examples, the potential of the algorithm for use in algorithms to identify index-1 saddle points of nonlinear functions is demonstrated.


Introduction
In many practical applications requiring the computation of eigenvalues, the matrix of interest is known only as a function of some parameters and is therefore often expressed as an interval matrix [4,10,18]. As a result, there is a need for methods that allow the calculation or estimation of the ranges of the eigenvalues of interval matrices. However in general, problems associated with the eigenvalues of interval matrices are difficult problems. For example, checking positive-(semi)definiteness [18,22] or regularity (existence of singular matrix) [20] of interval matrices are known to be NP-hard problems. Moreover, computing approximate solutions for the minimum and maximum eigenvalues of symmetric interval matrices can be NP-hard [8].
Eigenvalue bounding methods also play an important role in deterministic global optimization algorithms. They are used in order to create valid convex underestimators [1] of general nonlinear functions. Furthermore, methods for bounding the lowest and second lowest eigenvalues of a symmetric interval matrix can be used as a test for identifying domains in which a twice-continuously differentiable function contains (or does not contain) index-1 saddle points. This can be used within a global deterministic algorithm to speed up the location of index-1 saddle points of potential energy functions [19], a challenging problem with applications in chemical engineering and other areas [5].
A number of methods have been proposed in the literature to obtain lower and upper bounds on the smallest and largest eigenvalues, respectively, of interval matrices [1,7,9,[23][24][25]. Other methods have been devised to compute bounds for each individual eigenvalue [9,11,12,21]. An evolutionary method approach for inner bounds was presented by Yuan et al. [27]. Exact bounds for individual eigenvalues have been given by Deif [3] provided that the signs of the eigenvector entries remain constant over the interval matrix. This condition limits the applicability of this result.
The algorithms by Hladík et al. [10], Leng et al. [15], and Leng [14] can be used to calculate the real eigenvalue set of an interval matrix with any given precision. These algorithms begin with the calculation of an initial inclusion set and proceed by successive identification and removal of parts of the initial inclusion set which do not belong to the eigenvalue set. In particular, the algorithm by Hladík et al. has been shown to be fast and applicable to very large matrices (with small interval widths). However, when the ranges of individual eigenvalues overlap, the methods in [10,14,15] can only provide, at best, the bounds of the union of the overlapping ranges.
We present a branch-and-bound algorithm for the calculation of the bounds of any individual eigenvalue of symmetric interval matrices. The branching occurs on the interval entries of the input matrix. We use Ronh's theorem [9,24], which is an interval extension of Weyl's theorem [6] and local improvement methods in order to obtain valid bounds at each step. The algorithm can be used to calculate the bounds of a specific eigenvalue regardless of whether its range overlaps with that of other eigenvalues or not. Furthermore, the algorithm does not necessarily require the use of interval arithmetic.
The paper is organized as follows: In Section 2 we give a brief introduction to interval matrices and a few definitions. In Section 3 we present the pseudocode of the interval-matrix branchand-bound algorithm. In Section 4 we present the general bounding approach used in the main algorithm. In Section 5 we present two local search algorithms. One from the existing literature and one given here. These algorithms are used in the main algorithm in order to improve the inner bounds and speed up convergence. In Section 6 we present results from the application of the method and finally and in Section 7, we draw our conclusions.

Preliminaries
We denote interval variables with lower case letters inside square brackets, [x], and the corresponding lower and upper bounds of [x] as x and x, respectively. Symmetric interval matrices are denoted by capital letters inside square brackets. An interval matrix is simply a matrix with interval instead of scalar entries. For example, a 2 × 2 symmetric interval matrix is . For an n × n symmetric scalar matrix M, we denote by λ i (M) the ith largest eigenvalue of M. Therefore, we order the eigenvalues as λ n (M) ≤ λ n−1 (M) ≤ · · · ≤ λ 1 (M). We will make use of the following definitions:

The interval-matrix branch-and-bound algorithm
In this section we introduce the interval-matrix branch-and-bound algorithm which follows a 'classic' branch-and-bound scheme [13]. Given an n × n symmetric interval matrix , the order of the eigenvalue (e.g. kth largest eigenvalue), k, for which the bounds will be calculated, the maximum number of iterations, maxiters and the precision . We denote by L the list which contains sublists of the form is a symmetric matrix with l i and u i lower and upper bounds of λ k ([M i ]). Since the branching procedure can be represented by a binary tree, we will refer to the sublists in L as nodes. Furthermore, we denote the best lower and upper bounds by BLB and BUB, respectively. Finally, we denote by RLB the lowest lower bound of the nodes that have been removed due to the fact that the required precision has been achieved (e.g. the difference between the lower bound at a node and BUB is less than ).
The choice of branching strategy (in our case which entry we choose to branch on in step 10 of Algorithm 1) can have a strong influence on the performance of branch-and-bound algorithms (see, for example, [2]). As will be evident from Proposition 4.3 in Section 4, in order to achieve theoretical convergence for a given precision, , it is necessary to branch on all (off-diagonal) interval entries. A straightforward branching scheme that meets this requirement would be to branch on the interval with the maximum width. However, since the lower bound is given by λ k (C) − [E] ∞ , a more judicious choice might be to branch on the interval entry which reduces [E] ∞ the most. However, from our experiments with the algorithm, we have not observed any significant effects on performance due to these different branching schemes. Finally, note that in Algorithm 1, the node which is visited at each step is the one with the current lowest lower bound.

M∈[M]
λ k (M) can be calculated in an analogous way to Algorithm 1. In the following section, we discuss the bounding steps and branching of the algorithm in more detail. In the analysis which follows we will assume infinite precision arithmetic. Choose the first entry, L 1 , from list L.

General bounding approach
We can write any given interval matrix We can write Equation (1) as Thus, we can use the following for calculating bounds on λ k ([M]): where b is a lower bound of −ρ([E]). Since the value λ k (C) is attained for C ∈ [M] this means that λ k ([M]) cannot be greater than λ k (C), and thus λ k (C) is a valid upper bound. However, a better upper bound can be achieved with the use of local improvement step(s). We will see more on this matter later in this section. On the other hand, we can calculate a lower bound of −ρ([E]) using a number of methods (see, e.g. [1]) or even the exact value of ρ([E]) with the Hertz method [7] (O(2 n−1 )). We will use (This bound is the same as the one we would get by the interval Gerschgorin method [1]). Note that in steps 12 and 13 the new lower bounds l 1 , l 2 can actually be worse (lower) than l. In such case we simply replace them with l. The reason for this is that λ k (C 1 ) and/or λ k (C 2 ) could be less than λ k (C) and at the same Next, the following lemma tells us that there is no need to branch on the diagonal entries.
Then, given accuracy , by successive bisection every submatrix will have u − l ≤ after a total of (( Proof It is helpful to imagine the bisection process as a binary tree. What we would like is to know how many bisections we need to perform in order to have a full binary tree where for each leaf of the tree the corresponding submatrix In order to halve [E] ∞ once, we need to branch on each of the (n 2 − n)/2 off-diagonal elements. Which means that the depth of the tree will grow to d = (n 2 − n)/2. Therefore, in order to achieve the required accuracy, we will need to end up with a tree of depth d ≥ ((n 2 − n)/2) log 2 ((n − 1)w/2 ). The number of leaves of a full binary tree is 2 d and the corresponding number of bisections is 2 d − 1. Thus the required number of iterations would be At this point we make the following notes: first, although the number of iterations given by (9) is forbiddingly high, this is a worst-case scenario where no node fathoming takes place. As it can be seen in Section 6, the actual performance appears better than (9). This suggests and renders the algorithm practical for small-sized matrices with wide intervals. Furthermore, the actual dimension of the problem depends on the number of interval off-diagonal entries in the matrix and thus can be less than (n 2 − n)/2. Second, to the best of our knowledge, no other method exists that can solve the problem of calculating λ k and/or λ k for k ∈ {2, . . . , n − 1} for general interval symmetric matrices, for a given accuracy . Thus, we do not know of any better upper bound on the algorithmic complexity of the problem.

Local search algorithms
As mentioned in the previous section, we can improve the trivial upper bound, For this purpose we present two local search methods. The first is by Hladík et al. [11], and the second introduced here, based on Theorem 5.1 given in Section 5.2.

Hladík et al. [11] local search
In [11], Hladík et al. proposed a number of local search methods which they used along with eigenvalue bounding methods in order to obtain inner and outer approximations of eigenvalue ranges. Although in general the values obtained from the local search methods are approximate (conservative), they can be sometimes shown to be exact [11]. Here, we will make use of the method, from [11], shown in Algorithm 2 where with C we denote the centre matrix of [M] and with E being the radius matrix (e.g. e ij = (m ij − m ij )/2 ). By v k (M) we denote the eigenvector of M which corresponds to the kth largest eigenvalue.  6: For a local search of the maximum of λ k (i.e. valid lower bound of λ k ), we just replace λ k = −∞, < with > and − with + in Steps 2, 3, and 6, respectively.

A new local search method
It is known that the values λ n and λ 1 of a symmetric interval matrix [M] are attained at extreme matrices of [M] [7]. However, for the rest of the eigenvalue bounds (λ n , λ 1 , and λ k , λ k for k = 2, . . . , n − 1) of a general symmetric interval matrix, the situation is more complicated since boundary values can be attained in I( [M]) and/or B([M]). Nevertheless, for matrices with non-zero width off-diagonal elements, we can prove the following necessary conditions for ). Note that in the following theorem, the request for the diagonal elements to have zero widths is not an extra assumption and this stems from Lemma 4.2. Otherwise, it would be meaningless to refer to solutions in the interior, since that would never be true. Furthermore, for the remainder of this section, we assume eigenvectors have been normalized with respect to the Euclidean norm. are being the square of the ith entry of v k ) and consider the matrix in an attempt to improve the upper bound, in the same way as with Algorithm 2. This is detailed in Algorithm 3.

Algorithm 3
A new local improvement algorithm.
For a local search of the maximum of λ k (i.e. lower bound of λ k ), we replace in Steps 5 and 7.
Notice that Algorithm 2 can be applied to any kind of symmetric matrix while Algorithm 3 requires the off-diagonal entries to have non-zero widths. A common feature is that both algorithms, most of the times, terminate after one step and only rarely after two or three steps. Furthermore, Algorithm 2 always searches extreme matrices which makes it more suitable for obtaining an upper bound on λ n and a lower bound on λ 1 . A comparison between the bounds obtained by the two algorithms is given in Example 1 in Section 6.
In practice we use both local search algorithms and we simply keep the best result. However, we do not apply them in every bounding step of the main algorithm, in order to reduce computational time. The methods are applied at the initial step and then we proceed at each iteration by using the bound obtained by the centre matrix eigenvalue calculation. If this bound happens to be better than the current best upper (or lower if we are bounding λ k ) bound we then apply Algorithms 2 and 3 for further improvement.
We can make a stronger statement than the one of Theorem 5.1 by observing the following: If λ k = λ k+1 , λ k+1 = λ k+2 , and m = v k 2 ∞ + v k+1 2 ∞ < 1 then, for adequately small ρ, . More formally, we have The eigenpairs of M 1 are (v i , λ i + ρm) for i < k and (v k , λ k+i − ρ(1 − m)) for i = 0, 1, . . . , s. Following the same reasoning as in the proof of Theorem 5.1, for adequately small ρ the matrix Remark: We could make use of Theorem 5.2 in order to expand the local search method given in Algorithm 3 but this would make the algorithm more complicated without providing any significant benefit in performance. This is due to the fact that the method usually terminates because a bound on one of the interval elements of the interval matrix is reached rather than because λ k = λ k+1 . Nevertheless, we have stated Theorem 5.2 for completeness.

Results
In this section we present the results from the application of the algorithm to a number of randomly generated symmetric interval matrices. Given dimension n and radius R, we obtain an interval matrix Note that λ n and λ 1 can be computed much faster, as mentioned previously, by the method from [7]. Nevertheless, we compute the extreme eigenvalue bounds for completeness. The algorithm is implemented in Python 2.7 and the calculations are performed with an Intel Core i7-3770 CPU @ 3.40GHz. Note that the calculations were not performed in a numerically verified way (e.g. use of interval arithmetic) which may lead to numerical errors in the bounds. If required, an implementation based on numerically validated bounds can be obtained with minor effort. The bounds in the results and the interval matrices used in each example are outwardly rounded to the third decimal place (this is why the bounds for λ 1 in Table 1 appear to have width > 10 −1 yet are marked as converged).

Example 1 -comparison of the two local improvement algorithms
In this example we make a comparison between the two local improvement methods, Algorithms 2 and 3. We generate four groups of random 5 × 5 matrices consisting of a thousand matrices each and with radii R = 5, 10, 20, and 30, respectively. We apply each algorithm in order to obtain upper bounds on λ k for k = 1,2,3,4 and lower bounds on λ k for k = 2,3,4,5. In Figure 1, for each group of random matrices, we plot a histogram of the difference between the bounds from the two algorithms. Positive values indicate that the bounds determined by Algorithm 3 were better than those obtained by Algorithm 2, while negative values indicate the reverse. More explicitly, we denote by λ (2) k U and λ (2) k L the upper and lower bounds on λ k Notes: An asterisk indicates the bound widths have converged to 10 −1 . I 10 and I 10 are the percentages of improvement, for the bounds of λ k and λ k , respectively, after 10 steps (Equations (15) and (16)), S is the percentage of relative sharpness of the final bounds (Equation (17)) and OBI is the percentage of outer bounds improvement (Equation (18)). As we see from Figure 1, for smaller radii, Algorithm 2 performs significantly better than Algorithm 3. However, as the radius increases, the relative performance of Algorithm 3 improves, as indicated by the increasing median value. Although on average Algorithm 2 performs better, an overall improvement can be achieved with the combination of the two methods (use of the best result). This comes at no significant cost since both methods terminate after one step in most cases and occasionally after two or three steps.

Examples 2-6 -performance of the eigenvalue bounding algorithm
The next set of examples is used to investigate the performance of Algorithm 1. For Examples 2-6, we plot the interval eigenvalue bounds as obtained by the algorithm and give a summary table of algorithmic performance. To explain the quantities in the table, we denote by λ (i) k L and λ (i) k U the lower and upper bounds, respectively, on λ k at iteration i. We also denote by λ (i) k L and λ (i) k U the lower and upper bounds, respectively, on λ k at iteration i. When the iteration superscript is omitted, the quantity refers to the corresponding value at the final iteration of the algorithm. In each table, we give the following information: the final bounds for λ k , λ k L , and λ k U ; the final bounds for λ k , λ k L , and λ k U ; the computational times required to obtain these bounds; the percentage improvement, denoted by I 10 and I 10 , of the bounds on λ k and λ k , respectively, between the initial and the tenth iterations: and the relative sharpness, S, of the final bounds as a percentage: the percentage improvement of the outer bounds between the initial bounds and final outer bounds, OBI:

Example 2: n = 3, R = 10
In this example we calculate eigenvalue bounds for the following randomly generated, 3 × 3, symmetric matrix with R = 10:  Figure 2, we plot the inner, outer and initial outer eigenvalue ranges. Note that all the intervals in Figure 2   Notes: An asterisk indicates the bound widths have converged to 10 −1 . I 10 and I 10 are the percentages of improvement, for the bounds of λ k and λ k , respectively, after 10 steps (Equations (15) and (16)), S is the percentage of relative sharpness of the final bounds (Equation (17)) and OBI is the percentage of outer bounds improvement (Equation (18)).

Example 3: n = 4, R = 5
In this example we increase the dimension of the random test matrix to 4 × 4 and reduce the radius to R = 5 and obtain the following matrix:  Notes: An asterisk indicates the bound widths have converged to 10 −1 . I 10 and I 10 are the percentages of improvement, for the bounds of λ k and λ k respectively, after 10 steps (Equations (15) and (16)), S is the percentage of relative sharpness of the final bounds (Equation (17)) and OBI is the percentage of outer bounds improvement (Equation (18)

Example 4: n = 5, R = 5
In this example we further increase the dimension of the random test matrix to 5 × 5 and maintain the radius to R = 5 to generate: Results and eigenvalue ranges for this example are shown in Table 3 and Figure 4, respectively. Figures 5-7 are indicative of the algorithm's behaviour. In the first case ( Figure 5), nodes are being removed at a high rate and thus convergence is achieved very fast. In the second case ( Figure 6), nodes are removed adequately fast and convergence to the required tolerance can still be achieved by increasing the maximum iteration number. In the third case (Figure 7), almost no nodes are removed and the algorithm does not converge. Nevertheless, a significant improvement with respect to the initial bounds is achieved at the termination of the algorithm. In Figure 8 we  show the progress of the bounds on λ 4 when we do not make use of local search algorithms. A comparison of this figure with Figure 6 demonstrates that the use of the local search can significantly increase the convergence speed.  In practical applications, the matrix might contain only a few interval entries. In this example, a random 7 × 7 symmetric matrix with a small number of interval entries (12 off-diagonal and 2 Notes: An asterisk indicates the bound widths have converged to 10 −1 . I 10 and I 10 are the percentages of improvement, for the bounds of λ k and λ k , respectively, after 10 steps (Equations (15) and (16)), S is the percentage of relative sharpness of the final bounds (Equation (17)) and OBI is the percentage of outer bounds improvement (Equation (18) diagonal entries randomly chosen with R = 5) is considered. The matrix is shown in Appendix 1.
Results are given in Table 4 and eigenvalue bounds are plotted in Figure 9. Although this matrix is larger than that in Example 4, the CPU time required is similar. In this example we apply Algorithm 1 to a 10 × 10 randomly generated tridiagonal symmetric matrix, again with R = 5. The matrix is shown in Appendix 2. Results are given in Table 5 and eigenvalue bounds are plotted in Figure 10. It can be seen that computational performance decreases for this larger matrix. Nevertheless, tight bounds are obtained on most of the eigenvalues. The results show that all 10 eigenvalues overlap despite the use of a relatively small radius of 5, and the bounds obtained are much tighter than the worst-case range of [−36.027, 34.841]. Notes: An asterisk indicates the bound widths have converged to 10 −1 . I 10 and I 10 are the percentages of improvement, for the bounds of λ k and λ k respectively, after 10 steps (Equations (15) and (16)), S is the percentage of relative sharpness of the final bounds (Equation (17)) and OBI is the percentage of outer bounds improvement (Equation (18)

Examples 7-9: identification of non index-1 areas
In many rate-based physical processes, such as chemical reactions or nucleation, the identification of transition states on potential energy surfaces is an important and challenging problem. Mathematically the problem is described as follows: Given a C 2 function f : B ⊂ R n → R find x * ∈ B such that ∇f (x * ) = 0 and ∇ 2 f (x * ) has 1 negative eigenvalue and n − 1 positive eigenvalues. The point x * is called a transition state or index-1 saddle point.
One emerging approach to this problem is based on the use of deterministic global search methods in which one aims first to locate all critical points and then to identify within those the physically relevant ones, namely index-1 saddle points (and possibly minima that are also of interest) [16,17,26]. However, in such an approach, it is necessary to obtain full convergence for all critical points, including those that are not of practical relevance. In order to focus the computational effort on the location of index-1 saddle points, a number of practical tests that allow the identification of non index-1 areas have recently been proposed [19]. These tests can Notes: 'Dim.' refers to the dimensionality of the example, L is the maximum edge length for the random hyper-rectangles, columns 5 and 6 indicate the number of hyper-rectangles, out of 10,000, that are found not to contain an index-1 saddle point, using Rohn's method and Algorithm 1, respectively, column 7 refers to the percentage increase in the number of hyper-rectangles found with the proposed approach be used in each iteration of a global deterministic search method and allow the exclusion of such areas, thereby avoiding further investigation and speeding up the convergence process. One such test is based on the use of Rohn's method. For example, if we find that λ n−1 ≤ 0 for the interval Hessian matrix derived from the function f and calculated over a hyper-rectangular area B, then no index-1 saddle points exist in B and thus B can be excluded from the search. Motivated by the fact that within only a few iterations the algorithm leads to a significant improvement of the initial eigenvalue bounds, as seen from the I 10 values in Tables 1-5, it might be advantageous to use a few steps of the algorithm instead of a single step (Rohn's method) for the exclusion of non index-1 areas. In Table 6 we show the results for three selected test functions. For each test function, we randomly create 10,000 hyper-rectangles, within a given domain, with each edge having a length selected randomly in the interval [0, L]. For each hyper-rectangular area, we first calculate the corresponding interval Hessian matrix. We compare the application of two approaches to obtain an upper bound on λ n−1 and determine whether the region may contain an index-1 saddle point: Rohn's method (i.e. the initial step of Algorithm 1) or Algorithm 1 for a maximum of only five steps. For each approach, we report the number of hyper-rectangle found such that λ n−1 ≤ 0.
where n = 5 and x ∈ [−5, 5] 5 The results indicate that a significant increase in the number of regions identified not to contain an index-1 saddle point is achieved when using Algorithm 1. Furthermore, the solution given by the local search method can be used in a stopping criterion since if the lower bound of λ n−1 is found to be strictly positive, there is no reason to proceed. In practice we would apply the local search only once and prior to branch-and-bound procedure for bounding λ n−1 . Finally, in the same way described above, we can use the algorithm to identify convex areas (λ n ≥ 0) and index-1 areas (λ n < 0 and λ n−1 > 0).

Conclusions
We have presented a branch-and-bound algorithm for calculating bounds on all individual eigenvalues of a symmetric interval matrix. The algorithm is based on calculating successively tighter bounds by branching on the off-diagonal interval entries of the input matrix using Rohn's method for outer bounds and local search methods for inner bounds. In contrast to other methods, the algorithm provides valid and distinct bounds on each eigenvalue, regardless of whether the ranges of the eigenvalues overlap. Application to five examples, up to a 10 × 10 matrix, has shown that the algorithm can achieve significant reductions in the range of each eigenvalue compared to existing methods. The use of local search methods has been found to increase the convergence speed significantly. Two approaches have been used for local search: one developed previously [11] and one proposed here. The [11] method is found to perform best on average, but not systematically, making the combination of these two fast approaches desirable.
The algorithm is particularly effective for low-dimensional problems, where by dimension we mean the number of interval entries in the initial matrix. While the algorithm becomes more computationally demanding for larger problems, a few iterations always yield substantial reductions in the eigenvalue ranges and provide a low-cost approach to obtain good bounds. Furthermore, as shown in Examples 7-9, the proposed algorithm can be used as an effective improvement over Rohn's method as a test in deterministic global search methods for the location of index-1 saddle points.

Data statement
Supporting data for this work are available on request by writing to the corresponding author (c.adjiman@imperial.ac.uk).

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
Financial support from the Engineering and Physical Sciences Research Council (EPSRC) of the UK, via a Leadership Fellowship [EP/J003840/1], is gratefully acknowledged.