Single Precision Calculation of Iterative Refinement of Pairs
of a Real Symmetric-Definite Generalized Eigenproblemby Using a Filter Composed of a Single Resolvent

6. Relative Residual of Eigenpair• We use relative residual Θ given in (8) to evaluate thequality of the approximate pair (λ, v).
|| A v − λB v ||Θ ≡.
(8)|| λB v ||If this number is small, the approximation of the eigenpair is good.

HIROSHI MURAKAMI (Tokyo Metropolitan University)• We solve approximations of those pairs of real sym-def GEVP
whose eigenvalues are in a specified interval by using a filter.
• The filter is constructed by resolvents. We assume some direct
method is used to solve a large simultaneous linear equations which
gives an action of a resolvent.

sin φ ≤ Θ .

Figure 2: Conceptual shape of transfer-func magnitude |g(t)|(for lower-exterior eigenvalues)
9173. Designs of Present Filters• The filter is specified by three parameters ( n, µ, gs ).
– n : degree of Chebyshev polynomial.

FE partitionings for cubes is common : (N1, N2, N3) = (50, 60, 70).
A and B have size N =210, 000 and lower-bandwidth wL=3, 051.

– µ : normal coordinate of the edge of the stop-band.

• Calculations were made using IEEE 754 Single Precision(FP32, precision 7.2 digits).

– gs : bound of transfer-func magnitude in stop-band.

– a real polynomial P of a single resolvent whose shift ρ is real.
!
F ≡ P R(ρ) .

– The range of transfer-func in the pass-band is [ gp, 1 ].

– mixed type (a real polynomial P of the real part of a linear combination ofresolvents).
!
PF ≡ P c∞I + Re kj=1 cj R(ρj ) .

Table 2: EX-1: number of iterations vs. number of approximate pairs
and the max of relative residuals(the number of true pairs is 402)

• For each filtering, ratios of content rates of unrequired
eigenvectors to content rates of required eigenvectors
are reduced by factors of gs / gp or less.

• Fortran code written for double precision (FP64) was
just converted to single precision (FP32) as :REAL(KIND=8)→ REAL(KIND=4)COMPLEX(KIND=8) → COMPLEX(KIND=4)• With lower precision, there is less margin for accuracy.
But we tried, because in recent years, power saving bylower precision calculation has attracted attention.

10225

7. Experiments of Iterative Refinements (in S-P)

– the real part of a linear combination of resolvents whose shifts ρj are complex.
PF ≡ c∞I + Re nj=1 cj R(ρj ) .

– a real poly P of the imag-part of a single resolvent whose shift ρ′ is imaginary.
!
F ≡ P Im R(ρ′) .

(9)deg n = 4, gs = 10−3
IT # pairs Θmax
1(5)2.3E+002 (394) 2.3E-01
3402
1.6E-024402
9.6E-045402
3.8E-046402
3.8E-04

deg n = 4, gs = 10−4
IT # pairs Θmax
1(82) 5.0E-01
2402
6.8E-023402
2.3E-034402
4.1E-04
5402
3.9E-046402
4.1E-04

deg n = 4, gs = 10−5
IT # pairs Θmax
1 (139) 1.6E-01
2402
2.7E-023402
2.7E-034402
8.7E-04
5402
4.5E-046402
4.5E-04

deg n = 6, gs = 10−5
IT # pairs Θmax
1 (222) 2.3E-01
24021.1E-02
3402
3.8E-044402
3.8E-04

deg n = 8, gs = 10−5
IT # pairs Θmax
1 (264) 1.8E-01
2402
3.6E-033402
4.0E-044402
3.9E-04

deg n = 10, gs = 10−5
IT # pairs Θmax
1 (287) 1.7E-01
2402
2.2E-033402
3.9E-044402
4.0E-04

1826

1

7.1. Designs of filters used for experiments• Therefore, we desire to reduce the number of resolvents.

• From the specified parameter set ( n, µ, gs ),the shift ρ, the coefficient of resolvent γ and also gp arecalculated by (3).

! 1−1 1 ,coshσ←µsinhgs2nρ ← a − (b − a) σ,(3)γ←(b−a)(σ+µ),p
 gp ← gs cosh 2n sinh−1 (µ − 1)/(1 + σ) .

• The filter F is given by (4).

F ≡ gs Tn (2γ R(ρ) − I) .

(4)

• For interior pairs, the filter is a deg n Chebyshev-polyof an imag-part of a single resolvent whose shift is imag.

– Degree n is 4, and values of gs are 10−3, 10−4 and 10−5.
– Value of gs is 10−5, and degrees n are 6, 8 and 10.

– a polynomial of a single resolvent whose shift is real.
– a polynomial of the imag-part of a single resolvent
whose shift is imaginary.
But, transfer properties of these filters cannot be madevery good even their shifts and polynomials are tuned.
For example, it is not possible to make the steep discriminating property or the higher uniformity of transferrate magnitude in the pass-band.
(The filter whose shift is imaginary has better properties
compared with the one whose shift is real.)

LOG10 THETA-4-4-420• When the shift is imaginary,
the interval [ a, b ] for eigenvalue can be placed anywhere.

• But, when the filtering is just repeated, magnitudes of
eigenvectors are enhanced or depressed by powers of
their transfer-rate magnitudes.
⇒ Eigenvectors whose transfer-rate magnitudes are relatively small lose relative accuracy or may vanish.
5• From the specified parameter set ( n, µ, gs ),the shift ρ′, the coefficient of the resolvent γ ′ and also
gp are calculated by (5).

!11
−1
σ ← µ sinh 2n cosh gs ,!
 ρ′ ← a+b + b−a σ √−1,22(5)!
2µ +σ 2b−a′


γ ← 2

σ ,p



−1
(µ2 − 1)/(1 + σ 2) .
gp ← gs cosh 2n sinh

• In each iteration the same filter is used,
and the set of vectors is improved by the combinationof B-orthonormalization and filtering.

gs

gp

gs / gp

gp

gs / gp4

10−3

1.93×10−2

5.19×10−2

7.16×10−2

1.40×10−2

410−4

3.56×10−3

2.81×10−2

1.88×10−2

5.33×10−3

410−5

5.33×10−4

1.88×10−2

3.69×10−3

2.71×10−3

610−5

1.53×10−3

6.54×10−3

1.25×10−2

8.01×10−4810−5

2.55×10−3

3.92×10−3

2.11×10−2

4.73×10−4

1010−5

3.34×10−3

2.99×10−3

2.74×10−2

3.65×10−4– We use 100 times the machine epsilon as the absolutevalue of the threshold.
– Those singular vectors whose singular values are below the threshold are cut (removed).

-4-5100204060
EIGENVALUE

801000deg n = 8, gs = 10−5204060
EIGENVALUE

80100deg n = 10, gs = 10−5

• We try to solve those 801 pairs whose eigenvalues are in
the interior interval [ a, b ] = [ 100, 200 ].
• There are 1, 192 eigenvalues in [ a′, b′ ] = [ 75, 225 ] which is
the union of the pass-band and transition-bands.
• We use 1, 300 for m the number of vectors for filtering,which is more than 1, 192 and to be sufficient.
• Results of experiments are shown (Tab. 3, Fig. 9).

Table 3: EX-2: number of iterations vs. number of approx pairs andthe max of relative residuals(the number of true pairs is 801)

gs=1e-3
gs=1e-4
gs=1e-5-2-3-5-6
00.5

1

1.5

2

2.5

3

T

Figure 4: (for lower-exterior pairs) Transfer-func mag |g(t)| (n = 4)

deg n = 4, gs = 10−3
IT # pairs Θmax
1 (329) 1.9E-01
2801
5.3E-023801
9.4E-044801
3.1E-05

deg n = 4, gs = 10−4
IT # pairs Θmax
1 (598) 1.8E-01
2801
8.8E-033801
5.6E-054801
2.6E-05

deg n = 4, gs = 10−5
IT # pairs Θmax
1 (701) 3.0E-01
2801
2.4E-033801
2.9E-054801
2.6E-05

deg n = 6, gs = 10−5
IT # pairs Θmax
1 (800) 3.2E-01
2801
2.2E-043801
2.7E-054801
2.7E-05

deg n = 8, gs = 10−5
IT # pairs Θmax
1 (825) 3.3E-01
2801
8.3E-053801
3.6E-054801
3.5E-05

deg n = 10, gs = 10−5
IT # pairs Θmax
1 (828) 3.2E-01
2801
6.0E-053801
3.8E-054801
3.8E-0529

IT=1IT=2IT=3
IT=4-1n=10n= 8n= 6n= 4-3-1-1-2IT=1IT=2
IT=3
IT=4-3

IT=1
IT=2
IT=3
IT=4-3-4-4-4-5-5-5

LOG10 | G(T) |-6

the number of vectors m in X (i) and Y (i) is updated.

100-2120

140160180-6200deg n = 4, gs = 10−3-600.5

11.5T

22.50

LOG10 | G(T) |-1180-30-1-1-1-2
IT=1
IT=2
IT=3
IT=4-3-2
IT=1
IT=2
IT=3
IT=4-3

IT=1
IT=2
IT=3
IT=4-3-4-4-4-5-5-5120

140160
EIGENVALUE180deg n = 6, gs = 10−5200-6-6100120140160
EIGENVALUE

180deg n = 8, gs = 10−5

200100120140160
EIGENVALUE180

deg n = 10, gs = 10−530
Figure 9: EX-2: eigenvalue vs. relativeresidual(the num of true pairs is 801.
m=1, 300 vectors are filtered.)• A Chebyshev polynomial is used for the polynomial to
make the filter design simple and also the filtering easy
to implement by using the three-term recurrence.

-4

0.5

1

1.5

2

2.5

3

Figure 6: (for interior pairs) Transfer-func mag |g(t)| (n = 4, right-half)• But a filter with such easy and simple construction cannot have a good transfer properties.
⇒ Approximate pairs are not highly accurate or missed
especially with low precision calculation.

23• Expansion basis inside an element: tri-linear functions
in each direction.

31• However, it has been found that even the filter’s properties are not very good, repeating the combination oforthonormalization and filtering several times for a setof vectors improves the approximate pairs.

n=10n= 8
n= 6
n= 4

0-1

LOG10 | G(T) |(N1 ≤ N2 ≤ N3 ) .

200deg n = 4, gs = 10−5

• We used the filter which consists of a single resolvent tolessen both amounts of computation to factor matrices
and especially storage to hold the matrix factors.

-215• The number of exact eigenvalues in any interval can beobtained by counting.

140
160
EIGENVALUE• By using a filter, we approximate pairs of a real-sym-def
GEVP whose eigenvalues are in a specified interval.

T• The exact eigenvalues of this GEVP are given by simple
formulas.

120

8. Conclusion
gs=1e-3
gs=1e-4
gs=1e-50• To this GEVP, the filter diagonalization method isapplied to solve pairs with λ ∈ [ a, b ].

1000100-6Lower bandwidth of A and B : wL = 1 + N1 + N1 N2.

200

0-6-5

Figure 3: Concept of FE partitioning of a cube. Case (N1, N2, N3) = (3, 5, 6) .

180

3Figure 5: (for lower-exterior pairs) Transfer-func mag |g(t)| (gs = 10−5)• An EVP of 3-D Laplacian (7) is considered for a cubewith a side length π with zero-Dirichlet boundary.

• Sides of the cube are equi-divided into N1 + 1, N2 + 1,N3 + 1 sub-intervals, to obtain finite elements (Fig. 3).

140160
EIGENVALUE

deg n = 4, gs = 10−4

LOG10 THETA

LOG10 THETA-522−∆ Ψ(x, y, z) = λ Ψ(x, y, z) .
(7)
By FEM discretization, we obtain a real symmetricdefinite GEVP A v = λB v .

120-3-4

4. We construct Z the basis of approximate invariant-subspacerequired from linear combinations of columns of Y (IT).

-6100

EIGENVALUE

-2

160-1

3. for i = 1, 2, . . . , IT doX (i) ⇐ B-orthonormalization of Y (i−1)； Y (i) ⇐ F X (i)；enddo

• Size of both A and B : N = N1 N2 N30-2
LOG10 THETA

2. Y (0) ⇐ a set of m random vectors.

-2
LOG10 THETA

0eigenvalue may be placed anywhere, and the filter is (2).
! ′′(2)F ≡ gs Tn 2γ Im R(ρ ) − I .

Figure 1: Eigenvalue λ and normalized coordinate t(for lower-exterior eigenvalues)
8-50deg n = 6, gs = 10−5-1

1. Factor shift matrix C(ρ) = A−ρB to prepare the filter F.

5. GEVP for Experiments

7-2-421

2. Simple Type Filters for Present Experiments

Tn(x) : deg n Chebyshev poly.
I: identity.
Im : operator to take imag-part.
γ and γ ′: real consts.
gs : max value of transfer-func magnitude in stop-band.

-2280• For B-orthonormalization, B-SVD with threshold is used.

14• If an imaginary shift ρ′ is used, the interval [ a, b ] for-1-480

100

IT=1
IT=2
IT=3
IT=40-34060
EIGENVALUE80

1-32060

deg n = 4, gs = 10−5-3-540EIGENVALUE-1020-46(1)

0

IT=1
IT=2
IT=3
IT=40-2

for interior pairsn• Starting from a set of random vectors,
we iteratively improve the set of vectors spanning theapproximate invariant subspace.

• To B-orthonormalize a set vectors, we use B-SVD, thesingular value decomposition with metric B.

F ≡ gs Tn ( 2γ R(ρ) − I ) .

100

120

5. Approximate pairs of the original GEVP are the Ritz
pairs obtained from Rayleigh-Ritz procedure applied tothe basis Z.

• When interval [ a, b ] for eigenvalue is lower-exterior,the shift ρ may be a real number, and the filter is (1).

80(6)

Note, if the effective rank is found reduced in orthonormalization,

Filter is a polynomial of a single resolvent.
And a Chebyshev polynomial is used for the polynomial.
Transfer properties of this kind of filter is not very good.

60

deg n = 4, gs = 10−4-1

4.1. Iterative refinement by using a filter• Therefore, we make experiments to improve accuracyof approximate pairs by repeating the combination oforthonormalization and filtering several times.

40EIGENVALUE

0This method is similar to simultaneous inverse iteration
which has been used in the structure analysis.
The principle is (simultaneous) orthogonal iterations [1][2][3].

2027
Figure 8: EX-1: eigenvalue vs. relative
residual(the number of true pairs is 402.
m=800 vectors are filtered.)

for lower-ext pairs

13• To prevent loss of accuracy, in every step of iteration
we orthonormalize the set of vectors before it is filtered.

-50

IT=1
IT=2
IT=3
IT=40

Table 1: Properties of designed six filters for lower-exterior pairs andinterior pairs (µ = 1.5) ( gs/gp is the reduction rate per iteration.)

LOG10 | G(T) |• Even properties of the filter is not very good, approximations could be improved by iteration.
(If the same filter is used, matrix factors can be reused.)

-5

100

14. Refinements of Vectors by Using a Filter• If we use 2 to 4 resolvents, properties of the filter canbe improved very much. But, both amounts of computation to factor matrices and especially storage to holdtheir factors are proportional to the number of resolvents.

80

7.3. (EX-2): Solution of interior eigenpairs

1244060
EIGENVALUE

deg n = 4, gs = 10−3

19

• The filter F is given by (6).
!
F ≡ gs Tn 2γ ′ Im R(ρ′) − I .

-2-3

3.2. Filter with an imaginary shift• Two kinds of filters which consist of a single resolvent:-2-30• For both types of filters, the values gp and gs / gp areshown (Tab. 1) for six filter designs.

113-1-3-5• Properties of the filter is specified by a set of parameters( n, µ, gs ) and we set µ = 1.5 fixed.
• For both types of filters for lower-exterior pairs and interior pairs, we prepared six filter designs.

-2

IT=1
IT=2
IT=3
IT=4
IT=5
IT=60-1

LOG10 THETA• In this study, we assume : The simultaneous linear equations of a large size is solved by some direct method.
⇒ Both amounts of computation to factor matrix andespecially storage to hold the factors tend to be thebottle neck.

• When the interval [ a, b ] for eigenvalue is lower-exterior,and a ≤ λmin.

1IT=1
IT=2
IT=3
IT=4
IT=5
IT=60-1LOG10 THETA• For a vector x, the action of a resolvent y ← R(ρ) x is tosolve a simultaneous linear equations C(ρ) y = B x for y,here C(ρ) ≡ A − ρB is the shifted matrix with a shift ρ.
(The dominant part of the filtering calculation.)• For lower-exterior pairs, the filter is a deg n Chebyshev
polynomial of a single resolvent whose shift is real.

0LOG10 THETA

3.1. Filter with a real shift

1

IT=1
IT=2
IT=3
IT=4
IT=5
IT=6

LOG10 THETA• The filter consists of resolvents R(ρj ) ≡ (A − ρj B)−1B.
• Types of filters F by construction are:

• Results of experiments are shown (Tab. 2, Fig. 8).

• When φ is the angle between vectors Av and λBv, then

LOG10 THETA• Those pairs of a real sym-def GEVP Av = λBv whose
eigenvalues are in the interval [a, b] are solved approximately by using a filter.

• We use 800 for m the number of vectors for filtering,which is more than 764 and to be sufficient.

LOG10 THETA

1. Introduction• There are 764 eigenvalues in [ a, b′ ] = [ 0, 150 ] which is theunion of the pass-band and the transition-band.

LOG10 THETA

1• We try to solve those 402 pairs whose eigenvalues are in
the lower-exterior interval [ a, b ] = [ 0, 100 ].

– This value does not depend on the normalization of
the vector v. It does not depend on the simultaneous
scalings of both matrices A and B also.
– For the vector’s norm || · ||, we used 2-norm.

• In this study, the filter we used is a polynomial of a single resolvent
so to reduce both costs to factor the matrix and to hold the factors,
but properties of such kind of filter are not very good and residuals
of approximate pairs will not be very small.
• Even filter’s properties are not very good, the set of vectors spanning the approximate invariant-subspace is improved well by a few
iterations of the combination of orthonormalization and filtering.

7.2 (EX-1): Solution of lower-exterior eigenpairs• Even in single precision calculations, approximate pairsof an eigenproblem are improved to have small relative
residuals after several iterations.

-2

-3

-4

REFERENCES-5[1] Rutishauser, H.: Computational aspects of F. L. Bauer’s simultaneous iteration
method, Numer. Math., Vol.13, No.1, pp.4–13 (1969).

-6[2] Rutishauser, H.: Simultaneous iteration method for symmetric matrices, Numer.
Math., Vol.16, pp.205–223 (1970).

00.5

1

1.5
T

2

2.5

3

Figure 7: (for interior pairs) Transfer-func mag |g(t)| (gs = 10−5, righthalf)
24[3] Golub, G.H. and Van Loan, C.F.: Matrix Computations, 4th Ed., The John Hopkins Univ. Press (2013). (§8.2.4:‘Orthogonal Iteration’).
32200