{'fraction_non_alphanumeric': 0.038239650736727214, 'fraction_numerical': 0.010602635066656272, 'mean_word_length': 4.344791666666667, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Determinants of Interval Matrices

In this paper we shed more light on determinants of interval matrices. Computing the exact bounds on a determinant of an interval matrix is an NP-hard problem. Therefore, attention is first paid to approximations. NP-hardness of both relative and absolute approximation is proved. Next, methods computing verified enclosures of interval determinants and their possible combination with preconditioning are discussed. A new method based on Cramer's rule was designed. It returns similar results to the state-of-the-art method, however, it is less consuming regarding computational time. As a byproduct, the Gerschgorin circles were generalized for interval matrices. New results about classes of interval matrices with polynomially computable tasks related to determinant are proved (symmetric positive definite matrices, class of matrices with identity midpoint matrix, tridiagonal H-matrices). The mentioned methods were exhaustively compared for random general and symmetric matrices.


Introduction.
Interval determinants can be found in various applications. They were used e.g., in [24] for testing regularity of inverse Jacobian matrix, in [28] for workspace analysis of planar flexure-jointed mechanism, in [30] for computer graphics applications or in [40] as a testing tool for Chebyshev systems.
In this work we first address computational properties of determinants of general interval matrices. We are going to prove two new results regarding absolute and relative approximation of interval determinants. Next, we slightly mention known tools that can be used for computing interval determinants -interval Gaussian elimination, Hadamard inequality and Gerschgorin circles. We introduce our new method based both on Cramer's rule and solving interval linear systems. Regarding symmetric matrices, there are many results about enclosing their eigenvalues and they can be also used for computing interval determinants. All the methods work much better when combined with some kind of preconditioning. We briefly address that topic. We also prove that some classes of interval matrices have some tasks related to interval determinant computable in polynomial time (symmetric positive definite matrices, some matrices with identity midpoint matrix, tridiagonal H-matrices). At the end we provide thorough numerical testing of the mentioned methods on random general and symmetric interval matrices.
2. Basic notation and definitions. In our work it will be sufficient to deal only with square interval matrices. An interval matrix is defined by A = {A ∈ R n×n | A ≤ A ≤ A} for A, A ∈ R n×n such that A ≤ A (understood component-wise). To compute with intervals we use the standard interval arithmetic, for more details on the interval arithmetic see for example [25] or [27].
We denote intervals and interval structures in boldface (a, A, b). Real point matrices and vectors will be denoted in normal case (A, b). An interval coefficient of A lying at the position (i, j) is denoted by A ij .
An interval can be also defined by its midpoint a c ∈ R and radius a ∆ ∈ R as a = [a c − a ∆ , a c + a ∆ ]. Interval vectors and matrices are defined similarly. Notation mid(a), rad(a) can be sometimes used instead of a c , a ∆ respectively. The set of all real closed intervals is denoted by IR and the set of all square interval matrices of order n is denoted by IR n×n . When we need (in a proof) open intervals we write them with brackets, i.e. (a, a).
The magnitude is defined by mag(a) = max(|a|, |a|) which is sometimes confused with the absolute value |a| = {|a|, a ∈ a}. The width of an interval a is defined by w(a) = a − a. All these notions can be intuitively defined for vectors, we just use them component-wise. We will also use the interval vector Euclidean norm The relation a ≤ b holds when a ≤ b (similarly for <). When we compare two interval structures, the relation is applied component-wise. In the following text, E will denote a matrix consisting of ones of a corresponding size. The identity matrix of a corresponding size will be denoted I with e i denoting its i-th column. By A + we denote the Moore-Penrose pseudoinverse matrix to A and by A −T we denote the inverse matrix to A T . Spectral radius of A is denoted (A). Now, we define the main notion of this work.
Definition 2.1 (Interval determinant). Let A be a square interval matrix, then its interval determinant is defined by Computing the exact bounds, i.e., hull, of det(A) is a hard problem. That is why, we are usually satisfied with an enclosure of the interval determinant. Of course, the tighter the better.
Definition 2.2 (Enclosure of interval determinant). Let A be a square interval matrix, then an interval enclosure of its determinant is defined as any d ∈ IR such that det(A) ⊆ d.
3. What was known before. As it was said in the introduction, to the best knowledge of ours, there are only a few theoretical results regarding interval determinants. Some of them can be found in e.g., [20,31,36]. From linearity of a determinant with respect to matrix coefficients we immediately get the fact that the exact bounds on an interval determinant can be computed as minimum and maximum determinant of all 2 n 2 possible "edge" matrices of A.
Proposition 3.1. For a given square interval matrix A the interval determinant can be obtained as Regarding complexity of determinant computation we have the following theorem [20,36].
where A is rational nonnegative is NP-hard.

Approximations.
In the end of the previous section we saw that the problem of computing the exact bounds of an interval determinant is generally an NP-hard problem. One can at least hope for having some approximation algorithms. Unfortunately, we prove that this is not the case, neither for relative nor for absolute approximation.
Theorem 4.1 (Relative approximation). Let A be an interval matrix with A c nonnegative positive definite matrix and A ∆ = E. Let ε be arbitrary such that 0 < ε < 1. If there exists a polynomial time algorithm returning [a, a] such that then P = NP.
Proof. From [36] we use the fact that for a rational nonnegative symmetric positive definite matrix A, checking whether the interval matrix We show that if such algorithm existed, it would decide whether a given interval matrix is regular.  Proof. Let matrix A c consist of rational numbers with nominator and denominator representable with k bits (we can take k as the maximum number of bits needed for any nominator or denominator). Then nominators and denominators of coefficients in A c − E and A c + E are also representable using O(k) bits. For each row we can multiply these matrices with product of all denominators from both matrices in the corresponding row. Now, each denominator uses still k bits and each nominator uses O(nk) bits. We obtained a new matrix A . The whole matrix now uses O(n 3 k) bits which is polynomial in n.
We only multiplied by nonzero constants therefore the following property is holds After cancellation the new matrix A has integer bounds. Its determinant must also have integer bounds. Therefore deciding whether A is regular means deciding whether | det(A )| ≥ 1. We can multiply one arbitrary row of A by 2ε and get a new matrix A having det(A ) = 2ε det(A ). Now, we can apply the approximation algorithm and compute absolute approximation [a , a ] of the determinant of A . Let det(A ) ≥ 1. Then det(A ) ≥ 2ε and the lower bound of absolute approximation is On the other hand, if a > 0 then det(A )/2ε = det(A ) ≥ a > 0.
Hence, even det(A ) > 0 and since it is an integer it must be greater or equal to 1. The case of det(A ) ≤ −1 is handled similarly. Therefore, 0 5. Enclosures of determinants -general case.

Gaussian elimination.
To compute a determinant of an interval matrix, we can use the well known Gaussian elimination -after transforming a matrix to the row echelon form an enclosure of determinant is computed as the product of intervals on the main diagonal. For more detailed description of the interval Gaussian elimination see for example [1,16,27]. Gaussian elimination is suitable to be used together with some form of preconditioning (more details will be explained in section 5.6). We would recommend the midpoint inverse version as was done in [40].

Gerschgorin discs.
It is a well known result that a determinant of a real matrix is a product of its eigenvalues. To obtain an enclosure of an interval determinant, any method returning enclosures of eigenvalues of a general interval matrix can be used, e.g., [10,14,19,23]. Here we will employ simple but useful bounds based on the well known Gerschgorin circle theorem. This classical result claims that for a square real n × n matrix A its eigenvalues lie inside the circles in complex plane with centers A ii and radius j =i |A ij |. When A is an interval matrix, to each real matrix A ∈ A there corresponds a set of Gerschgorin discs. Shifting coefficients of A shifts or scales the discs. All discs in all situations are contained inside discs with centers mid(A ii ) and radii rad(A ii ) + j =i mag(A ij ) as depicted in Figure 1.
As in the case of real Gerschgorin discs, it is also well known that in the union of k intersecting circles there somewhere lie k eigenvalues. By intersecting circles we mean that their projection on the horizontal axis is a continuous line. That might complicate the situation a bit. In the intersection of k discs there lie k eigenvalues and their product contributes to the total determinant. That is why, we can deal with each bunch of intersecting discs separately. We compute a verified interval enclosure of a product of k eigenvalues regardless of their position inside this bunch. The computation of the verified enclosure will depend on the number of discs in the bunch (odd/even) and on whether the bunch contains the point 0. In Figure 2 all the possible cases and resulting verified enclosures are depicted. The resulting determinant will be a product of intervals corresponding to all bunches of intersecting discs. The formulas of enclosures are based on the following simple fact. The eigenvalue lying inside an intersection of circles can be real or complex (c + bi). In the second case the conjugate complex number c − bi is also an eigenvalue. Their product b 2 + c 2 can be enclosed from above by a 2 as depicted in Figure  3. The whole reasoning is based on Pythagorean theorem and geometric properties of hypotenuse.

Hadamard inequality.
A simple but rather crude enclosure of interval determinant can be obtained by the well known Hadamard inequality. For an n × n real matrix This inequality is simply transformable for the interval case. Since the inequality holds for every A ∈ A we have It is a fast and simple method. A drawback is that the obtained enclosure is often quite wide. A second problem is that it is impossible to detect the sign of the determinant, which might be sometimes useful.

Cramer's rule.
Our next method is based on Cramer's rule. It exploits methods for computing enclosure of a solution set of a square interval linear system. There are plenty of such algorithms, i.e., [11,25,27,39]. Here we use the method "\" built in Octave interval package. When solving a real system Ax = e 1 using Cramer's rule we obtain where det(A 2:n ) emerges by omitting the first row and column from A and x 1 is the first coefficient of the solution of Ax = e 1 . We have reduced our problem of determinant computation to a problem with lower dimension and we can repeat the same procedure iteratively until the determinant in the numerator is easily computable. For an interval matrix A we actually get where x 1 is an interval enclosure of the first coefficient of the solution of Ax = e 1 , computed by some of the cited methods. Notice that we can use arbitrary index i instead of 1. The method works when all enclosures of x 1 in the recursive calls do not contain 0.

Monotonicity checking. The derivative of a real nonsingular matrix
Provided the interval matrix A is regular and B is an interval enclosure for the set {A −T | A ∈ A}, then 0 ∈ det(A) and the signs of det(A c )B ij give information about monotonicity of the determinant. As long as 0 is not in the interior of B ij , then we can do the following reasoning. If det(A c )B ij is a nonnegative interval, then det(A) is nondecreasing in A ij , and hence its minimal value is attained at In this way, we split the problem of computing det(A) into two sub-problems of computing the lower and upper bounds separately. For each subproblem, we can fix those interval entries of A at the corresponding lower or upper bounds depending on the signs of B ij . This makes the set A smaller in general. We can repeat this process or call another method for the reduced interval matrix.
Notice that there are classes of interval matrices the determinant of which is automatically monotone. They are called inverse stable [33]. Formally, A is inverse stable if |A −1 | > 0 for each A ∈ A. This class also includes interval M-matrices [3], inverse nonnegative [21] or totally positive matrices [6] as particular subclasses that are efficiently recognizable; cf. [13]. 5.6. Preconditioning. In the interval case by preconditioning we mean transforming an interval matrix into a better form as an input for further processing. It is generally done by multiplying an interval matrix A by a real matrix B from left and by a real matrix C from right and we get some new matrix BAC. Regarding determinants, from properties of the interval arithmetics we easily obtain det(B)·det(A)·det(C) ⊆ det(BAC) and we will further use the fact There are many possibilities how to choose the matrices B, C for a square interval matrix. As in [7], we can take the midpoint matrix A c and compute its LU decomposition P A c = LU . When setting B ≈ L −1 P, C = I, we get Another option is using an LDL decomposition. A symmetric positive definite matrix A can be decomposed as A = LDL T , where L is upper triangular with ones on the main diagonal and D being diagonal matrix. By setting B ≈ L −1 , C ≈ B T and obtain det(A) ⊆ det(L −1 AL −T ). In interval linear system solving, there are various preconditioners utilized depending on criteria used [12,17]. The most common choice is taking B ≈ A −1 c , C = I when A c is regular. Then Unlike the previous real matrices, the matrix A −1 c does not have to have its determinant equal to ±1. We need to compute a verified determinant of a real matrix. In [29] there are many variants of algorithms for computation of verified determinants of real matrices. We use the one by Rump [38].
6. Enclosures of determinants -special cases. Even though we are not going to compare all of the mentioned methods in this section, for the sake of completeness we will mention some cases of matrices that enable the use of other tools. For special classes of interval matrices we prove new results stating that it is possible to compute exact bounds of their determinants in polynomial time.
6.1. Symmetric matrices. Many problems in practical applications are described using symmetric matrices. We specify what we mean by an interval symmetric matrix by the following definition.

Definition 6.1 (Symmetric interval matrix). For a square interval matrix A we define
Next we define its eigenvalues. For symmetric interval matrices there exist various methods to enclose each i-th set of eigenvalues. A simple enclosure can be obtained by the following theorem in [14,37].
There exist various other approaches for computing enclosures of the eigenvalues, including [18,22]. There are several iterative improvement methods [4,15], too. For the exact minimum and maximum extremal eigenvalues, there is a closed-form expression [9], which is however exponential.
6.2. Symmetric positive definite matrices. Let A S be a symmetric positive definite matrix, that is, every A ∈ A S is positive definite. Checking positive definiteness of a given symmetric interval matrix is NP-hard [20,34], but there are various sufficient conditions known [35].
The matrix with maximum determinant can be found by solving the optimization problem max log det(A) subject to A ∈ A S since log is an increasing function and det(A) is positive on A S . This is a convex optimization problem that is solvable in polynomial time using interior point methods; see Boyd & Vandenberghe [5]. Therefore, we have: Proposition 6.4. The maximum determinant of a symmetric positive definite matrix is computable in polynomial time.

Matrices with
results in an interval matrix the center of which is the identity matrix I. This motivates us to study such matrices more in detail. Suppose that A is such that A c = I. Such matrices have very useful properties. For example, solving interval linear systems is a polynomial problem [32]. Also checking regularity of A can be performed effectively just by verifying ρ(A ∆ ) < 1; see [27]. Proof. We will proceed by mathematical induction. Case n = 1 is trivial.
We will proceed by mathematical induction. Case n = 1 is trivial. For the general case, we express the determinant of A ∈ A as in (5.1) det(A) = det(A 2:n )/x 1 .
By induction, the smallest value of det(A 2:n ) is attained for A 2:n = A 2:n . Since A c = I and A is regular det(A) > 0, det(A 2:n ) > 0, therefore x 1 > 0 and as it is the first coefficient of the solution of Ax = e 1 , its largest value is attained for A = A; see [32]. Therefore A = A simultaneously minimizes the numerator and maximizes the denominator. Example 6.6. If the condition ρ(A ∆ ) < 1 does not hold, then the claim is not true in general. Consider We have (A ∆ ) = 3 and det(A) = −2, however, det(A) = [ −6, 14]. The minimum bound is attained e.g., for the matrix  Computing the maximum determinant of A is a more challenging problem, and it is an open question whether is can be done efficiently in polynomial time. Obviously, the maximum determinant of A is attained for a matrix A ∈ A such that A ii = A ii for each i. Specifying the off-diagonal entries is, however, not so easy.

Tridiagonal H-matrices. Consider an interval tridiagonal matrix
Suppose that it is an interval H-matrix, which means that each matrix A ∈ A is an H-matrix. Interval H-matrices are easily recognizable, see, e.g., Neumaier [26,27].
Without loss of generality let us assume that the diagonal is positive, that is, a i > 0 for all i = 1, . . . , n. Otherwise, we could multiply the corresponding rows by −1. Recall that the determinant D n of a real tridiagonal matrix can be computed by a recursive formula as follows D n = a n D n−1 − b n c n D n−2 . Complexity of the determinant computation for general tridiagonal matrices remains an open problem, similarly as solving an interval system with tridiagonal matrix [20]. Nevertheless, not all problems regarding tridiagonal matrices are open or hard, e.g., checking whether a tridiagonal matrix is regular can be done in polynomial time [2].
7. Comparison of methods. In this section the described methods are compared. All these methods were implemented for Octave and its interval package by Oliver Heimlich [8]. This package also contains function det, which computes an enclosure of the determinant of an interval matrix by LU decomposition, which is basically the same as the already described Gaussian elimination method and that is why we do not explicitly compare the methods against this function. All tests were run on an 8-CPU machine Intel(R) Core(TM) i7-4790K, 4.00GHz. Let us start with general matrices first. 7.1. General case. For general matrices the following methods are compared: • GE -interval Gaussian elimination • HAD -interval Hadamard inequality • GERSCH -interval Gerschgorin circles • CRAM -our method based on Cramer's rule The suffix "inv" is added when the preconditioning with midpoint inverse was applied and "lu" is added when the preconditioning based on LU decomposition was used. We use the string HULL to denote the exact interval determinant.
Example 7.1. To obtain a general idea how the methods work, we can use the following example. Let us take the midpoint matrix A c and inflate it into an interval matrix using two fixed radii of intervals - The resulting enclosures of the interval determinant by all methods are shown in Table 1 Based on this example it is not worth to test all methods, because some of them do not work well in comparison to others or do not work well without preconditioning. That is why we later test only -GEinv, CRAMinv, HADinv and GERSCHinv.
We can perceive the method GEinv used in [40] as the "state-of-the-art" method. Therefore, every other method will be compared to it. Primarily, for a given matrix A and a method() we compute the ratio of widths of interval enclosures of det(A) computed by both methods as rat(A) = w(method(A)) w(GEinv(A)) .
We test all methods for sizes n = 5, 10, 15, 20, . . . , 50 and random interval square matrices with given fixed radii of intervals (10 −3 or 10 −5 ). We test on 100 matrices for each size. For each size and method average ratio of computed enclosures, average computation time and its variance is computed. It can happen that an enclosure returned by a method is infinite. Such case is omitted from the computation of average or variance.
The remaining part to be described is generation of random matrices. First, a random midpoint matrix with coefficients uniformly within bounds [−1, 1] is generated. Then, it is inflated into an interval matrix with intervals having their radius equal to 10 −3 or 10 −5 respectively.
Let us begin with the average ratios of widths. They are presented Table 2. When the ratio is a number less then 1000, it is displayed rounded to 2 digits. When it is greater, only the approximation 10 x is displayed. Computation times are displayed in Table 3. For each size of matrix the average computation time is displayed; the numbers in brackets are standard deviations. To more clearly see the difference in computation time between the two most efficient methods GEinv and CRAMinv see Figure 4.

Symmetric matrices.
We repeat the same test procedure with the best methods for interval symmetric matrices. Since these matrices have their eigenvalues real we can add the methods using real bounds on real eigenvalues. Symmetric matrices are generated in a similar way as before, only they are shaped to be symmetric. We compare the preconditioned methods with midpoint inverse GEinv, GERSCHinv, HADinv and CRAMinv. We add one new method EIG based on computation of enclosures of eigenvalues using Theorem 6.3. The method GEinv stays the reference method, i.e, we compare all methods with respect to this method.
The enclosures widths for symmetric matrices are displayed in Table 4. We can see that as in the general case CRAMERinv does slightly worse than GEinv. Another thing we can see is that EIG is worse than both CRAMERinv and GEinv.   Table 4 Ratios of enclosures for symmetric matrices with radii 10 −5 and 10 −3 .
The computation times are displayed in Table 5. We can see that EIG shows low computational demands compared to the other methods. One can argue that we can use filtering methods to get even tighter enclosures of eigenvalues. However, they work well in specific cases [15] and the filtering is much more time consuming.

Conclusion.
In the paper we showed that, unfortunately, even approximation of exact bounds of an interval determinant is NP-hard problem (for both relative and absolute approximation). On the other hand, we proved that there are some special types of matrices where interval determinant can be computed in polynomial time -symmetric positive definite, certain matrices with A c = I or tridiagonal H-matrices.  Table 5 Times of computation (in seconds) for symmetric matrices with radii 10 −5 .
We discussed four methods GE -the "state-of-the-art" Gaussian elimination, GERSCH -our generalized Gerschgorin circles for interval matrices, HAD -our generalized Hadamard inequality for interval matrices and CRAM -our designed method based on Cramer's rule. We introduced a method that can possibly improve an enclosure based on monotonicity checking. All methods combined with preconditioning were tested on random matrices of various sizes. For interval matrices with radii less than 10 −3 the methods GEinv and CRAMinv return similar results. The larger the intervals the slightly worse CRAMinv becomes. However, its computation time is much more convenient (it is possible to compute a determinant of an interval matrix of order 50 by CRAMinv at the same cost as an interval matrix of order 20 by GEinv). Matrices of order larger than 5 need some form of preconditioning otherwise GE and CRAM return infinite intervals. In our test cases the lu preconditioning did not prove to be suitable. The methods HAD and GERSCH always return finite intervals, but these intervals can be huge. Both methods work better with the inv preconditioning. The HADinv returns much tighter intervals than GERSCH, however, it can not distinguish the sign of determinant since the enclosure is symmetric around 0.
The analysed properties of the methods do not change dramatically when dealing with symmetric matrices. The newly added method EIG showed constant and not so huge overestimation and much smaller computation times. The possible improvement of EIG enclosures for symmetric matrices (by applying suitable forms of filtering and eigenvalue enclosures) might be matter of further research. There are many more options for future research -studying various matrix decompositions and preconditioners or studying other special classes of matrices.