Veriﬁed solutions of sparse linear systems
Takeshi Ogita
Division of Mathematical Sciences
Tokyo Woman’s Christian University
Japanjoint work with
Shin’ichi Oishi and Siegfried M. Rump
SCAN2012, Novosibirsk
September 26, 2012
1 / 50Outline
F: a set of ﬁxed precision ﬂoating-point numbers, e.g., IEEE 754 binary64
Let us consider Ax = b where A ∈ Fn×n, b ∈ Fn.
Task
By the use of ﬂoating-point arithmetic (with intervals),
we aim to• prove A is nonsingular, and• compute a forward error bound of an approximate solution x
e of Ax = bs.t.
|(A−1b)i − x
ei| ≤ i for 1 ≤ i ≤ n.

2 / 50Brief assumptions and conditions• The matrix A is large, sparse and moderately ill-conditioned.
• The veriﬁcation process should be as fast as possible.
• Obtained error bounds should be tight (meaningful).
⇓Some information on A−1 is necessary. (To estimate kA−1k is essential.)⇓
One of the Grand Challenges in Interval Analysis[1] A. Neumaier: Grand Challenges and Scientiﬁc Standards in Interval
Analysis, Reliable Computing, 8 (2002), 313–320.
“Apart from a paper by Rump, nothing bas been done on the interval side.”
3 / 50Notation• For x = (x1, x2, . . . , xn)T ∈ Rn
|x| = (|x1|, |x2|, . . . , |xn|)T ∈ Rn
kxk∞ = max |xi|1≤i≤n

• For A = (aij ) ∈ Rm×n
|A| = (|aij |) ∈ Rm×n√√
kAk2 = ρ(AT A) = λmax(AT A)∑
kAk∞ = max|aij |
1≤i≤m

1≤j≤n
4 / 50• For A = (aij ), B = (bij ) ∈ Rm×n
A≤B⇐⇒

aij ≤ bij

for all (i, j)• o: zero vector• e: vector of all ones
• O: matrix of all zeros• I: identity matrix• u: rounding error unit (unit round-oﬀ), u ≈ 10−16 in IEEE 754 binary64• κ(A) = kAk · kA−1k: condition number
5 / 50Diﬃcult points for sparse matrices
For dense linear systems there are several eﬃcient methods for thispurpose (e.g. by Rump (1980), Oishi-Rump (2002), Hansen-Bliek-RohnNing-Kearfott ([Neumaier] 1999)).
• Common basis: use of an approximate full inverse of either A or its
LU factors.
• Cost: comparable with a standard numerical algorithm, Gaussianelimination with partial pivoting.
• Applicability: κ(A) . 1/u ∼ 1016 in binary64.
• Model implementation: verifylss in INTLAB, a Matlab toolbox for
reliable computing.
6 / 50For sparse cases things are much diﬀerent: Still diﬃcult in terms of bothcomputational complexity and memory requirements.
• Diﬃculty: destruction of the sparsity of A if using full inverses.
• Exception: diagonally dominant and M -matrix or alike.

7 / 50A (n = 48)

A −1

0

0

10

10

20

20

30

30

40

40

0

10

20
30
40
nnz(A) = 400

0

10
20
30
40
−1
nnz(A ) = 2304 = 48 2

Figure 1: Destruction of sparsity of A (n = 48).

8 / 50Figure 2: Destruction of sparsity of A (n = 1600).

9 / 50More precisely

Prof. Rump formulated the following challenge:

Derive a veriﬁcation algorithm which computes an inclusion of thesolution of a linear system with a general symmetric sparse matrixof dimension 10000 with condition number 1010 in IEEE 754 double
precision, and which is no more than 10 times slower than the best
numerical algorithm for that problem.
[2] S. M. Rump: Veriﬁcation methods: Rigorous results using ﬂoatingpoint arithmetic, Acta Numerica, 19 (2010), 287–449.
• κ(AT A) = κ(A)2• Treatable range in ﬂ-pt: κ(A) . u−1 ≈ 1016 in binary64• If κ(A) is small, then a super-fast veriﬁcation method for s.p.d. matrices(to be explained) can be used after calculating AT A.
10 / 50

In this talk we aim to do the following things:
1. survey existing veriﬁcation methods for sparse linear systems.
• monotone (including M-matrix) [e.g. heat equation]• H-matrix [e.g. ﬂuid dynamics, electromagnetics]• symmetric and positive deﬁnite [e.g. structure analysis]• general symmetric
2. try to partially solve the problem for general symmetric matrices:• A is large, e.g. n ≥ 10000, and sparse.
√• A is moderately ill-conditioned, e.g. u−1 < κ(A) < u−1.

11 / 50

Basic principles of veriﬁed numerical computations
1. Utilize results by standard (non-interval) numerical algorithms with pureﬂoating-point arithmetic as much as possible.
• Quality of such results are usually good.
• There are many fast and reliable (but not veriﬁed) numerical libraries
such as BLAS/LAPACK and sparse routines.
2. Use interval arithmetic only if absolutely necessary.
• To avoid slowing down computational speed.
• To avoid explosions of interval width.
=⇒

Leave the use of interval arithmetic as late as possible. [Wilkinson]
12 / 50

Current status of fast veriﬁed solutions of linear systems
dense

directsolvergeneralRump (1980), Oishi-Rump (2002)
s.p.d.
Rump (1993), Rump-Ogita (2007)
H-matrix
Ning-Kearfott (1997)
sparsedirectgeneralRump (1994)
solvers.p.d.
Rump (1993), Rump-Ogita (2007)
symmetric
Rump (1995)
strictly
diagonallyanydominant(trivial)(including monotone∗
Ogita-Oishi-Ushiro (2001)
iterative
H-matrix
Ogita-Oishi (2006)
solver)TN matrix
similar to monotone
others–∗) It is not trivial to determine whether a give matrix is monotone.
13 / 50

Dense matrices (for reference)

14 / 50

Veriﬁcation methods for dense matrices:(1) Krawczyk-Rump
Theorem (Rump, 1980)
Let A ∈ Rn×n, R ∈ Rn×n, b ∈ Rn and xe ∈ Rn be given. Let [] ∈ IRn
be closed and bounded with [] 6= ∅. Let int([]) denote the interior of[]. If[y] := R(b − Aex) + (I − RA)[] ⊆ int([]),then A is nonsingular andA−1b ∈ x
e + [y].

The 1st stage of INTLAB’s verifylss for dense linear (interval) systems.

15 / 50

Veriﬁcation methods for dense matrices:(2) Hansen-Bliek-Rohn-Ning-Kearfott

Theorem (Ning-Kearfott, 1997)
Let an H-matrix A ∈ Rn×n and b ∈ Rn be given. Let y, z ∈ Rn be
deﬁned by y := M(A)−1|b| and zi := [M(A)−1]ii. Let p, q ∈ Rn be
deﬁned by pi := [M(A)]ii − zi and qi := yi/zi − |bi|. Then A−1b ∈ [x]wherebi + [−qi, qi][xi] :=.
Aii + [−pi, pi]







• The 2nd stage of verifylss for dense linear (interval) systems.
• The results may be of better quality than those of the Rump’s approach
for ill-conditioned linear systems; normally the quality is similar.
16 / 50

Veriﬁcation methods for dense matrices:(3) Yamamoto
Theorem (Yamamoto, 1984)
Let A, R ∈ Rn×n and b ∈ Rn.
If kI − RAk∞ < 1, then



−1

|A




kR(b − Aex)k∞b−xe| ≤ |R(b − Aex)| +|I − RA|e.
1 − kI − RAk∞



• It is easy to implement the method.
• The results are usually as good as those of the Rump’s approach.

17 / 50

Veriﬁcation methods for dense matrices:(4) Oishi-Rump
Key estimation: kI − RAk∞
1. P A ≈ LU . [ 23 n3 ﬂops]
2. XL ≈ L−1 and XU ≈ U −1. [ 23 n3 ﬂops in total]
3. R := XU XLP (not explicitly compute it).
4. Use a priori error bounds by backward error analysis.
Evaluation in O(n2) ﬂops:
kI − RAk∞ ≤ c1uk|XU |(|XL|(|L|(|U |e)))k∞ + c2u,
where c1, c2 are some computable factors and u is the underﬂow unit.
• The same computational eﬀort for calculating an approximate solution.
18 / 50

Sparse matrices

19 / 50

Veriﬁcation methods for sparse matrices:
strictly diagonally dominant
Suppose A = (aij ) ∈ Rn×n is strictly (row) diagonally dominant.
e := A − D.
Let D := diag(a11, . . . , ann) and A−1
Setting R := D−1 = diag(a11, . . . , a−1nn ) yields

e ∞ < 1,kI − RAk∞ = kI − D−1Ak∞ = kD−1Ak
since

∑

j6=i |aij | < |aii| for all i.

20 / 50

Veriﬁcation methods for sparse matrices:monotone (including M-matrix)
monotone = inverse nonnegativeDeﬁnition 1. [monotone] A matrix A ∈ Rn×n is called monotone if
Av ≥ o for v ∈ Rn implies v ≥ o.
Lemma 2. A is monotone if and only if A is nonsingular with A−1 ≥ O.
Deﬁnition 3. [M-matrix] Let A = (aij ) ∈ Rn×n with aii > 0 andaij ≤ 0 for i 6= j. Then A is called an M-matrix if A is nonsingular andA−1 ≥ O.

21 / 50

Theorem (Ogita-Oishi-Ushiro, 2001)
Let A ∈ Rn×n with A being monotone and b, ye ∈ Rn.
If ke − Aey k∞ < 1, then



−1

kA




ke
y k∞
.
k∞ ≤
1 − ke − Aey k∞• To solve Ay = e, the same solver for solving Ax = b can be applied.
• ke − Aey k∞ < 1 is suited as a stopping criterion for iterative solvers.
• It is not trivial to determine whether A is monotone.

22 / 50

Proof of the theorem for monotone matrices
Since A−1 ≥ O, we havekA−1k∞ = k|A−1|ek∞ = kA−1ek∞≤ kA−1e − yek∞ + key k∞
≤ kA−1k∞ke − Aey k∞ + key k∞.
This yields(1 − ke − Aey k∞)kA−1k∞ ≤ ke
y k∞.
If ke − Aey k∞ < 1, then
−1

kA

ke
y k∞
k∞ ≤
.
1 − ke − Aey k∞23 / 50

Numerical results (1)• A, b: from discretizing 2-D Poisson’s equation by FEM• The problem size n is varied from 10, 000 to 250, 000.
• Solver: MICCG methodkb − Aexk2– stopping criteria:≤ 10−12,kbk2

ke − Aey k∞ ≤ 10−3

24 / 5023y(h : un hosen; T1 = 0)x

: k1 = 1:0; f1 = 0

 = 4: k2 = 0:1; f2 = 0
2

2

1
10

5

: k3 = 1:0; f3 = 20



in Ω
 div{−k · grad(u)} = f{−k · grad(u)} × n = 0
on Γ2 {−k · grad(u)} × n = h(u − T∞) on Γ3

25 / 50

Table 1: Computing time and relative error bound kA−1b− x
ek∞/kA−1bk∞dim(A) (n) approx. solution [s] veriﬁcation [s] rel. error bound
10,000
3.3
1.7
4.1 × 10−10
40,00027.1
10.2
2.5 × 10−9
90.7
32.3
7.1 × 10−9
90,000160,000
216.2
77.0
1.6 × 10−8
250,000
458.5
146.8
3.3 × 10−8
Intel Celeron 566MHz CPU [Computing, Suppl. 15 (2001)]
Veriﬁcation process can be faster than approximation one!

26 / 50

Veriﬁcation methods for sparse matrices:H-matrix
For A = (aij ) ∈ Rn×n, the comparison matrix M(A) = (baij ) of A isdeﬁned as{|aij | (i = j)
b
aij =
.
−|aij | (i 6= j)Deﬁnition 4. [H-matrix] A is called an H-matrix if M(A) is an Mmatrix.
Lemma 5. A is an H-matrix if and only if there exists a vector v > osuch that M(A)v > o.
Lemma 6. If A is an H-matrix, then |A−1| ≤ M(A)−1.
From Lemma 6, it follows thatkA−1k∞ ≤ kM(A)−1k∞.
27 / 50

How to determine whether A is an H-matrix?
Suppose we do not know whether A is an H-matrix.
b := M(A). (At least A
b is an L-matrix for any A with nonzero
Put A
diagonals; baij > 0 for i = j (baii > 0) and b
aij ≤ 0 for i 6= j.)There are some possibilities:
1. Use an approximation ve of the eigenvector corresponding to the minimumb (that is the Perron vector of A
b−1 if A is an H-matrix)
eigenvalue of A
[Rump, 2012], orb = e. [Neumaier, 1999]
2. Use an approximate solution of Av
bv k∞ < 1 implies Ae
bv > o.
⇒ ke − Ae

28 / 50

Veriﬁcation methods for sparse matrices:
another approach
For symmetric A
λi(A): eigenvalues of A
−1

kA

1
k2 =
min |λi(A)|For non-symmetric A√
σi(A): singular values of A ( = λi(AT A))
−1

kA

1
k2 =
min σi(A)29 / 50

Veriﬁcation methods for sparse matrices:
symmetric and positive deﬁnite
Rump’s algorithm
1. Set α := ψu · tr(A). (ψ: computable)
2. Execute a Cholesky factorization for A − 2αI ≈ LLT .
3. If succeeded, then λmin ≥ α.
An a priori error estimate by a backward error analysis:

kLLT − (A − 2αI)k2 ≤ ψu · tr(A − 2αI) ≤ ψu · tr(A) = α=⇒

INTLAB function: verifylss, isspd
30 / 50

Property• Only one ﬂ-pt Cholesky factorization chol(A − 2αI) is nesessary.
(Direct sparse solvers can be used.) Super-fast!
• If chol(A − 2αI) runs to completion, then it is veriﬁed that “A ispositive deﬁnite”. (and λmin(A) ≥ α > 0)
(It is veriﬁed rigorously.)• Even if chol(A − 2αI) failed, it is not veriﬁed that “A is not positivedeﬁnite”.
(A may be positive deﬁnite, although it is unlikely.)

31 / 50

Numerical results (2)Test matrices: University of Florida Sparse Matrix Collection
Computer environment:
CPU: Intel Dual-Core Xeon 2.80GHz × 4 processors
Memory: 32GB
OS: Red Hat Enterprise Linux WS
Software: Matlab Version 7.1.0.183 (R14) Service Pack 3
nameship 003
shipsec1
cfd2
af shell(3,4,7,8)
apache2n bw w/wo RCM time (sec)
121,7283659/3659
260
140,874
5238/5238
538
123,440
2179/4333
127
504,855
2470/4909
633
715,1762993/65837
1176
32 / 50

Veriﬁcation methods for sparse matrices:
general symmetric
Rump’s algorithm
1. Estimate the smallest magnitude eigenvalue (denoted by τe1).
2. Set α := 0.9 · |eτ1|.
3. Execute an LDLT factorization for A − αI ≈ L1D1LT1 .
4. Compute β1 ≥ kL1D1LT1 − (A − αI)k2.
5. Check the inertia of D1.
6. Execute an LDLT factorization for A + αI ≈ L2D2LT2 .
7. Compute β2 ≥ kL2D2LT2 − (A + αI)k2.
8. Check the inertia of D2.
9. Compute a lower bound of min |λi(A)|: σ ≥ α − max{β1, β2}.
=⇒a little unstable: κ(A) ≤ κ(A ± αI)
33 / 50

at least n2 eigenvaluesje1 jat least n1 eigenvalues

je1j

0
2

2

 
1

1

R

Figure 3: Lower bound of the smallest magnitude eigenvalue
A similar approach for bounding eigenvalues can be found in[3] N. Yamamoto: A simple method for error bounds of eigenvalues ofsymmetric matrices, Linear Alg. Appl., 324 (2001), 227–234.

34 / 50

Veriﬁcation methods for sparse matrices:non-symmetric
The following three approaches are known (Rump):
1. B = AT A and apply the super-fast method for s.p.d. matrices.
2. A = LDM T⇒

σ1(A) ≥ σ1(L) · σ1(D) · σ1(M ).

eDeMfT + ∆ (due to rounding errors) and• In practice, A = L
e · σ1(D)e · σ1(Mf) − k∆k2.
σ1(A) ≥ σ1(L)
[
]
T
O A
3. G :=
and apply any method for symmetric matrices to G.
A O• {λi(G), 1 ≤ i ≤ 2n} = {±σj (A), 1 ≤ j ≤ n}. ⇒ κ(G) = κ(A)• For small α, an LDLT factorization for G − αI is a little unstable.
35 / 50

A new approach for sparse matrices

36 / 50

Lower bound of the smallest singular value• Present status: Few methods of obtaining σ ≤ σ1(A) are known
except some methods by Rump based on LDLT factorization.
• Special case: A super-fast veriﬁcation method for SPD matrices by
Rump using Cholesky factorization.
– applicable up to κ(A) ∼ u−1/ψ where ψ := max nnz(L(i, :)) for ai
Cholesky factor L.
• Suboptimal approach: use of AT A or AAT .
– An obvious drawback: it squares the condition number of A, so that√applicable up to κ(A) ∼ 1/ u ∼ 108.

37 / 50

Preliminaries


Theorem 7. [eigenvalue perturbation] Let A and B
symmetric n × n matrices. Then it holds for i = 1, 2, . . . , n


be

real

|λi(A) − λi(B)| ≤ kA − Bk2.







Theorem 8. Let A = AT ∈ Rn×n. For some α ∈ R, supposeA − αI = XDX T
where X is some nonsingular matrix and D ∈ Rn×n. Then the inertia
of D is equivalent to a triplet of the number of eigenvalues of A which
are larger than, smaller than or equal to α.


38 / 50

Theorem 9. [Lehmann bounds] Let A = AT ∈ Rn×n. Let λi, 1 ≤i ≤ n, be eigenvalues of A withλ1 ≤ · · · ≤ λn .
Suppose ν ∈ R satisfies λk < ν ≤ λk+1 for some k. Let X be a real n×k
matrix of full rank. Put A1 = X T X, A2 = X T AX, A3 = X T A2X,
B1 = νA1 − A2 and B2 = ν 2A1 − 2νA2 + A3. Let µj , 1 ≤ j ≤ k, be
generalized eigenvalues of (B1, B2) with
µ1 ≤ · · · ≤ µk .
If B1 is positive definite, then it holds for j = 1, . . . , k that1
λk−j+1 ≥ ν − .
µj39 / 50

Principle of the proposed algorithm
We try to derive a veriﬁcation algorithm which is• fast (comparable with the cost for LDLT factorization)• stable (applicable for cases κ(A) > 1010)
For this purpose,• Find two approximate eigenvalues τek , τek+1 where the gap |eτk+1| − |eτk |is suﬃciently large.
• Use block LDLT factorizations and their a priori error estimates.
• Apply Lehmann bounds for τej , j = 1, . . . , k.
40 / 50

0

j

e1 j

j

ek j

j

ek+1 jR

Figure 4: Distribution of absolute values of eigenvalues.

(2)

at least n− eigenvalues−|τk+1 |
β2
−α

ν2

γmin

(1)



k approximate eigenvaluesτp(1)0τp(k)σσ

at least n+ eigenvalues|τk+1 |

γmax

ν1

β1

R
α

Figure 5: Distribution of eigenvalues around zero.

41 / 50

Rounding error analysis on block LDLT factorizations
Block LDLT factorization: P AP T = LDLT , where


D=




D1
D2

...



L11


 , L =  L.21 L.22
...

 .
.
D`
L`1 L`2 · · ·





L``

Each Di and Lii is a 1 × 1 or 2 × 2 block, with Lii being 1 or the 2 × 2
identity matrix, respectively. The rest of L is partitioned accordingly.

42 / 50

There are several methods with diﬀerent pivoting strategies:• Bunch–Parlett (1971)• Bunch–Kaufman (1977)• others
For the rounding error analysis, we need a backward error bound for thesolution of linear systems involving 2 × 2 pivots.

43 / 50

We assume that the 2 × 2 linear system Ey = f is solved successfully witha computed solution ye satisfying(E + ∆E)ey = f,|∆E| ≤ c|E|(1)

for some constant c > 0.
Under some conditions we can prove that the condition (1) is rigorously
satisﬁed with

(GEPP)
 4γ21
c =γ298 (the explicit inverse without scaling)6
 1γ
6 556 (the explicit inverse with scaling)where γm = mu/(1 − mu) [≈ mu for not so large m].

44 / 50

Theorem 10. [Ogita-Rump] Let A = AT ∈ Fn×n. Let L = (lij ),D = (dij ) and P be a computed block LDLT factors of A. Suppose
the condition (1) is satisfied for some constant c > 0. For 1 ≤ i, j ≤ ndefines(i, j) := |{k ∈ N : 1 ≤ k < min(i, j) and lik dkk ljk 6= 0}|
and denote{

γs(i,j)+1 if s(i, j) 6= 0
αij :=.
0otherwise
Put 2 = max{c, α} where α = maxi,j αij for 1 ≤ i, j ≤ n. Then it
holds that

|P AP T − LDLT | ≤ 2(P |A|P T + |L||D||L|T ).



45 / 50

Proposed algorithm
1. Find two approximate eigenvalues τek , τek+1 where the gap |eτk+1| −|eτk | is suﬃciently large.
2. Take α in (|eτk |, |e
τk+1|).
3. Execute a block LDLT for P (A − αI)P T ≈ L1D1LT1 .
4. Compute β1 ≥ kL1D1LT1 − P (A − αI)P T k2.
5. Check the inertia of D1.
6. Execute a block LDLT for P (A + αI)P T ≈ L2D2LT2 .
7. Compute β2 ≥ kL2D2LT2 − P (A + αI)P T k2.
8. Check the inertia of D2.
9. Apply Lehmann bounds for τej , j = 1, . . . , k.
β1, β2: computed by a priori error estimates

46 / 50

Numerical results (for reference)
We evaluate the performance of the proposed algorithm.
CPU: 2.66 GHz Intel Core 2 Duo, Memory: 8GB
We implement a hybrid algorithm:
Stage-1: Rump’s algorithm
Stage-2: Proposed algorithm

Example: Random sparse symmetric matrices having 5 clusteredeigenvalues in [1, 1.1] and n − 5 eigenvalues in [102, cnd] in magnitude.
|λ(A)| = {1, 1.025, 1.05, 1.075, 1.1, 102, . . . , cnd}=⇒

κ(A) = cnd, symmetric and indeﬁnite.
47 / 50

Example: random sparse, density = 5/n, min |λ(A)| = 1
Computing time (sec) and lower bound of the smallest magnitude
eigenvaluen condStage-1 Stage-2
lower bound---------------------------------------------10,000 1e12
0.86
-0.916
10,000 1e13
failed1.72
0.998
10,000 1e14
failed1.41
0.974
20,000 1e15
failedfailed-20,000 1e12
1.80-0.58920,000 1e13
failed5.40
0.999
20,000 1e14
failed6.54
0.996
20,000 1e15
failedfailed--48 / 50

n cond
Stage-1 Stage-2
lower bound---------------------------------------------50,000 1e12
8.03
-0.38750,000 1e13
failed74.89
0.999
50,000 1e14
failed48.88 *
0.997
50,000 1e15
failedfailed-100,000 1e12
8.81
-0.426
100,000 1e13
failed153.620.999
100,000 1e14
failed
485.85 *0.999
100,000 1e15
failedfailed--49 / 50

Conclusions• There exists eﬃcient (practically useful) veriﬁcation methods for sparse
matrices having special property.
• Veriﬁed numerical computation for general sparse linear systems is still
diﬃcult.
• Nevertheless, (we) never ever give up!

Thanks for your kind attention!

50 / 50