# An orthogonal high relative accuracy algorithm for the symmetric eigenvalue problem

Ηλεκτρονική - Συσκευές

13 Οκτ 2013 (πριν από 4 χρόνια και 9 μήνες)

81 εμφανίσεις

An orthogonal high relative accuracy
algorithm for the symmetric
eigenvalue problem
Froil¶an Mart¶³nez-Dopico,Julio Moro &
Juan M.Molera
Departamento de Matem¶aticas
Spain
Hagen,July 2000
1
History
(Brief) History of Eigensolvers for Dense Symmetric
Matrices:
1.
Classical methods (QR,Divide-and-conquer,...) do
not give high relative accuracy
2.
Jacobi method works well for positive de¯nite ma-
trices (Demmel&Veseli¶c)
3.
J-Orthogonal method (Slapni·car&Veseli¶c) gives high
relative accuracy for general symmetric matrices
2
Outline of the talk
IDEA:To get the eigenvalues of a general sym-
metric matrix,get the singular values ¯rst,and
then,put the signs
[¤;W] = SSVD(U;§;V )
SUMMARY:
1.
Theory
2.
Algorithm
3.
Error analysis
4.
Numerical experiments
3
Eigenvalues and eigenvectors from SVD
Let A 2 R
(n£n)
be symmetric and invertible,with
²
A = U §V
T
§ = diag[¾
1

2
;:::;¾
n
]
We order:¾
1
>:::> ¾
i
>:::> ¾
r
;r · n,and,
U = [U
1
:::U
i
:::U
r
] V = [V
1
:::V
i
:::V
r
]
U
i
;V
i
2 R
n£n
i
Cluster:f¾
i
;U
i
;V
i
g
²
If all the s-values are di®erent,we knowthat
the e-values are
¸
i
= (v
T
i
u
i

i
and the e-vectors
v
i
or u
i
Av
i
= ¸
i
v
i
4
Eigenvalues and eigenvectors from SVD
How to get the eigensystem from U;§ and V
when there are clusters?
V
T
A
2
V = U
T
A
2
U = §
2
From here we get
V
T
i
U
j
= 0 i 6= j
Now
V
T
(AV ) = V
T
(U §) =
=
2
4
V
T
1
U
1
.
.
.
V
T
r
U
r
3
5
2
4
¾
1
I
n
1
.
.
.
¾
r
I
n
r
3
5
=
=
2
4
¾
1
¢
1
.
.
.
¾
r
¢
r
3
5
¢
i
= V
T
i
U
i
2 R
n
i
£n
i
are orthogonal and symmetric
¢
i
= W
i
J
i
W
T
i
= [W
+
i
W
¡
i
]J
i
[W
+
i
W
¡
i
]
T
J
i
=
"
+I
m
+
i
0
0 ¡I
m
¡
i
#
5
Eigenvalues and eigenvectors from SVD
V
T
A V =
2
4
W
1

1
J
1
)W
T
1
.
.
.
W
r

r
J
r
)W
T
r
3
5
=
= W
2
6
6
6
6
6
4
¾
1
I
m
+
1
¡¾
1
I
m
¡
1
.
.
.
¾
r
I
m
+
r
¡¾
r
I
m
¡
r
3
7
7
7
7
7
5
W
T
with W =
2
4
W
1
.
.
.
W
r
3
5
and m
+
i
+m
¡
i
= n
i
,m
+
i
¡m
¡
i
= trace(¢
i
);
The (multiplicities of the) e-values corresponding to
each s-value cluster ¾
i
are:
¸
i
= §¾
i
m
§
i
=
n
i
§trace(¢
i
)
2
and e-vectors
V
i
W
§
i
6
Algorithm SSVD (Signed SVD)
Algorithm:
input:A = A
T
2 R
(n£n)
;9A
¡1
Step 0:
Get a rank revealing decomposition of
A,A = XDX
T
Step 1:
Get SVD of A,A = U§V
T
7
Algorithm SSVD (Signed SVD)
Step 2:
for k=1:n
2.a Decide which is the i-th cluster
¾
i
= ¾
k¡n
i
+1
¼:::¼ ¾
k
Get U
i
and V
i
2.b ¢
i
= V
T
i
U
i
If n
i
= 1 then
2.c.1 ¸
i
= ¢
i
¾
i
2.d.1 e-vector q
i
= v
i
else
2.c.2 t
i
= trace¢
i
,
m
¡
i
=
n
i
¡trace¢
i
2
for j = 1:m
¡
i
¸
j
= ¡¾
i
for j = m
¡
i
:n
i
¸
j
= ¾
i
2.d.2 Diagonalize ¢
i
¢
i
= [W
+
i
W
¡
i
]J
i
[W
+
i
W
¡
i
]
T
e-vectors V
i
W
§
i
endif
end
output:for i = 1:r A(V
i
W
§
i
) = §¾
i
(V
i
W
§
i
)
8
Error analysis
Step 1:Get SVD of A = XDX
T
with high relative
accuracy
[
b
U;
b
§;
b
V ] = SGEPVS
[1]
(X;D) (1)
Algorithm:SGEPVS
Given A = A
T
= XDX
T
2 R
(n£n)
1.
Perform QR factorization of XD = QRP;
A = QRPX
T
2.
Multiply to get W = RPX
T
;A = QW
3.
Compute SVD of W =
¹
U§V
T
using one-sided
Jacobi;A = Q
¹
U§V
T
4.
Multiply U = Q
¹
U;A = U§V
T
Theorem [1]
For the algorithmin (1) there exist U
0
and V
0
orthog-
onal matrices such that
1.
k
b
U ¡U
0
k
2
= O(²);k
b
V ¡V
0
k
2
= O(²)
2.
XDX
T
= (I +E
L
)U
0
b
§V
0
(I +E
R
)
with kE
L
k
2
= O(² ·
2
(X)) and kE
R
k
2
= O(² ·
2
(X) ·
2
(R
0
))
and R
0
= D
¡1
R is as well conditioned as possible,with
D diagonal
[1] Demmel,Gu,Eisenstat,Slapni·car,Veseli¶c,Drma·c,LAA 299
(1999)21-80
9
Error analysis
Theorem [1-4]
Let A = XDX
T
= U§V
T
and
b
U;
b
§;
b
V the SVD ob-
tained from (1),then
1.
¯
¯
¯
¯
¾
i
¡b¾
i

i
¯
¯
¯
¯
· O(² ·
2
(X) ·
2
(R
0
))
2.
q
ksin£(U
i
;
b
U
i
)k
2
F
+k sin£(V
i
;
b
V
i
)k
2
F
·
·
p
n
i
O(² ·
2
(X) ·
2
(R
0
))
"
1 +
1
relgap(
b
§
i

¹
i
)
#
3.
9P 2 C
n
i
£n
i
unitary,such that
q
kU
i
P ¡
b
U
i
k
2
F
+kV
i
P ¡
b
V
i
k
2
F
·
·
p
n
i
O(² ·
2
(X) ·
2
(R
0
))
"
1 +
1
relgap
b
(
b
§
i

¹
i
)
#
where
1.
relgap(
b
§
i

¹
i
) = min
¹2§
¹
i
;b¹2
b
§
i
j¹ ¡ b¹j

¹
i
=
[
j6=i
§
j
2.
relgap
b
(
b
§
i

¹
i
) = min
½
relgap(
b
§
i

¹
i
);min
¹2§
i
;b¹2
b
§
i
¹ + b¹
¹
¾
[2] Eisenstat,Ipsen,SIAM JNA 6
(1995)32
[3] Li,SIAM JMA 2
(1999)471
[4] Dopico,BIT 40
(2000)395
10
Error analysis
Step 2:Get the e-values and e-vectors
[
b
¤;
c
W] = SSVD(
b
U;
b
§;
b
V ) (2)
2.a Decide which is the i-th cluster:
fb¾
i
;
b
U
i
;
b
V
i
g
If
¯
¯
¯
¯

i
¡b¾
i+1

i
¯
¯
¯
¯
· O(² ·
2
(X) ·
2
(R
0
)) then b¾
i
¼ b¾
i+1
2.b-c t
i
= trace(
b
V
T
i
b
U
i
),m
¡
i
=
n
i
¡t
i
2
Theorem
¯
¯
¯
fl(trace(fl(
b
V
T
i
b
U
i
)) ¡trace(V
T
i
U
i
)
¯
¯
¯
·
·
p
n
i
O(² ·
2
(X) ·
2
(R
0
))
"
1 +
1
relgap
b
(
b
§
i

¹
i
)
#
²
It is enough to have
¯
¯
¯
fl(trace(fl(
b
V
T
i
b
U
i
)) ¡trace(V
T
i
U
i
)
¯
¯
¯
< 1 because
trace(V
T
i
U
i
) = ¡n
i
;¡n
i
+2;:::;n
i
¡4;n
i
¡2;n
i
In this case:
¯
¯
¯
¯
¯
¸
i
¡
b
¸
i
b
¸
i
¯
¯
¯
¯
¯
· O(² ·
2
(X) ·
2
(R
0
))
²
If high relative accuracy is deserved (O(² ·
2
(X) ·
2
(R
0
)) =
10
¡p
) then it is enough to have relgap
b
& 10
¡p+1
to
get high relative accuracy in the eigenvalues
11
Error analysis
Step 2.d Diagonalize
b
¢
i
=
b
V
T
i
b
U
i
[
c
W
§
i
;
b
J
i
] = SSYEV
³
sym(
b
V
T
i
b
U
i
)
´
c
W
§
i
=
b
V
i
c
W
§
i
Theorem
°
°
°
fl(
b
V
i
c
W
§
i
) ¡Q
§
i
°
°
°
F
·
·
p
n
i
O(² ·
2
(X) ·
2
(R
0
))
"
1 +
1
relgap
b
(
b
§
i

¹
i
)
#
where Q
§
i
is a basis of the direct sum of invariant sub-
spaces corresponding to the positive (negative) eigen-
values in the cluster ¾
i
.
12
Numerical Experiments
DOUBLE PRECISION ²
D
= 1:11 ¤ 10
¡16
²
DSYEVJ [J-Orthogonal method,Slapni·car&Veseli¶c]
[G;J;P] = DGJGT(A);PAP
T
= GJG
T
[Q
(D)
J

(D)
J
] = DJGJF(G;J);A = Q
(D)
J
¤
(D)
J
Q
(D)T
J
²
DRRSSV [The algorithm in this talk]
[U;§;V ] = DGEPVS(A);A = U§V
T
[Q
(D)
S

(D)
S
] = DSSVD(U;§;V );A = Q
(D)
S
¤
(D)
S
Q
(D)T
S
SINGLE PRECISION ²
S
= 5:96 ¤ 10
¡8
²
[Q
(S)
J

(S)
J
] = SSYEVJ(A)
²
[Q
(S)
S

(S)
S
²
[Q
(S)
JAC

(S)
JAC
] = SJAC(A)
Relative errors
Eigenvalues max
i
¯
¯
¯
¯
¯
¯
¸
(D)
Ji
¡¸
(S)
i
¸
(D)
Ji
¯
¯
¯
¯
¯
¯
Eigenvectors kv
(D)
Jr
¡v
(S)
r
k
2
13
Numerical Experiments
Matrices in rank revealing form
A = XDX
T
2 R
(n£n)
n = 50
X;D from LAPACK's:SLATM1,SLATMR,SGESVD
k(X) = 10
2
;10
3
;:::;10
6
s-values aritmetically distr.MODE
X
= 4,
k(D) = 10
2
;10
4
;:::;10
16
s-values randomly distr.MODE
D
= 5,
10 samples for each k(X);k(D) Left Jacobi on W
Figure 1:fg50xl45.eps
14
Numerical Experiments
Matrices in rank revealing form
A = XDX
T
2 R
(n£n)
n = 100
X;D from LAPACK's:SLATM1,SLATMR,SGESVD
k(X) = 10
2
;10
3
;:::;10
6
s-values geometrically distr.MODE
X
= 3,
k(D) = 10
2
;10
4
;:::;10
16
s-values geometrically distr.MODE
D
= 3,
10 samples for each k(X);k(D) Left Jacobi on W
Figure 2:fg100xl33.eps
15
Numerical Experiments
A = DA
0
D 2 R
(n£n)
n = 50
A
0
;D from LAPACK's:SLATM1,SLATMR
k(A
0
) ¼ 1 matrix elements randomly distr.MODE
A
0
= 6,
k(D) = 10
2
;10
4
;:::;10
12
s-values aritmetically distr.MODE
D
= 4,
10 samples for each k(D) Left Jacobi on W
Figure 3:fg50al64.eps
16
Numerical Experiments
A = DA
0
D 2 R
(n£n)
n = 100
A
0
;D from LAPACK's:SLATM1,SLATMR
k(A
0
) ¼ 1 matrix elements randomly distr.MODE
A
0
= 6,
k(D) = 10
2
;10
4
;:::;10
10
s-values randomly distr.MODE
D
= 5,
10 samples for each k(D) Left Jacobi on W
Figure 4:fg100al65.eps
17
Numerical Experiments
Matrices with clusters
A = XD
Bl
X
T
2 R
(n£n)
n = 100
D
Bl
= diag[I
N=4
;(1=
p
k(D))I
N=2
;(1=k(D))I
N=4
]
X;D
Bl
from LAPACK's:SLATM1,SLATMR
k(X) = 1 +²
0:3
k(D) = 10
2
;10
4
;:::;10
10
s-values randomly distr.MODE
D
= 5,
10 samples for each k(D) Left Jacobi on W
Figure 5:fg100xbl3l55.eps
18
Numerical Experiments
Matrices with clusters
A = XD
Bl
X
T
2 R
(n£n)
n = 100
D
Bl
= diag[I
N=4
;(1=
p
k(D))I
N=2
;(1=k(D))I
N=4
]
X;D
Bl
from LAPACK's:SLATM1,SLATMR
k(X) = 1 +²
0:5
k(D) = 10
2
;10
4
;:::;10
16
s-values randomly distr.MODE
D
= 5,
10 samples for each k(D) Left Jacobi on W
Figure 6:fg100xbl5l55.eps
19
Numerical Experiments
Matrices in rank revealing form
A = XDX
T
2 R
(n£n)
n = 50
X;D from LAPACK's:SLATM1,SLATMR,SGESVD
k(X) = 10
2
s-values geometrically distr.MODE
X
= 3,
k(D) = 10
2
;10
4
;:::;10
16
s-values geometrically distr.MODE
D
= 3,
20 samples for each k(D) Left Jacobi on W
Figure 7:v50xdx.eps
20
Numerical Experiments
Matrices in rank revealing form
A = XDX
T
2 R
(n£n)
n = 50
X;D from LAPACK's:SLATM1,SLATMR,SGESVD
k(X) = 10
2
s-values geometrically distr.MODE
X
= 3,
k(D) = 10
2
;10
4
;:::;10
16
s-values geometrically distr.MODE
D
= 3,
20 samples for each k(D) Left Jacobi on W
Figure 8:v50xdxm.eps
21
Numerical Experiments
Matrices in rank revealing form
A = XDX
T
2 R
(n£n)
n = 10
X;D from LAPACK's:SLATM1,SLATMR,SGESVD
D = diag[I
8
;²;¡²]
k(X) = 1 +
p
² Tol
Cluster
= n ¤ ²
10 samples Left Jacobi on W
Figure 9:v10rgap.eps
22
Numerical Experiments
Matrices in rank revealing form
A = XDX
T
2 R
(n£n)
n = 10
X;D from LAPACK's:SLATM1,SLATMR,SGESVD
D = diag[I
8
;²;¡²]
k(X) = 1 +
p
² Tol
Cluster
=
p
n ¤ ²
10 samples Left Jacobi on W
Figure 10:v10rgapsq.eps
23
Numerical Experiments
A =
2
6
6
6
6
6
6
4
0 w 0 0 0 0
w 0 1 0 0 0
0 1 0 w 0 0
0 0 w 0 1 0
0 0 0 1 0 w
0 0 0 0 w 0
3
7
7
7
7
7
7
5
,w = 10
¡3
,(Veseli¶c 1996)
Exact e-values ¼ f1;1;¡1;¡1;10
¡9
;¡10
¡9
g
Relative errors
Eigenvalues max
i
¯
¯
¯
¯
¯
¸
(D)
Ji
¡¸
(S)
i
¸
(D)
Ji
¯
¯
¯
¯
¯
=
= 0:887 ¤10
¡7
¡6
(SY EV J)
Eigenvectors kV
(D)
J
¡V
(S)
k
2
=
= 7:7146 ¤10
¡5