Cholesky 分解ノート

1 Cholesky .. 2008 6 9 , 2013 8 25 .. Cholesky .. Sylvester . L . U . U . 1 . , . Cholesky LU . ( . Hermite Hermite ) . LU . LU Cholesky . Cholesky . A . A = LDU (L , D , U ). LDU . A = AT = U T DT LT = U T DLT. LDU . L = UT , U = LT.. A = LDLT . 1. D D1/2 L . (LD1/2 L ). A = LLT.. A . P AP T = LDLT (P , L , D ). ( [2], [3] ) A 0 .. A = LDLT (L , D ). ( Cholesky ) A . A = LLT (L ). ( Cholesky ) . Trefethen and Bau III [9] [7], [5] . 2 Cholesky .. A = (aij ) M (N ; R) L = ( ij ) M (N ; R) .. A = LLT. A Cholesky ( Cholesky factorization ) 1 . Cholesky (1). Cholesky LU . 1. decompose . factorization .. 2.. (LU ) N N . (1) A GL(N ; R) P GL(N ; R), L GL(N ; R), U GL(N ; R) . (1) P A = LU.. (2) a 0 A GL(N ; R) . L GL(N ; R), U GL(N ; R) . (2) A = LU.. (3) (1), (2) LU . LU = L DU (L , U , D ).

2 ( LDU ) . a MATLAB det A(1 : r, 1 : r) A r .. A = (aij ) GL(N ; R) = 0 . LDU L = ( ij ) GL(N ; R), U = (uij ) GL(N ; R), D = diag (d1 , d2 , , dN ) GL(N ; R) .. A = LDU.. LDU = A = AT = (LDU )T = U T DT LT = U T DLT (U T , LT ). LDU . U = LT . A = LDLT Sylvester D di . (L 1 A(L 1 )T = D ) . ( ). D1/2 = diag d1 , d2 , , dN. D ( (D1/2 )2 = D) . b = LD1/2. L.. bL. L bT = (LD1/2 )(LD1/2 )T = LD1/2 (D1/2 )T LT = LD1/2 D1/2 LT = LDLT = A. 3. Cholesky . A = LLT L D . ( ) L .. Cholesky . Cholesky Cholesky .. A = LLT.. N. (3) aij = ik jk (i, j = 1, 2, , N ). k=1. L = ( ij ) . i<j = ij = 0. (3) .. min(i,j). aij = ik jk (i, j = 1, 2, , N ). k=1. A i j (3) .. min(i,j). aij = ik jk (1 j i N ). k=1. i . a11 = 211 , a21 = 21 11 , a22 = 221 + 222 , a31 = 31 11 , a32 = 31 21 + 32 22 , a33 = 231 + 232 + 233 , ai1 = i1 11 , ai2 = i1 21 + i2 22 , , aii = 2i1 + 2i2 + + 2ii.

3 AN 1 = N 1 11 , aN 2 = N 1 21 + N 2 22 , , aN N = 2N 1 + 2N 2 + + 2N N . 4.. 11 = a11 , . 21 = a21 / 11 , 22 = a22 221 , . 31 = a31 / 11 , 32 = (a32 31 21 ) / 22 , 33 = a33 231 232 , .. i1 = ai1 / 11 , i2 = (ai2 i1 21 ) /ai2 , ( ).. j 1. ij = aij ik jk / jj (j = 1, 2, , i 1), k=1.. ii = aii 2i1 2i2 2i,i 1 , .. N 1 = aN 1 / 11 , N 2 = (aN 2 N 1 21 ) /aN 2 , ( ).. j 1. N j = aN j N k jk / jj (j = 1, 2, , N 1), k=1.. N N = aN N 2N 1 2N 2 2N,N 1. ( ). A = (aij ) (i j ) L = ( ij ) . (i j ) ( 0 ) . Cholesky (ijk ) . for (i = 1; i <= N; i++) {. /* L[i][j] (i>j) */. for (j = 1; j < i; j++) {. s = a[i][j];. for (k = 1; k < j; k++). s -= L[i][k] * L[j][k];. L[i][j] = s / L[j][j];. }. /* L[i][i] */. s = a[i][i];. for (k = 1; k < i; k++). s -= sqr(L[i][k]);. L[i][i] = sqrt(s);. }.. 11 , 21 , 22 , 31 , , 33 , , N 1 , N 2 , , N N ( ).

4 Jik 11 , 21 , , N 1 , 22 , 23 , , N 2 , , N 1,N 1 , N,N 1 , N N ( ) ( ) . 5. Cholesky (jik ) . for (j = 1; j <= n; j++) {. /* L[j][j] */. s = a[j][j];. for (k = 1; k < j; k++). s -= sqr(L[j][k]);. if (s < 0) {. fprintf(stderr, "s < 0\n");. exit(0);. }. L[j][j] = sqrt(s);. /* L[i][j] (i>j) */. for (i = j + 1; i <= n; i++) {. s = a[i][j];. for (k = 1; k < j; k++). s -= L[i][k] * L[j][k];. L[i][j] = s / L[j][j];. }. }. }.. ( [5], Trefethen and Bau III [9], Higham [1] . ) A = LLT A = RT R . ( R = LT ) . i j . Cholesky A = RT R . RT R = A i j (i, j) . r1i r1j + r2i r2j + + rii rij = aij v u ( ). u . i 1 . i 1. t rii = aii 2. rkj , rij = aij rki rkj /rii (i > j). k=1 k=1. (A aij ).. L A = LLT A Cholesky . L . 6. 1 /*. 2 * 3 */. 4. 5 #include < >. 6 #include < >. 7 #include < >. 8 #include < > /* INT_MAX */.

5 9. 10 /* [0,1) */. 11 double drandom(). 12 {. 13 return random() / (double) INT_MAX;. 14 }. 15. 16 /* */. 17 void mul_mm(int n, matrix ab, matrix a, matrix b). 18 {. 19 int i, j, k;. 20 double s;. 21 for (i = 1; i <= n; i++). 22 for (j = 1; j <= n; j++) {. 23 s = 0;. 24 for (k = 1; k <= n; k++). 25 s += a[i][k] * b[k][j];. 26 ab[i][j] = s;. 27 }. 28 }. 29. 30 /* */. 31 void print_matrix(int n, matrix a). 32 {. 33 int i, j;. 34 for (i = 1; i <= n; i++) {. 35 for (j = 1; j <= n; j++). 36 printf("%f ", a[i][j]);. 37 printf("\n");. 38 }. 39 }. 40. 41 void clear_m(int n, matrix a). 42 {. 43 int i, j;. 44 for (i = 1; i <= n; i++). 45 for (j = 1; j <= n; j++). 46 a[i][j] = 0;. 47 }. 48. 49 /* */. 50 double sqr(double x) { return x * x; }. 51. 52 /* Cholesky */. 53 void Cholesky (int n, matrix L, matrix a).]

6 54 {. 7. 55 int i, j, k;. 56 double s;. 57 for (i = 1; i <= n; i++) {. 58 for (j = 1; j < i; j++) {. 59 s = a[i][j];. 60 for (k = 1; k < j; k++). 61 s -= L[i][k] * L[j][k];. 62 L[i][j] = s / L[j][j];. 63 }. 64 /* L[i][i] */. 65 s = a[i][i];. 66 for (k = 1; k < i; k++). 67 s -= sqr(L[i][k]);. 68 if (s < 0) {. 69 fprintf(stderr, "s < 0\n");. 70 exit(0);. 71 }. 72 L[i][i] = sqrt(s);. 73 }. 74 }. 75. 76 int main(). 77 {. 78 int i, j, N;. 79 matrix a, L, Lt;. 80 printf("N="); scanf("%d", . 81 a = new_matrix(N+1, N+1);. 82 L = new_matrix(N+1, N+1);. 83 Lt = new_matrix(N+1, N+1);. 84 for (i = 1; i <= N; i++) {. 85 for (j = 1; j <= i; j++). 86 Lt[j][i] = L[i][j] = drandom();. 87 for (j = i + 1; j <= N; j++). 88 Lt[j][i] = L[i][j] = ;. 89 }. 90 mul_mm(N, a, L, Lt);. 91 printf("L=\n").)}

7 92 print_matrix(N, L);. 93 printf("a=\n");. 94 print_matrix(N, a);. 95. 96 clear_m(N, L);. 97 Cholesky (N, L, a);. 98 printf("L=\n");. 99 print_matrix(N, L);. 100. 101 return 0;. 102 }. 8. cholesky1 .. oyabun% ./cholesky1. N=4. L=. a=. L=. oyabun%.. jik Cholesky .. 1 /* Cholesky */. 2 void Cholesky (int n, matrix L, matrix a). 3 {. 4 int i, j, k;. 5 double s;. 6 for (j = 1; j <= n; j++) {. 7 /* L[j][j] */. 8 s = a[j][j];. 9 for (k = 1; k < j; k++). 10 s -= sqr(L[j][k]);. 11 if (s < 0) {. 12 fprintf(stderr, "s < 0\n");. 13 exit(0);. 14 }. 15 L[j][j] = sqrt(s);. 16 for (i = j + 1; i <= n; i++) {. 17 s = a[i][j];. 18 for (k = 1; k < j; k++). 19 s -= L[i][k] * L[j][k];. 20 L[i][j] = s / L[j][j];. 21 }. 22 }. 23 }. 9. jik .. oyabun% ./cholesky2. N=4. L=. a=. L=. oyabun%.. Cholesky .

8 U = LT L L[i][j] u[j][i] . A . ( ) . 1 /* Cholesky A=U^T U (A ) */. 2 void Cholesky (int n, matrix U, matrix a). 3 {. 4 int j, i, k;. 5 double s;. 6 for (i = 1; i <= n; i++) {. 7 /* U[i][i] */. 8 s = a[i][i];. 9 for (k = 1; k < i; k++). 10 s -= sqr(U[k][i]);. 11 if (s < 0) {. 12 fprintf(stderr, "s < 0\n");. 13 exit(0);. 14 }. 15 U[i][i] = sqrt(s);. 16 for (j = i + 1; j <= n; j++) {. 17 s = a[i][j];. 18 for (k = 1; k < i; k++). 19 s -= U[k][j] * U[k][i];. 20 U[i][j] = s / U[i][i];. 21 }. 22 }. 23 }. 10. Cholesky (2). Trefethen and Bau III [9] .. ( ) A (1, 1) 1 Gauss 1 . ( ) ( )( ). 1 wT 1 0 1 wT. A= =. w K w I 0 K wwT.. ( ) ( )( ). 1 wT 1 0 1 wT. =. 0 K wwT 0 K wwT 0 I.. ( )( )( ). 1 0 1 0 1 wT. A= . w I 0 K wwT 0 I.. ( ). a11 wT. A=. w B.. = a11 . ( )( )( ). 0 1 0T 0. A= = R1T A1 R1.

9 W/ I 0 K wwT /a11 w/ I.. A1 = R2T A2 R2 , A2 = R3T A3 R3 , , AN 1 = RN. T. AN RN. AN = I (N ) . A = R1T R2T RN. T. RN R2 R1 . R = RN R2 R1 . A = RT R.. Cholesky . R=A. for k = 1 to m for j=k+1 to m Rj,j:m = Rj,j:m Rk,j:m Rkj /Rkk . Rk,k:m = Rk,k:m / Rkk . 11.. N 3 /6 + O(N 2 ), N 3 /3 + O(N 2 ) Gauss .. N O(N ) .. LU L U . 1 Cholesky 1 . N (O(N 2 ) N . ) Cholesky .. Trefethen and Bau III [9] ( . ) Cholesky . Cholesky Gauss Cholesky ( ) . R 2 . R 2 = RT 2 = A 1/2. ( ).. p (1 p ) N . 1 . A 1/2 R = RT N A 1/2 . N. (Cf. LU partial pivoting ( L = 1 ) U =. 2N 1 A ).. ( Cholesky ) A N . (A), (B) Cholesky . M R e . eT R. e = A + A, A . (4) A M (N ; R) R = O( M ) ( M 0). A .. O( ) Higham [1] Theorem . eT | |R|. | A| N +1 |R e n . nu n = (u the unit round o ). 1 nu IEEE 754 u = M . 12.

10 E R . R.. e R . R. = O( (A) M ). R .. Ax = b . 1 LU .. A Ax = b A A = LLT Cholesky . LLT x = b y . Ly = b . LT x = y x O(N 2 ) . ( ) .. ( Cholesky 1 ) . 1 Ax = b (A), (B) . Cholesky . x e . A . A M (N ; R) (A + A)e x = b, = O( M ). A .. Higham [1] , Theorem . (O( ) ) . eT | |R|. | A| 3N +1 |R e . A 2 3N +1 N (1 N N +1 ) 1 A 2.. ( 1 . ). 13. 3 Cholesky . ( . ).. ( Cholesky ) N A = (aij ) k (k = 1, 2, , N ) 0 . A = LDLT (L , D ).. 1 . d1.. 21 1. 0 . d2. 0 .. L= . 31 32 1 , D= .. N 1 N 2 N,N 1 1. 0 dN. LDLT (i, j) .. d1 j1.. d2 j2 .. min{i,j}. d . j 1 j,j 1 . ( i,1 i,2 i,i 1 1 0 0) = dk ik jk ( ii = 1). dj . k=1. 0 .. 0. A i j ii = 1.. i 1.. dk 2ik i di + (j = i). k=1. aij = dk ik jk =.. i 1. k=1 . dk ik jk + di ji (j = i + 1, i + 2, , N ).. k=1. i = 1, 2, , N .. i 1. di = aii dk 2ik , k=1.