D D 1/2 L (LD 1/2 L ) A = LL T A P AP T = LDL T (P, L, D ) ( [2], [3] ) A 0 A = LDL T (L, D ) ( Cholesky ) A A = LL T (L ) ( Cholesky ) Trefethen and

Transcription

1 Cholesky , Cholesky Sylvester L U U 1, Cholesky LU ( Hermite Hermite ) LU LU Cholesky Cholesky A A = LDU (L, D, U ) LDU LDU A = A T = U T D T L T = U T DL T L = U T, U = L T A = LDL T. 1

2 D D 1/2 L (LD 1/2 L ) A = LL T A P AP T = LDL T (P, L, D ) ( [2], [3] ) A 0 A = LDL T (L, D ) ( Cholesky ) A A = LL T (L ) ( Cholesky ) Trefethen and Bau III [9] [7], [5] 2 Cholesky 2.1 A = (a ij ) M(N; R) L = (l ij ) M(N; R) A = LL T A Cholesky (Cholesky factorization) Cholesky (1) Cholesky LU 1 decompose factorization 2

3 2.1 (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 ) ( LDU ) a MATLAB det A(1 : r, 1 : r) A r A = (a ij ) GL(N; R) 0 LDU L = (l ij ) GL(N; R), U = (u ij ) GL(N; R), D = diag (d 1, d 2,, d N ) GL(N; R) A = LDU LDU = A = A T = (LDU) T = U T D T L T = U T DL T (U T, L T ) LDU U = L T. A = LDL T Sylvester D d i (L 1 A(L 1 ) T = D ) ( D 1/2 = diag d1, d 2,, ) d N D ( (D 1/2 ) 2 = D) L = LD 1/2 L L T = (LD 1/2 )(LD 1/2 ) T = LD 1/2 (D 1/2 ) T L T = LD 1/2 D 1/2 L T = LDL T = A. 3

4 2.3 Cholesky A = LL T L D ( ) L 2.4 Cholesky 2.2 Cholesky Cholesky A = LL T (3) a ij = N l ik l jk (i, j = 1, 2,, N) k=1 L = (l ij ) i < j = l ij = 0 (3) a ij = min(i,j) k=1 l ik l jk (i, j = 1, 2,, N) A i j (3) a ij = i a 11 = l 2 11, min(i,j) k=1 a 21 = l 21 l 11, a 22 = l l 2 22, l ik l jk (1 j i N). a 31 = l 31 l 11, a 32 = l 31 l 21 + l 32 l 22, a 33 = l l l 2 33, a i1 = l i1 l 11, a i2 = l i1 l 21 + l i2 l 22,, a ii = l 2 i1 + l 2 i2 + + l 2 ii,. a N1 = l N1 l 11, a N2 = l N1 l 21 + l N2 l 22,, a NN = l 2 N1 + l 2 N2 + + l 2 NN. 4

5 l 11 = ± a 11, l 21 = a 21 /l 11, l 22 = ± a 22 l 2 21, l 31 = a 31 /l 11, l 32 = (a 32 l 31 l 21 ) /l 22, l 33 = ± a 33 l 2 31 l 2 32,. l i1 = a i1 /l 11, l i2 = (a i2 l i1 l 21 ) /a i2, ( ) j 1 l ij = a ij l ik l jk /l jj (j = 1, 2,, i 1),. k=1 l ii = ± a ii l 2 i1 l2 i2 l2 i,i 1, l N1 = a N1 /l 11, l N2 = (a N2 l N1 l 21 ) /a N2, ( ) j 1 l Nj = a Nj l Nk l jk /l jj (j = 1, 2,, N 1), k=1 l NN = ± a NN l 2 N1 l2 N2 l2 N,N 1 ( ). A = (a ij ) (i j ) L = (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); } l 11, l 21, l 22, l 31,, l 33,, l N1, l N2,, l NN ( ) jik l 11, l 21,, l N1, l 22, l 23,, l N2,, l N 1,N 1, l N,N 1, l NN ( ) ( ) 5

6 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 = LL T A = R T R ( R = L T ) i j Cholesky A = R T R R T R = A i j (i, j) r 1i r 1j + r 2i r 2j + + r ii r ij = a ij ( ) i 1 r ii = aii rkj 2, r i 1 ij = a ij r ki r kj /r ii (i > j). k=1 (A a ij ) k=1 2.5 L A = LL T A Cholesky L 6

7 2.5.1 cholesky1.c 1 /* 2 * cholesky1.c 3 */ 4 5 #include <stdio.h> 6 #include <math.h> 7 #include <matrix.h> 8 #include <limits.h> /* INT_MAX */ 9 10 /* [0,1) */ 11 double drandom() 12 { 13 return random() / (double) INT_MAX; 14 } /* */ 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 } /* */ 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 } 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 } /* */ 50 double sqr(double x) { return x * x; } /* Cholesky */ 53 void cholesky(int n, matrix L, matrix a) 54 { 7

8 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 } int main() 77 { 78 int i, j, N; 79 matrix a, L, Lt; 80 printf("n="); scanf("%d", &N); 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] = 0.0; 89 } 90 mul_mm(n, a, L, Lt); 91 printf("l=\n"); 92 print_matrix(n, L); 93 printf("a=\n"); 94 print_matrix(n, a); clear_m(n, L); 97 cholesky(n, L, a); 98 printf("l=\n"); 99 print_matrix(n, L); return 0; 102 } 8

9 2.5.2 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

10 2.5.4 jik oyabun%./cholesky2 N=4 L= a= L= oyabun% 2.6 Cholesky U = L T 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

11 2.7 Cholesky (2) Trefethen and Bau III [9] ( ) A (1, 1) 1 Gauss 1 ( ) ( ) ( ) 1 w T w T A = = w K w I 0 K ww T ( 1 w T 0 K ww T ) = ( K ww T ) ( 1 w T 0 I ) A = ( 1 0 w I ) ( A = K ww T ( a 11 α = a 11 ( ) ( α T A = w/α I 0 K ww T /a 11 w w T B ) ) ( ) ( 1 w T 0 I α 0 w/α I ) ). = R T 1 A 1 R 1 A 1 = R T 2 A 2 R 2, A 2 = R T 3 A 3 R 3,, A N 1 = R T NA N R N A N = I (N ) A = R T 1 R T 2 R T NR N R 2 R 1. R = R N R 2 R 1 Cholesky R = A for k = 1 to m for j=k+1 to m R j,j:m = R j,j:m R k,j:m R kj /R kk R k,k:m = R k,k:m / R kk A = R T R 11

12 2.8 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 2.9 Trefethen and Bau III [9] ( ) Cholesky 2.7 Cholesky Gauss Cholesky ( 2.2 ) R 2 R 2 = R T 2 = A 1/2 ( ) p (1 p ) N 1 N A 1/2 R = R T N A 1/2. (Cf. LU partial pivoting ( L = 1 ) U = 2 N 1 A ) 2.2 (Cholesky ) A N (A), (B) Cholesky ε M R (4) δa M(N; R) s.t. RT R = A + δa, δa A = O(ε M) (ε M 0). 2.3 O( ) Higham [1] p.197 Theorem 10.3 δa γ N+1 R T R γ n γ n = nu (u the unit round off) 1 nu IEEE 754 u = ε M 12

13 R R R R = O(κ(A)ε M ) R 2.10 Ax = b 1 LU A Ax = b A A = LL T Cholesky LL T x = b y Ly = b L T x = y x O(N 2 ) ( ) 2.4 (Cholesky 1 ) 1 Ax = b (A), (B) Cholesky x A M(N; R) s.t. (A + A) x = b, A A = O(ε M) 2.5 Higham [1] p.198, Theorem 10.4 (O( ) ) A γ 3N+1 R T R A 2 γ 3N+1 N(1 Nγ N+1 ) 1 A 2 ( 1 ) 13

14 3 Cholesky ( ) 3.1 ( Cholesky ) N A = (a ij ) δ k (k = 1, 2,, N) 0 A = LDL T (L, D ) 1 d 1 l L = l 31 l 32 1 d 2, D = l N1 l N2 l N,N 1 0 d N 1 LDL T (i, j) (l i,1 l i,2 l i,i ) d 1 l j1 d 2 l j2. d j 1 l j,j 1 d j 0. 0 = min{i,j} k=1 d k l ik l jk ( l ii = 1). A i j l ii = 1 i 1 i d i + d k l 2 ik (j = i) k=1 a ij = d k l ik l jk = i 1 k=1 d k l ik l jk + d i l ji (j = i + 1, i + 2,, N) k=1 i = 1, 2,, N i 1 d i = a ii d k l 2 ik, k=1 ) l ji = 1 i 1 (a ij d k l ik l jk d i k=1 14 (j = i + 1,, N)

15 {d i }, {l ji } i j d i Gauss i d (i) ii 0 d 1 = a 11 = δ 1 0, l j1 = a j1 /d 1 (j = 2,, N), d 2 = a 22 d 1 l 2 21 = δ 2 /δ 1 0, l j2 = (a 2j d 1 l 21 l j1 )/d 2 (j = 3,, N), d 3 = a 33 d 1 l 2 31 d 2 l 2 32 = δ 3 /δ 2 0, l j3 = (a 3j d 1 l 31 l j1 d 2 l 32 l j2 )/d 3 (j = 4,, N), A Cholesky LU ( [6] ) 4 ( ) Cholesky LU A Cholesky A = L T L 1 Ax = b Ly = b, L T x = y 15

16 1 /* */ 2 /* Ux=b */ 3 void solve_uxb(int n, matrix u, vector b) 4 { 5 int i, j; 6 for (i = n; i >= 1; i--) { 7 for (j = i + 1; j <= n; j++) 8 b[i] -= u[i][j] * b[j]; 9 b[i] /= u[i][i]; 10 } 11 } /* L^T x=b solve_uxb() u[i][j] L[j][i] */ 14 void solve_ltxb(int n, matrix L, vector b) 15 { 16 int i, j; 17 for (i = n; i >= 1; i--) { 18 for (j = i + 1; j <= n; j++) 19 b[i] -= L[j][i] * b[j]; 20 b[i] /= L[i][i]; 21 } 22 } /* Lx=b */ 25 void solve_lxb(int n, matrix L, vector b) 26 { 27 int i, j; 28 for (i = 1; i <= n; i++) { 29 for (j = 1; j < i; j++) 30 b[i] -= L[i][j] * b[j]; 31 b[i] /= L[i][i]; 32 } 33 } /* U^T x=b */ 36 void solve_utxb(int n, matrix U, vector b) 37 { 38 int i, j; 39 for (i = 1; i <= n; i++) { 40 for (j = 1; j < i; j++) 41 b[i] -= U[j][i] * b[j]; 42 b[i] /= U[i][i]; 43 } 44 } 16

17 oyabun%./cholesky3 N=5 L= a= L= x= oyabun% A Cholesky A.1 C 1 /* 2 * band_cholesky.c --- A Cholesky 3 */ 4 5 #include <stdio.h> 6 #include <limits.h> 7 #include <math.h> 8 #include <matrix.h> 9 10 double drandom() { return random() / (double) INT_MAX; } 11 double sqr(double x) { return x * x; } 12 double max(double a, double b) { return (a > b)? a : b; } 13 double min(double a, double b) { return (a < b)? a : b; } /* Cholesky A=U^T U (A, U ) 16 * 17 * m n A=(a_{ij}) (1 i n, 1 j n) 18 * ASB=(asb_{ij) (1 i n, 1 j m+1) 19 * 20 * 21 * asb[i][j-i+1] = a[i][j] (1 i n, i j min(i+m,n)) 22 * 23 * U ASB 24 * 17

18 25 */ 26 void band_cholesky(int n, int m, matrix ASB) 27 { 28 int j, i, k; 29 double s; 30 for (i = 1; i <= n; i++) { 31 /* U[i][i] */ 32 s = ASB[i][1]; 33 for (k = max(1,i-m); k < i; k++) 34 s -= sqr(asb[k][i-k+1]); 35 if (s < 0) { 36 fprintf(stderr, "s < 0\n"); 37 exit(0); 38 } 39 ASB[i][1] = sqrt(s); 40 /* U[i][j] (i < j) */ 41 for (j = i + 1; j <= min(i + m, n); j++) { 42 s = ASB[i][j-i+1]; 43 for (k = max(1,j-m); k < i; k++) 44 s -= ASB[k][j-k+1] * ASB[k][i-k+1]; 45 ASB[i][j-i+1] = s / ASB[i][1]; 46 } 47 } 48 } /* */ 51 /* U^T x=b */ 52 void solve_band_utxb(int n, int m, matrix UB, vector b) 53 { 54 int i, j; 55 for (i = 1; i <= n; i++) { 56 for (j = max(i-m,1); j < i; j++) 57 b[i] -= UB[j][i-j+1] * b[j]; 58 b[i] /= UB[i][1]; 59 } 60 } /* Ux=b */ 63 void solve_band_uxb(int n, int m, matrix UB, vector b) 64 { 65 int i, j; 66 for (i = n; i >= 1; i--) { 67 for (j = i + 1; j <= min(i+m,n); j++) 68 b[i] -= UB[i][j-i+1] * b[j]; 69 b[i] /= UB[i][1]; 70 } 71 } /* */ 74 void print_symmetric_band_matrix1(int n, int m, matrix asb) 75 { 76 int i, j; 77 for (i = 1; i <= n; i++) { 78 for (j = 1; j <= n; j++) 79 if (abs(j - i) > m) 80 printf("%f ", 0.0); 18

19 81 else { 82 if (j >= i) 83 printf("%f ", asb[i][j-i+1]); 84 else /* i j */ 85 printf("%f ", asb[j][i-j+1]); 86 } 87 printf("\n"); 88 } 89 } /* */ 92 void print_upper_triangular_band_matrix1(int n, int m, matrix ub) 93 { 94 int i, j; 95 for (i = 1; i <= n; i++) { 96 for (j = 1; j <= n; j++) 97 if (abs(j - i) > m) 98 printf("%f ", 0.0); 99 else { 100 if (j >= i) 101 printf("%f ", ub[i][j-i+1]); 102 else 103 printf("%f ", 0.0); 104 } 105 printf("\n"); 106 } 107 } void print_vector1(int n, vector x) 110 { 111 int i; 112 for (i = 1; i <= n; i++) 113 printf("%f ", x[i]); 114 printf("\n"); 115 } int main() 118 { 119 int n, m, i, j, k; 120 double s; 121 matrix ASB, UB; 122 vector x, b; /*, */ 125 n = 8; m = 3; /* ASB, UB n (m+1) */ 128 ASB = new_matrix(n + 1, (m+1) + 1); 129 UB = new_matrix(n + 1, (m+1) + 1); 130 x = new_vector(n + 1); 131 b = new_vector(n + 1); /* U */ 134 for (i = 1; i <= n; i++) 135 for (j = i; j <= min(i+m, n); j++) 136 UB[i][j-i+1] = drandom(); 19

20 137 printf("u=\n"); 138 print_upper_triangular_band_matrix1(n, m, UB); /* A=U^T U */ 141 for (i = 1; i <= n; i++) 142 for (j = i; j <= min(i+m,n); j++) { /* */ 143 s = 0; 144 /* U[k][i]*U[k][j] 145 max(i-m,1) k i max(j-m,1) k j i j max(j-m,1) k i 147 */ 148 for (k = max(j-m,1); k <= i; k++) 149 s += UB[k][i-k+1] * UB[k][j-k+1]; 150 ASB[i][j-i+1] = s; 151 } 152 printf("a:=u^t U\n"); 153 print_symmetric_band_matrix1(n, m, ASB); /* x */ 156 for (i = 1; i <= n; i++) 157 x[i] = i; /* b=a x */ 160 for (i = 1; i <= n; i++) { 161 s = ASB[i][1] * x[i]; 162 for (j = i + 1; j <= min(i+m,n); j++) 163 s += ASB[i][j-i+1] * x[j]; 164 for (j = max(i-m,1); j < i; j++) 165 s += ASB[j][i-j+1] * x[j]; 166 b[i] = s; 167 } /* Cholesky A=U^T U */ 170 band_cholesky(n, m, ASB); 171 printf("cholesky U=\n"); 172 print_upper_triangular_band_matrix1(n, m, ASB); /* U^T U x = b */ 175 solve_band_utxb(n, m, UB, b); 176 solve_band_uxb(n, m, UB, b); 177 printf("x=\n"); 178 print_vector1(n, b); return 0; 181 } 20

21 oyabun%./band_cholesky U= A:=U^T U Cholesky U= x= oyabun% A.2 C++ & Profil Profil 1 /* 2 * band_cholesky_c++.c --- A Cholesky 3 */ 4 5 #include <iostream.h> 6 #include <iomanip.h> 7 #include <limits.h> 8 #include <math.h> 9 #include <Matrix.h> 10 #include <Vector.h> double drandom(void) { return random() / (double) INT_MAX; } 13 double sqr(double x) { return x * x; } 14 // double max(double a, double b) { return (a > b)? a : b; } 15 // double min(double a, double b) { return (a < b)? a : b; } /* Cholesky A=U^T U (A, U ) 21

22 18 * 19 * m n A=(a_{ij}) (1 i n, 1 j n) 20 * ASB=(asb_{ij) (1 i n, 1 j m+1) 21 * 22 * 23 * asb[i][j-i+1] = a[i][j] (1 i n, i j min(i+m,n)) 24 * 25 * U ASB 26 * 27 */ 28 void band_cholesky(int n, int m, MATRIX &ASB) 29 { 30 double s; 31 for (int i = 1; i <= n; i++) { 32 /* U[i][i] */ 33 s = ASB(i,1); 34 for (int k = max(1,i-m); k < i; k++) 35 s -= sqr(asb(k,i-k+1)); 36 if (s < 0) { 37 cerr << "s < 0" << endl; 38 exit(0); 39 } 40 ASB(i,1) = sqrt(s); 41 /* U[i][j] (i < j) */ 42 for (int j = i + 1; j <= min(i + m, n); j++) { 43 s = ASB(i,j-i+1); 44 for (int k = max(1,j-m); k < i; k++) 45 s -= ASB(k,j-k+1) * ASB(k,i-k+1); 46 ASB(i,j-i+1) = s / ASB(i,1); 47 } 48 } 49 } /* */ 52 /* U^T x=b */ 53 void solve_band_utxb(int n, int m, const MATRIX &UB, VECTOR &b) 54 { 55 for (int i = 1; i <= n; i++) { 56 for (int j = max(i-m,1); j < i; j++) 57 b(i) -= UB(j,i-j+1) * b(j); 58 b(i) /= UB(i,1); 59 } 60 } /* Ux=b */ 63 void solve_band_uxb(int n, int m, const MATRIX &UB, VECTOR &b) 64 { 65 for (int i = n; i >= 1; i--) { 66 for (int j = i + 1; j <= min(i+m,n); j++) 67 b(i) -= UB(i,j-i+1) * b(j); 68 b(i) /= UB(i,1); 69 } 70 } /* */ 73 void print_symmetric_band_matrix1(int n, int m, const MATRIX &asb) 22

23 74 { 75 for (int i = 1; i <= n; i++) { 76 for (int j = 1; j <= n; j++) 77 if (abs(j - i) > m) 78 cout << 0.0 << " "; 79 else { 80 if (j >= i) 81 cout << asb(i,j-i+1) << " "; 82 else /* i j */ 83 cout << asb(j,i-j+1) << " "; 84 } 85 cout << endl; 86 } 87 } /* */ 90 void print_upper_triangular_band_matrix1(int n, int m, const MATRIX &ub) 91 { 92 for (int i = 1; i <= n; i++) { 93 for (int j = 1; j <= n; j++) 94 if (abs(j - i) > m) 95 cout << 0.0 << " "; 96 else { 97 if (j >= i) 98 cout << ub(i,j-i+1) << " "; 99 else 100 cout << 0.0 << " "; 101 } 102 cout << endl; 103 } 104 } void print_vector1(int n, const VECTOR &x) 107 { 108 for (int i = 1; i <= n; i++) 109 cout << x(i) << " "; 110 cout << endl; 111 } int main() 114 { 115 int n, m, i, j, k; 116 double s; cout << setprecision(4); 119 cout << setiosflags(ios::fixed); 120 /*, */ 121 n = 8; m = 3; /* ASB, UB n (m+1) */ 124 MATRIX ASB(n,m+1), UB(n,m+1); 125 VECTOR x(n), b(n); /* U */ 128 for (i = 1; i <= n; i++) 129 for (j = i; j <= min(i+m, n); j++) 23

24 130 UB(i,j-i+1) = drandom(); 131 cout << "U=" << endl; 132 print_upper_triangular_band_matrix1(n, m, UB); /* A=U^T U */ 135 for (i = 1; i <= n; i++) 136 for (j = i; j <= min(i+m,n); j++) { /* */ 137 s = 0; 138 /* U[k](i)*U[k](j) 139 max(i-m,1) k i max(j-m,1) k j i j max(j-m,1) k i 141 */ 142 for (k = max(j-m,1); k <= i; k++) 143 s += UB(k,i-k+1) * UB(k,j-k+1); 144 ASB(i,j-i+1) = s; 145 } 146 cout << "A:=U^T U" << endl; 147 print_symmetric_band_matrix1(n, m, ASB); /* x */ 150 for (i = 1; i <= n; i++) 151 x(i) = i; /* b=a x */ 154 for (i = 1; i <= n; i++) { 155 s = ASB(i,1) * x(i); 156 for (j = i + 1; j <= min(i+m,n); j++) 157 s += ASB(i,j-i+1) * x(j); 158 for (j = max(i-m,1); j < i; j++) 159 s += ASB(j,i-j+1) * x(j); 160 b(i) = s; 161 } /* Cholesky A=U^T U */ 164 band_cholesky(n, m, ASB); 165 cout << "Cholesky U=" << endl; 166 print_upper_triangular_band_matrix1(n, m, ASB); /* U^T U x = b */ 169 solve_band_utxb(n, m, UB, b); 170 solve_band_uxb(n, m, UB, b); 171 cout << "x=" << endl; 172 print_vector1(n, b); return 0; 175 } 24

25 oyabun% g++ -W -I/usr/local/include -o band_cholesky_c++ band_cholesky_c++.c -lprofil -lmatrix oyabun%./band_cholesky_c++ U= A:=U^T U Cholesky U= x= oyabun% B Trefethen and Bau III [9] ( ) F (A) F fl: R F (B) x R ε R s.t. ε ε M, fl(x) = x(1 + ε). F,,, {+,,, }, x, y F ε R s.t. ε ε M, x y = (x y)(1 + ε). 25

26 C 2 C.1 R N R N (, ) R N (, ) x = (x i ) R N, y = (y i ) R N (x, y) R N = (x, y) = N x i y i. i=1 y y T x : (x, y) R N = y T x. C.2 R N a: R N R N R 1. x R N, y R N, z R N a(x + y, z) = a(x, z) + a(y, z), a(x, y + z) = a(x, z) + a(y, z). 2. x R N, y R N, λ R a(λx, y) = λa(x, y), a(x, λy) = λa(x, y). N A = (a ij ) M(N; R) a(x, y) def. = (Ax, y) R N R N R N a a ij = a(e j, e i ) ( N ) N a x j e j, y i e i = j=1 i=1 A = (a ij ) N N x j y i a(e j, e i ) = j=1 i=1 ( N N ) a ij x j y i = (Ax, y) R N. i=1 j=1 C.3 R N x R N y R N a(x, y) = a(y, x) a R N a(x, y) = (Ax, y) R N a(x, y) = (Ax, y) R N = y T Ax = (A T y) T x = (x, A T y) R N = (A T y, x) R N, 26

27 a(y, x) = (Ay, x) R N a x R N y R N (A T y, x) R N = (Ay, x) R N A = A T. R N A (Ax, y) C.4 R N R N a(, ) x R N a(x, x) 0 ( x R N a(x, x) = 0 x = 0 ). a(, ) R N a A a(x, y) = (Ax, y) R N C.5 2 a: R N R R N 2 2 A = (a ij ) N N a(x) = a ij x i x j, i=1 j=1 a(x) = (Ax, x) R N (Ax, x) R N = (x, A T x) R N = (A T x, x) R N A A T (Ax, x) R N = 1 ( ) 1 2 ((Ax, x) R N + (AT x, x) R N ) = 2 (A + AT )x, x A (A + A T )/2 ( 2 x (Ax, x) A ) (Ax, x) = 1 2 (A + AT )x = Ax R N 27

28 R N 2 R N (x, y) (x, x) a(x, y) = 1 2 (a(x + y, x + y) a(x, x) a(y, y)) = 1 (a(x + y) a(x) a(y)) 2 2 C.6 A = (a ij ) M(N; R) 2 f(x) = (Ax, x) = x T Ax P GL(N; R), x = P y f(x) = f(p y) = (P y) T A(P y) = y T P T AP y = (P T AP y, y) = (A y, y). A := P T AP 2 (A ) A diag (λ 1,, λ N ) f(x) = N λ i yi 2. i=1 λ i A C.1 (Sylvester (Sylvester s law of inertia)) A n P D := P T AP n D A D Sylvester ( 2013/8/26 10:52:43) D.1 Sylvester 2 ( ) 28

29 2 (2 ) A P P T AP A P T AP A P 1 AP ( ) (P T = P 1 ) ( ) 1,2 2 Taylor 2 2 ( [4]) f n a f (a) = 0 f(a + h) = f(a) + h 2 + O( h 2 ) (h 0). 2 Morse Lemma 2 Sylvester D.2 3 ( ) A A 0 π(a), ν(a), ζ(a) A n π(a), ν(a), ζ(a) 0, π(a) + ν(a) + ζ(a) = n ( ) ( Sylvester ) A n U (U 1 = U T ) λ 1,, λ n U T AU = diag (λ 1,, λ n ). x R n y := U 1 x x = Uy (Ax, x) = (AUy, Uy) = (U T AUy, y) = n λ j yj 2. 2 (Ax, x) ( ) j=1

30 ( ) (i) U T AU diag (λ 1,, λ n ) (ii) U (Ax, x) U ( U P ) U λ j A λ j (j = 1,, n) A Sylvester D.1 (Sylvester, (Sylvester s law of inertia)) A P P T AP = diag (λ 1,, λ n ) n p := {j {1,, n}; λ j > 0}, n n := {j {1,, n}; λ j < 0}, n z := {j {1,, n}; λ j = 0} ( ) P A ( 0 A ) D.2 ( Sylvester ) A P P T AP = diag (λ 1,, λ n ) n p := {j {1,, n}; λ j > 0}, n n := {j {1,, n}; λ j < 0}, n z := {j {1,, n}; λ j = 0} n p = π(a), n n = ν(a), n z = ζ(a). 30

31 D.1 P U T AU = diag[λ 1,, λ n ] U P T AP = U T AU = U 1 AU A λ 1,, λ n A n p, n n, n z π(a), ν(a), ζ(a) [8] Sylvester D.3 ( Sylvester ) A P B := P T AP π(a) = π(b), ν(a) = ν(b), ζ(a) = ζ(b) D.2 B ( ) D.2 D.3 ( ) D.3 n A, P λ 1,, λ n s.t. P T AP = diag[λ 1,, λ n ] B := P T AP ( ) λ 1,, λ n n p = π(b) = π(a), n n = ν(b) = ν(a), n z = ζ(b) = ζ(a). D.2 A, P B := P T AP B U : λ 1,, λ n R s.t. U T BU = diag[λ 1,, λ n ]. diag (λ 1,, λ n ) = U T BU = U T P T AP U = (P U) T A(P U) = P T AP, P := P U. P U P D.2 (B ) λ 1,, λ n 0 π(a), ν(a), ζ(a) π(b) = π(a), ν(b) = ν(a), ζ(b) = ζ(a). D.3 Sylvester D.1 D.1 ( ) x = Sy, x = Rz (Ax, x) = α 1 y α r y 2 r α r+1 y 2 r+1 α p y 2 p (α j > 0) = β 1 z β s z 2 s β s+1 z 2 s+1 β q z 2 q (β j > 0) 31

32 p q p = q, r = s y j, z k x y j z k 1 r < s y 1,, y p y 1,, y r, z s+1,, z q 1 ( p 1 y 1,, y p r + (q s) 1 r + (q s) = (r s) + q < q p ) y k, r + 1 k p y 1,, y r, z s+1,, z q y 1 = = y r = z s+1 = = z q = 0, y k = 1 x 1,, x n 1 (Ax, x) = α α r 0 2 α r+1 y 2 r+1 α k 1 2 α p y 2 p α k < 0, (Ax, x) = β 1 z β s z 2 s β s β q r > s y 1,, y r y r+1,, y p, z 1,, z s y k (1 k r) y r+1,, y p, z 1,, z s y r+1 = = y p = z 1 = = z s = 0, y k = 1 x 1,, x n 1 (Ax, x) = α 1 y α k α r y 2 r α r α p 0 2 α k > 0, (Ax, x) = β β s 0 2 β s+1 y 2 s+1 β q β 2 q 0 r = s. ( Ax, x) ( ) D.2 ( ) D.2 A n P := {V ; V R n, x V \ {0} (Ax, x) > 0}, N := {V ; V R n, x V \ {0} (Ax, x) < 0}, Z := {V ; V R n, x V \ {0} (Ax, x) = 0} {0} P, N, Z P, N, Z =. N p := max{dim V ; V P}, N n := max{dim V ; V N }, N z := max{dim V ; V Z} ( : dim{0} = 0 ) : N p = π(a), N n = ν(a), N z = ζ(a). 32

33 A v 1,, v n Av j = λ j v j V p := Span{v j ; λ j > 0}, V n := Span{v j ; λ j < 0}, V z := Span{v j ; λ j = 0} V p P, V n N, V z Z. dim V p = π(a), dim V n = ν(a), dim V z = ζ(a) N p, N n, N z ( ) N p π(a), N n ν(a), N z ζ(a). W p P, W n N, W z Z W p W n = W N W z = W z W p = {0} W p, W n, W z ( ) dim W p = N p, dim W n = N n, dim W z = N z W p W n W z R n n = π(a) + ν(a) + ζ(a) N p + N n + N z n. N p = π(a), N n = ν(a), N z = ζ(a) ( ). P x = P y (Ax, x) = n µ j yj 2 p q 0, (Ax, x) = α 1 y α p y 2 p α p+1 y 2 p+1 α p+q y 2 p+q, α j > 0 (1 j p + q) (a) p {P y; y p+1 = = y n = 0} (Ax, x) > 0 p π(a). (b) q {P y; y 1 = = y p = y p+q+1 = = y n = 0} (Ax, x) < 0 q ν(a). (c) n (p+q) {Sy; y 1 = = y p+q = 0} (Ax, x) = 0 n (p+q) ζ(a). (c) p + q + ζ(a) n (a) π(a) p, (b) ν(a) q j=1 n = π(a) + ν(a) + ζ(a) p + q + ζ(a) n. p = π(a), q = ν(a). 33

34 [1] Nicholas J. Higham. Accuracy and Stability of Numerical Algorithms. SIAM, [2] [3].., 2003., I, II,, (1993, 1994)..., [4]. 1 (2013 ). jp/~mk/lecture/tahensuu1-2013/tahensuu1-new-text.pdf, [5].., [6].. 7., [7],,..., [8],.., [9] Lloyd Nicholas Trefethen and David Bau III. Numerical Linear Algebra. SIAM,