2009 4

(1)

荻田武史

東京女子大学現代教養学部数理科学科

アルゴリズムによる計算科学の融合と発展筑波大学計算科学研究センター

(2)

概要

目的

LU

_分解，

QR

_分解，

Cholesky

_{分解，特異値分解などの行} 列分解のための高精度なアルゴリズムを提案する．

A: n

× n

_実行列．提案方式は

A

_{が非常に悪条件な場合も取り扱える．} 悪条件

: A

_{の条件数が非常に大きい．}

=

⇒

_例えば，

IEEE 754

_{倍精度の範囲で条件数が}

10

100_を超えるような問題も取り扱える．

=

⇒

_{直接的な応用}

:

_{連立一次方程式，行列式．その他，様々な} 分野に応用可能．

(3)

条件数

正方行列

A

_{に対して，}

A

_{の条件数を下記のように定義する：}

κ(A) :=

∥A∥ · ∥A

−1

∥.

=

⇒ κ(A)

_{は問題の難しさの指標となる}

=

⇒

_{例えば，連立一次方程式}

Ax = b

_{を考える（}

x

∗

:= A

−1

b

_）．右辺に摂動：

Ay = b +

∆b

_（

y

∗

:= A

−1

(b + ∆b)

_）．このとき，解の変化量

∆x = y

∗

− x

∗_{は下記をみたす：}

∥∆x∥

∥x

∗

_∥

≤

κ(A)

∥∆b∥

∥b∥

(4)

IEEE 754

規格の倍精度演算

=

⇒ 1

_{回の演算の相対精度}

: u

≈ 1.11 × 10

−16

=

⇒

条件数が

10

16より大きい（

κ(A) > u

−1）場合，通常の倍精度演算によって得られた計算結果は

1

_{桁も合っていない}_{場合がある．}

=

⇒

_{基本精度の限界} _{（基本精度：単精度や倍精度）}

=

⇒

_{このような問題を}

“

_{悪条件問題}

”

_{と呼ぶことにする．}

=

⇒

悪条件問題に対し，

“

結果精度の高い

”

行列分解アルゴリズムを提案する．

(5)

Notation

A = (a

_ij

), C = (c

_ij

)

∈ R

n×n

• |A|

：

A

_{のすべての要素に絶対値を付けた行列，}

|A| = (|a

_ij

|)

_．

• A ≤ C

：すべての

(i, j)

_{に対して，}

a

_ij

≤ c

_ij

(6)

LU

_{分解（部分軸交換付）}

A

∈ R

n×n_，

P A

≈ LU

_：











P

∈ R

n×n

:

_置換行列

L

∈ R

n×n

:

_単位_{下三角行列}

U

∈ R

n×n

:

_{上三角行列} 例）連立一次方程式

Ax = b

を解く．

→

P Ax = P b

→

LU

ex = b

→

Ly = P b

_を解く．

→

U

ex = y

_を解く．

LU

_{分解の精度をどのように推定するか？}

(7)

残差

∥P A − LU∥

で推定？

例えば，

LU

_{分解の残差}

∥P A − LU∥

_{あるいは相対残差は？}

∥P A − LU∥

∥A∥

< ε

=

⇒

_これでは_不十分

(Matlab demo.)

例）標準的なガウスの消去法を用いた場合：

(8)

メインアイディア

正確な

LU

_分解

P A = b

L b

U

_{を考えると，通常，}

A

_{の悪条件性は上三} 角行列

U

b

_{に反映される．}

κ(A)

∼ κ( b

U )

もし，

(1)

_{を満たすような}

U

b

−1_{の良い近似}

X

_U _{が得られれば，}

κ(A)

と比べて，

κ(AX

_U

)

_{は相対的に小さくなる．例えば，}

κ(AX

_U

)

≪

u

−1_．

∥ b

U X

_U

− I∥ < 1

(1)

=

⇒

A

_{に対する前処理となる（}

Crout

_{型でも同様）}

(9)

A

_を

X

_U _{によって前処理した後は}

∥X

L

AX

U

− I∥ < 1

(2)

あるいは

∥AX

U

X

L

− I∥ < 1

(3)

を満たすような

bL

の近似逆行列

X

_L_{を計算するのは難しくない．}

=

⇒

重要かつ難しい点：どのように悪条件な

A

_の前処理

X

_U _{を計算するか？}

(10)

仮定

κ( b

U )

≫ u

−1_，

X := b

b

U

−1_．

X

₁

= ﬂ ( b

X)

_とする（_{基本精度における最良近似}_）．

=

⇒

X

₁

= b

U

−1

+ ∆

_，

|∆

_ij

| ∼ u| b

U

_ij−1

|

_{（丸め誤差）} このような場合でも，

∥ b

U X

₁

− I∥ < 1

_{は必ずしも成立するわけで} はない：

∥ b

U X

₁

− I∥ = ∥ b

U ( b

U

−1

+ ∆)

− I∥ = ∥ b

U ∆

∥

≈ ∥ b

U

∥∥∆∥ ∼ u · κ( b

U )

> 1

=

⇒

U

b

−1_{の近似には，なんらかの高精度な形式が必要}

(11)

前処理としての近似逆行列

丸め誤差によって，

X

₁_は

U

b

_の_逆行列_{としては機能しない．しかし，}

A

_{の前処理としては機能する（}

U

b

_{の条件数とは無関係に）．}

=

⇒

κ( b

U X

₁

) =

O(u)

· κ( b

U ),

u

≈ 10

−16

=

⇒

条件数を小さくできる！

=

⇒

_{これを利用した積型の反復アルゴリズムを開発：}

P AX

_U

= b

L b

U X

₁

X

₂

. . . X

_k

≈ bL ≈ L

（反復停止条件：

∥P AX

U

− L∥ < ε

tol）

(12)

“

_高精度

”

_な逆

LU

_分解

A

∈ R

n×n_，

P A

≈ LU

_：











P

∈ R

n×n

:

_置換行列

L

∈ R

n×n

:

_単位_{下三角行列}

U

∈ R

n×n

:

_{上三角行列} 許容誤差

ε

_tol

< 1

_に対して

∥L

−1

P AU

−1

− I∥ ≤ ε

_tol _あるいは

∥P AU

−1

− L∥ ≤ ε

_tol

(13)

提案方式の方針

悪条件問題を取り扱うためには，なんらかの高精度演算が必要．

=

⇒

_{すべての計算を多倍長精度で計算するのは}_{効率が良くない}_．

=

⇒

内積計算や行列乗算に特定すれば，比較的高速なアルゴリズムが利用可能．

=

⇒

_{できる限り，高速な}_{通常の浮動小数点演算}_{を利用したい．}

(14)

高精度な内積計算

x, y

∈ F

n_{に対し，内積}

x

T

y =

n

X

i=1

x

_i

y

_i の計算は，数値線形代数の基礎なので，たくさんの高精度なアルゴリズムが提案されている．

(15)

高精度な総和計算／内積計算の文献（

一部

）

1965 Møller (BIT)

1970 Nickel (ZAMM) (Kahan-Babuˇska’s algorithm) 1971 Malcolm (Comm. ACM)

1972 Pichat (Numer. Math.)

1973 KieÃlbaszi´nski (Math. Stos.)

1974 Neumaier (ZAMM) (improved Kahan-Babuˇska’s algorithm) 1977 Bohlender (IEEE Trans. Comput.)

1982 Leuprecht, Oberaigner (Computing)

1986 Kulisch, Miranker (SIAM Review, originally 1970’s in Karlsruhe) 1991 Priest (Proc. IEEE Symposium)

1999 Anderson (SIAM J. Sci. Comput.)

2002 Li et al. (ACM Trans. Math. Softw., XBLAS) 2003 Demmel, Hida (SIAM J. Sci. Comput)

(16)

高精度な行列乗算

仮定任意に高精度な内積の計算結果が得られる．

F:

（基本精度における）浮動小数点数の集合

u:

_{浮動小数点演算の相対精度（}

IEEE 754

_倍精度：

u

≈ 10

−16_）できる限り，通常の浮動小数点数／浮動小数点演算を用いたい．

=

⇒

A

_1:p

:=

P

p_i=1

A

_i_，

B

_1:q

:=

P

q_i=1

B

_i _{のように行列を表現す} る．ただし，

A

_i

, B

_i

∈ F

n×n_かつ

|A

i

| ≤ u |A

i−1

| ,

|B

i

| ≤ u |B

i−1

|

のような関係が成り立っているものとする．

(17)

仮定行列乗算の汎用的な関数として下記のようなものが利用可能とする：

C

_1:L

= [A

_1:p

B

_1:q

]

L_K

,

C

_1:L

:=

L

X

i=1

C

_i

, C

_i

∈ F

n×n は，下記をみたすとする（

c

₁

, c

₂

:

O(1)

_の定数）_：

|C

1:ℓ

− A

1:p

B

1:q

| ≤ c

1

u

L

|A

1:p

B

1:q

| + c

2

u

K

|A

1:p

||B

1:q

|.

これは，

A

_1:p

B

_1:q_を

K

_{倍長精度で計算し，その結果を}

L

_{倍長精度に} 丸めたものに対応する（

K

≥ L

_{でないと意味がない）．}

(18)

Rump

_{は，悪条件な行列の逆行列を高精度に計算するアルゴリズム}

を開発した

[unpublished, 1980’s], [JJIAM, to appear]

_．

function R = AccInv(A,m)

% right inverse version

n = dim(A);

I = eye(n);

R = A \ I;

% pure fl-pt (Solve AR = I for R)

for k=2:m

C = AccMM(A,R,1);

% accurate dot product

T = C \ I;

% pure fl-pt (Solve CT = I for T)

R = AccMM(R,T,k+1);

% accurate dot product

end

(19)

アルゴリズム：高精度な逆

LU

_分解

[

_{高精度演算が必要な箇所のみ記載．それ以外は，基本精度}

]

0:

X

_1:0

= I

_，

k = 1

_とする．

1:

B

_k

← [A · X

_1:k₋₁

]

1_k_．

[

_{高精度に計算}

/

_{基本精度に丸める}

]

2:

B

_k_の

LU

_分解（

P

_k

B

_k

≈ L

_k

U

_k_）．

3:

U

_k−1

≈ T

_kを計算．

4:

X

_1:k

← [X

_1:k₋₁

· T

_k

]

k_k_．

[

/

_{高精度に格納}

]

5:

_もし，

∥U

_k

∥∥T

_k

∥ ≤ ε

_tol

· u

−1_ならば

L := L

_k

, X

_U

:= X

_1:k_および

P := P

_k_を返す．

6:

k

← k + 1

_として，

1

_に戻る．

(20)

アルゴリズム：高精度な逆

QR

_分解

[

_{高精度演算が必要な箇所のみ記載．それ以外は，基本精度}

]

0:

X

_1:0

= I

_，

k = 1

_とする．

1:

B

_k

← [A · X

_1:k₋₁

]

1_k_．

[

/

_{基本精度に丸める}

]

2:

B

_k_の

QR

_分解（

B

_k

≈ Q

_k

R

_k_）．

3:

R

−1_k

≈ T

_kを計算．

4:

X

_1:k

← [X

_1:k₋₁

· T

_k

]

k_k_．

[

/

_{高精度に格納}

]

5:

_もし，

∥R

_k

∥∥T

_k

∥ ≤ ε

_tol

· u

−1_ならば

Q := Q

_k

, X

_R

:= X

_1:k_を返す．

6:

k

← k + 1

_として，

1

_に戻る．

(21)

難しい点

•

アルゴリズムは簡潔だが，解析が難しい

•

実際，厳密な解析は不可能（期待値や確率的な話になる）

(22)

¶ ³

Proposition 1.

A

∈ F

n×n_とする．

X

_1:k _を逆

LU

_分解（逆

QR

分解）アルゴリズムによって得られた

n

× n

上三角行列とする．このとき，ある仮定の下で，

k

≥ 1

_{に対して下記が成立する：}

κ(AX

_1:k

) = 1 +

O(

u

k

)

· κ(A)

(4)

µ ´

Remark 1.

A =

P

m_i=1

A

_i

, A

_i

∈ F

n×n_{のような入力でも取り扱え} るようにアルゴリズムを拡張するのは難しくない（もし，高精度な内積計算／行列乗算のアルゴリズムが利用可能なら）．

(23)

∥ b

RX

_1:k

− I∥ = ∥AX

_1:k

− b

Q

T

∥ < 1

_{を達成するためには，下記のよ}

うな

k

_で十分．

k =

»

log[ε

_tol

· κ(A)

−1

]

log u

¼

.

しかしながら，一般的に

κ(A)

_{は事前にはわからないため，}

k

_をあらかじめ決めることはできない．そこで，下記のような停止条件を利用：

∥R

k

∥∥T

k

∥ < ε

tol

· u

−1

(5)

_が成立

=

⇒

∥AX

_1:k

− b

Q

T

∥ < ε

_tol

(24)

提案アルゴリズムの特徴

1.

_{問題の難しさに応じて，必要なだけ反復的に計算精度を増} やすことができる（アダプティブ性）．

2.

高精度演算は，行列乗算（つまり内積計算）のみ．

3.

_{それ以外の処理は，}

BLAS

_，

LAPACK

_{などを利用可能で，す} べて通常の浮動小数点演算で実行（高速性）．

4.

_{高精度な内積計算／行列乗算が高速化されると，提案アル} ゴリズムも高速化される．

2

_{については，これも}

Level-3 BLAS

_{を利用した高速な方式を開発} 中（早大・尾崎氏との共同研究）．

(25)

数値実験

逆

LU

_分解／逆

QR

_{分解のふるまいについて数値実験する．}

•

倍精度を基本精度とする（

u = 2

−53

≈ 1.1 · 10

−16_）

•

テスト行列：

Rump

_の行列（

INTLAB

_の

randmat(n,cnd)

）

• n = 100

，

cnd = 10

100（

A

∈ F

100×100，

κ(A)

≈ 1.75 · 10

107）

(26)

Table 1: Rump

_の行列（

n = 100 and κ(A)

≈ 1.75 · 10

107_）に対

する逆

LU

_{分解の結果}

k

κ(U

_k

)

κ(T

_k

)

κ(AX

_1:k

)

u

k

κ(A)

1

3.50 · 10

18

3.50 · 10

18

2.37 · 10

93

1.94 · 10

91

2

5.28 · 10

18

5.28 · 10

18

2.18 · 10

78

2.16 · 10

75

3

4.01 · 10

18

4.01 · 10

18

1.79 · 10

64

2.39 · 10

59

4

4.85 · 10

18

4.85 · 10

18

3.45 · 10

48

2.66 · 10

43

5

1.99 · 10

18

1.99 · 10

18

6.77 · 10

33

2.95 · 10

27

6

1.16 · 10

18

1.16 · 10

18

1.30 · 10

18

3.27 · 10

11

7

2.73 · 10

17

2.73 · 10

17

3.96 · 10

2

< 1

8

1.91 · 10

2

1.91 · 10

2

8.68 · 10

1

< 1

(27)

Table 2: Rump

_の行列（

n = 100 and κ(A)

≈ 1.75 · 10

107_）に対する逆

QR

_{分解の結果}

k

κ(R

_k

)

κ(T

_k

)

κ(AX

_1:k

)

u

k

κ(A)

1

3.27 · 10

19

3.27 · 10

19

2.00 · 10

93

1.94 · 10

91

2

1.86 · 10

19

1.86 · 10

19

1.06 · 10

77

2.16 · 10

75

3

7.97 · 10

17

7.97 · 10

17

1.61 · 10

62

2.39 · 10

59

4

2.20 · 10

17

2.20 · 10

17

4.23 · 10

46

2.66 · 10

43

5

2.31 · 10

17

2.31 · 10

17

2.00 · 10

32

2.95 · 10

27

6

4.04 · 10

17

4.04 · 10

17

6.69 · 10

16

3.27 · 10

11

7

2.18 · 10

18

2.18 · 10

18

1.39 · 10

3

< 1

8

1.39 · 10

3

1.39 · 10

3

1.00 · 10

0

< 1

(28)

応用

(1):

_{連立一次方程式}

連立一次方程式の近似解の計算に提案アルゴリズムを適用：

1. A

T の逆

LU

分解（

Crout

型なら

A

でよい）

P A

T

X

_U

≈ L

⇔

X

_UT

AP

≈ L

T

2. ey = [X

_UT

· b]

1_m_を計算

3.

_{三角方程式}

L

T

z =

ey

_{を解く（得られた近似解を}

ez

_とする）

4. ex = P ez

_を計算

(29)

数値実験

(1): (

_{スケール化}

)Hilbert

_行列

H

_n

H

_n_{は要素が整数の行列（}

n

≤ 21

_{のとき，倍精度で正確に表現可能）}

• n = 15 (κ(H

15

)

≈ 6.12 × 10

20

)

•

右辺：

b = H

₁₅

e

∈ F

15_，

e := (1, . . . , 1)

T

• H

15

x = b

の真の解：

H

₁₅−1

b = e = (1, . . . , 1)

T

(Matlab demo)

(30)

数値実験

(2): Rump

_の行列

様々な条件数に対し，提案アルゴリズムを用いて得られた連立一次方程式の近似解の相対精度（ノルム評価）をみる．

• Rump

の行列

randmat

（

n = 500

，条件数は

10

20から

10

100）

• b = (1, . . . , 1)

T

•

提案アルゴリズムの許容誤差：

ε

_tol

= 10

−6

•

多倍長精度との速度比較のため，

GMP

_{ベースのガウスの消去法} も実行（真の解は既知とするトライアンドエラー方式）

(31)

function [p,rel_err] = test_gmp_lin(A,b,xt,tol)

%

xt: given exact solution of Ax = b

% tol: tolerance for relative error

d = 53; norm_xt = norm(double(xt));

while 1

xv = gmp_lin(A,b,d);

% solve Ax=b using GMP

% normwise relative error

rel_err = norm(double(xt-xv))/norm_xt;

if rel_err < tol, break, end

d = 2*d;

% d = 53, 106, 212, ...

end

(32)

Table 3: Rump

_{の行列の結果（}

n = 500

_，

ε

_tol

= 10

−6_）提案アルゴリズム

GMP-based GEPP

κ(A)

ε

₁

t

₁

k

ε

₂

t

₂

d/53

1.31

20

2.84 · 10

−13

1.28

2

5.83 · 10

−14

23.35

2

1.73

33

2.44 · 10

−15

3.95

3

1.11 · 10

−16

65.76

4

2.70

43

4.17 · 10

−15

8.87

4

1.11 · 10

−16

69.23

4

1.13

49

1.14 · 10

−13

9.23

4

1.11 · 10

−16

71.47

4

4.97

59

3.05 · 10

−15

16.27

5

1.04 · 10

−7

70.98

4

7.05

73

5.01 · 10

−15

25.94

6

1.11 · 10

−16

154.70

8

3.85

81

7.54 · 10

−11

26.39

6

1.11 · 10

−16

157.12

8

2.74

92

5.38 · 10

−13

39.88

7

1.11 · 10

−16

158.12

8

−15 −16

(33)

応用例

(2):

_{行列式の計算}

良い方法（

Higham in ASNA

_{）として知られているのは，}

LU

_分解して

U

_ii_{の積を計算：}

det(A)

≈ det(P

T

LU ) = det(P ) det(U ).

（もちろん，これは近似計算なので，結果の保証がない）行列式は，その符号判定の情報が重要な例が多い．

(34)

行列式の精度保証

X

_L

, X

_U

: L, U

の近似逆行列．

B := X

_L

P AX

_U

⇒ det(B) = det(P ) det(A) det(X

_U

)

⇒

det(A) = det(P )

det(B)

/ det(X

_U

)

(6)

=

⇒ B:

_{単位行列（対角行列）に近い}

, det(B)

≈ 1

=

⇒

Gershgorin

_{の定理により，}

B

_{の固有値の範囲がわかる}

=

⇒ det(B)

_{の包含が可能}

(35)

しかしながら，

A

が悪条件な場合，

L

も

U

も良い

LU

分解因子とならない．

=

⇒

_{その場合，}

B = X

_L

P AX

_U _{は単位行列（対角行列）とかけ離} れてしまう．

=

⇒ Gershgorin

_{の定理が意味をなさない．} この問題も，提案アルゴリズムで解決できる！

(Matlab demo)

(36)

今後の課題

今回のコンセプトを用いると

•

高精度な逆

Cholesky

_分解

•

高精度な逆

LDL

T_分解

•

高精度な特異値分解も構成できる（多少の修正は必要．解析も別）．（固有値分解は少し難しく，現在，検討中）

(37)

2009 4

荻田 武史

概要

LU

QR

Cholesky

A: n

× n

A

: A

=

⇒

IEEE 754

10

=

⇒

:

条件数

A

A

κ(A) :=

∥A∥ · ∥A

∥.

=

⇒ κ(A)

=

⇒

Ax = b

x

:= A

b

Ay = b +

∆b

y

:= A

(b + ∆b)

∆x = y

− x

∥∆x∥

∥x

∥

≤

κ(A)

∥∆b∥

∥b∥

IEEE 754

=

⇒ 1

: u

≈ 1.11 × 10

=

⇒

10

κ(A) > u

1

=

⇒

=

⇒

“

”

=

⇒

“

”

Notation

A = (a

), C = (c

)

∈ R

• |A|

A

|A| = (|a

|)

• A ≤ C

(i, j)

a

≤ c

LU

分解（部分軸交換付）

荻田武史

_∥

_{分解（部分軸交換付）}