Gauss

15 1 LU LDL T 6 : 1g00p013-5

1 6 1.1....................................... 7 1.2.................................. 8 1.3.................................. 8 2 Gauss 9 2.1..................................... 10 2.2 Gauss............................... 10 2.3.................................. 11 2.4 Gauss........................ 12 2.5...................................... 15 3 LU 16 3.1..................................... 17 3.2 LU................................ 17 3.3 LU........................... 21 3.4...................................... 23 4 Cholesky 24 4.1..................................... 25 4.2 Cholesky............................... 25 1

4.3 Cholesky.......................... 26 4.4...................................... 27 5 LDL T 28 5.1..................................... 29 5.2 LDL T.............................. 29 5.3................................ 31 5.4 LDL T......................... 32 5.5...................................... 35 6 36 6.1..................................... 37 6.2............................. 37 6.3......................... 37 6.4...................................... 39 7 40 7.1..................................... 41 7.2....................................... 41 7.3...................................... 42 8 43 8.1..................................... 44 8.2....................................... 44 8.3....................................... 44 8.4...................................... 45 46 48 2

A 49 3

6.1......................... 37 4

7.1 LU LDL T.................. 41 5

1 6

1.1 1 LDL T A A 1 1 LU LDL T 7

1.2 1 LU LDL T 1.3 2, Gauss 3,Gauss LU 4,LDL T Cholesky. 5,LDL T 6, 7, 8,. 8

2 Gauss 9

2.1 Gauss LU 2.2 Gauss x 1, x 2, x 3 a 11 x 1 + a 12 x 2 + a 13 x 3 = b 1 (2.1) a 21 x 1 + a 22 x 2 + a 23 x 3 = b 2 (2.2) a 31 x 1 + a 32 x 2 + a 33 x 3 = b 3 (2.3) 2x 1 + 5x 2 + 7x 3 = 23 (2.4) 4x 1 + 13x 2 + 20x 3 = 58 (2.5) 8x 1 + 29x 2 + 50x 3 = 132 (2.6) Gauss (2.4) (2.5) (2.6) 2x 1 + 5x 2 + 7x 3 = 23 (2.7) 13x 2 + 20x 3 = 12 (2.8) 29x 2 + 50x 3 = 40 (2.9) (2.8) (2.9). 2x 1 + 5x 2 + 7x 3 = 23 (2.10) 3x 2 + 6x 3 = 12 (2.11) 10 4x 3 = 4 (2.12)

Gauss (2.12) x 3 = 1 (2.11) x 2 = (12 6x 3 )/3 = (12 6 1)/3 = 2 (2.10) x 1 = (23 5x 2 7x 3 )/2 = (23 5 2 7 1)/2 = 3 a ii 0 2.3 Gauss a 11 0 a 21 x 1 + a 22 x 2 + a 23 x 3 = b 2 a 11 x 1 + a 12 x 2 + a 13 x 3 = b 1 a 31 x 1 + a 32 x 2 + a 33 x 3 = b 3 a 21 0 ( a 31 ) 11

0 a 11 0 x 1 a 12,,a 1n a 12,,a 1n a 11 a 12,,a 1n a 11 2.4 Gauss a 11 x 1 + a 12 x 2 + a 13 x 3 = b 1 a 21 x 1 + a 22 x 2 + a 23 x 3 = b 2 a 31 x 1 + a 32 x 2 + a 33 x 3 = b 3 k = 1, 2,..., n 1 k i (i = k + 1,..., n) (A b) = (A (1) b (1) ) (A (1) b (1) ) (A (2) b (2) ) (A (k) b (k) ) (A (k+1) b (k+1) ) (A (n) b (n) ) A (n) = u 11 u 12 u 1n u 22 u 2n.... 0 u nn, b(n) = y A (k) = (a (k) ij ), b (k) = (b (k) i ) Gauss 12

k = 1, 2,..., n 1 Gauss i = k + 1, k + 2,..., n m ik = a (k) ik /a(k) kk j = k + 1, k + 2,..., n b (k+1) i a (k+1) ij = b (k) i = a (k) ij m ik b (k) k m ik a (k) kj 1 a (1) 11 x 1 + a (1) 12 x 2 + + a (1) 1,n 1x n 1 + a (1) 1n x n = b (1) 1 a (2) 22 x 2 + + a (2) 2,n 1x n 1 + a (2) 2n x n = b (2) 2. a (n 1) n 1,n 1x n 1 + a (n 1) n 1,nx n a (n) nn x n = b (n 1) n 1 = b (n) n x n = b(n) n a (n) nn (2.13) 2 x n 1 x k = b(k) k a (k) k,k+1 x k+1 a (k) k,k+2 x k+2 a (k) kn x n a (k) kk (2.14) x 1 C 13

for (i=1;i<=n;i++){ pivot = a[i][i]; Gauss C for (j=i+1;j<=n;j++){ p=a[j][i]/pivot; a[j][i]=0.0; for (k=i+1;k<=n;k++) a[j][k] -= p*a[i][k]; b[j] -= p*b[i]; x[n] = b[n]/a[n][n]; for (i=n-1;i>=1;i--){ x[i]=b[i]; for (j=i+1;j<=n;j++) x[i] -= a[i][j]*x[j]; x[i] /= a[i][i]; / (n i) + (n i)(n i + 1) + (n i) + 1 = (n i)(n i + 3) + 1 / (n i)(n i + 1) + 1 + (n i 1) = (n i)(n i + 2) / n 1 1 + [(n i)(n i + 3) + 1] = 1 + (n 2 2ni + i 2 + 3n 3i + 1) i=1 n 1 i=1 14

n 1 = 1 + (n 2 n 1 + 3n + 1) 1 (2n + 3) i=1 = 1 + (n 2 + 3n 1)(n 1) (2n + 3) + (n 1)n(2n 1) 6 / n 1 i=1 (n i)(n i + 2) = n 1 i=1 = (n 2 + 2n) (n 2 2ni + i 2 + 2n 2i) n 1 i=1 n 1 1 (2n + 2) = (n 2 + 2n)(n 1) (2n + 2) = 2n3 + 3n 2 5n 6 i=1 i=1 = 2n3 + 6n 2 3n 6 i + n 1 i 2 i=1 (n 1)n 2 + i + n 1 i 2 i=1 (n 1)n 2 (n 1)n(2n 1) 6 n 3 /3 n 3 /3 1 Gauss O(n 3 ) 2.5 Gauss LU 15

3 LU 16

3.1 Gauss LU LU 3.2 LU Ax = b x = A 1 b Gauss b Gauss M 1 = 1 m 21 1 0 m 31 0 1...... m n1 0 0 1 A (1) = A M 1 A (1) = A (2) = a (1) 11 a (1) 12 a (1) 1n 0 a (2) 22 a (2) 2n...... 0 a (2) n2 a (2) nn 17

M 1 A (1) k 1 0 1.... 0 0 0 1 M K =, m jk = a (k) jk.. m k+1,k 1 /a(k) kk (3.1).. m k+2,k 0.......... 0 0 m n,k 0 1 k = 1, 2,..., n 1 Gauss i = k + 1, k + 2,..., n m ik = a (k) ik /a(k) kk j = k + 1, k + 2,..., n b (k+1) i a (k+1) ij = b (k) i = a (k) ij m ik b (k) k m ik a (k) kj M k A (k) = A (k+1) (3.2) M n 1 M n 2...M 2 M 1 A = A (n) 18

= U = a (1) 11 a (1) 12 a (1) 13 a (1) 1n a (2) 22 a (2) 23 a (2) 2n 0 a (1) 33 a (3) 3n.... a (n) nn (3.3) U 0 M n 1 M n 2...M 1 b = b (n) = b (1) 1 b (2) 2. b (n) n (3.4) Ax = b 3.3 U Ux = b (n) (3.5) M k 1 0 1.... 0 Mk 1 0 0 1 =.. m k+1,k 1.. m k+2,k 0.......... 0 0 m n,k 0 1, m jk = a (k) jk /a(k) kk (3.6) M k I L Mk 1 L = M 1 1 M 1 2...M 1 n 1 (3.7) 19

1 m 21 1 0 m 31 m 32 1 L =......... 1 m n1 m n2 m n3 m n,n 1 1, m jk = a (k) jk /a(k) kk (3.8) L 3.3 3.7 Gauss A A = LU (3.9) L U A LU A LU L,U Ax = b Ax = b LUx = b y Ly = b (3.10) Ux = y (3.11) x Ly = b m kk = 1 y 1 = b 1 k = 2, 3,..., n y k = b k k 1 j=1 m kj y j Ux = y 20

x n = y n /a (n) nn k = n 1, n 2,..., 1 n x k = y k a (k) kj x j j=k+1 /a (k) kk 3.3 LU 3 3 A A = 2 5 7 4 13 20 8 29 50 Gauss 2 5 7 0 3 6 0 0 4 0 0 2 5 7 A = 2 3 6 4 3 4 b A L U 21

A L U 2 5 7 1 0 0 2 5 7 4 13 20 = 2 1 0 0 3 6 8 29 50 4 3 1 0 0 4 (3.12) A L U LU LU for( i = 1; i <= n; i++ ){ for( j = 1; j <= n; j++ ){ LU C if( i == j ){ l[i][j] = 1; else{ l[i][j] = 0; u[i][j] = a[i][j]; for( i = 1; i <= n; i++ ){ for( j = i+1; j <= n; j++ ){ l[j][i] = u[j][i] / u[i][i]; for( k = i; k <= n; k++ ){ u[j][k] = u[j][k] - u[i][k] * l[j][i]; Gauss a[1][1] 0 a[1][1] a[2][1] a kk a[k][k] 0 0 a kk 22

3.4 Gauss LU Cholesky 23

4 Cholesky 24

4.1 LU A Cholesky 4.2 Cholesky A A = S T S 1 A LDL T Cholesky A > 0 a (k) kk D 1 2 = a (1) 11 0 a (2) 22... 0 a (n) nn (4.1) S = D 1 2 L T A A = S T S (4.2) A S s 11 s 12 s 1n s 22 s 2n S =.... 0 s nn S T S ij a ij (4.3) s 1i s ij + s 2i s 2j + + s ii s ij = a ij, (i j) (4.4) 25

s ii = a ii s ij = (a ij i 1 k=1 i 1 k=1 s 2 ki, (s 11 = a 11 ) s ki s kj )/s ii, (s 1j = a 1j /s 11 ) (4.5) s 11 = a 11 s 12, s 13,..., s 1n, s 22, s 23,... S S S T y = b, Sx = y Ax = b A 4.3 Cholesky Cholesky s ii = a ii s ij = (a ij i 1 k=1 i 1 k=1 s 2 ki, (s 11 = a 11 ) s ki s kj )/s ii, (s 1j = a 1j /s 11 ) (4.6) C 26

for (i = 1; i <= n; i++) { s = a[i][i]; for (k = 1; k < i; k++) s -= sqr(l[i][k]); if (s < 0) { Cholesky C exit(0); L[i][i] = sqrt(s); for (j =i + 1; j <= n; j++) { s = a[j][i]; for (k = 1; k < i; k++) s -= L[j][k] * L[i][k]; L[j][i] = s / L[i][i]; 4.4 A Cholesky LDL T 27

5 LDL T 28

5.1 A Cholesky LDL T 5.2 LDL T A LU U U U = DV (5.1) D = a (1) 11 0 a (2) 22 0 a (3) 33... a (n) nn (5.2) V = 1 v 12 v 13 v 1n 1 v 23 v 2n 1 v 3n 0.... 1, v kj = a (k) kj /a(k) kk, k < j (5.3) A A t = A (n k) (n k) Gauss 1 29

A (k) s = a (k) kk a (k) k+1,k. a (k) nk a (k) k,k+1 a (k) kn a (k) k+1,k+1 a (k) k+1,n..... a (k) n,k+1 a (k) nn (5.4) A (k) s = A (k)t s (5.5) a (k) ij = a (k) ji, k i, j n (5.6) k 1 a ij = a (k 1) ij a(k 1) i,k 1 a (k 1) k 1,k 1 a (k 1) = a (k 1) a (k) = a (k) αβ βα ij ji a (k 1) k 1,j (5.7) A = A T a (1) αβ = a(1) βα A v kj = a(k) kj a (k) kk = a(k) jk a (k) kk = m jk (5.8) V L V = L T (5.9) A A A = LDL T (5.10) A LDL T 30

5.3 A LDL T D ii d i = a (i) ii L ij l ij A = LDL T LDL T A k=1,2,...,n i l kj d j l ij = a ki i = 1, 2,..., k 1 j=1 k l ki d i l ki = a kk i=1 l ii = 1 d 1 = a 11 k = 2, 3,..., n LDL T i = 1, 2,..., k 1 i 1 l ki = a ki l kj l ij d j /d i d k = a kk j=1 k 1 lkid 2 i i=1 0 d 1 = a 11 l 21, d 2, l 31, l 32, d 3,... D L A LDL T Ax = LDL T x = Ly, y = DL T x Ly = b, DL T x = y Ax = b 31

5.4 LDL T LDL T 1 1 2x 1 + x 2 + x 3 = 8 (5.11) x 1 + 3x 2 + 2x 3 = 11 (5.12) x 1 + 2x 2 + 4x 3 = 16 (5.13) 2 1 1 x 1 1 3 2 = x 2 1 2 4 x 3 8 11 16 (5.14) Ax = b (5.15) A LDL T A = LDL T (5.16) 1 d 11 1 l 21 l 31 = l 21 1 d 22 1 l 32 (5.17) l 31 l 32 1 d 33 1 d 11 d 11 l 21 d 11 l 31 = d 11 d 11 l21 2 + d 22 d 11 l 21 l 31 + d 22 l 32 (5.18) d 11 d 11 l 21 l 31 + d 22 l 32 d 11 l31 2 + d 22 l32 2 + d 33 1 2 1 0.5 0.5 = 0.5 1 2.5 1 0.6 (5.19) 0.5 0.6 1 2.6 1 32

LU AX = b LDL T x = b LDL T x = b (5.20) DL T x = y (5.21) Ly = b (5.22) 1 0.5 1 0.5 0.6 1 = y 1 y 2 y 3 8 11 16 (5.23) y 1 y 2 y 3 = 8 7 7.8 (5.24) DL T x = y 2 2.5 2.6 1 0.5 0.5 1 0.6 1 x 1 x 2 x 3 = 8 7 7.8 (5.25) 2 1 1 2.5 1.5 2.6 x 1 x 2 x 3 = 8 7 7.8 (5.26) 33

x 1 x 2 x 3 = 2 1 3 (5.27) L D d 11 = a 11, l j1 = a j1 /d 1 (j = 2,..., n) d 22 = a 22 d 11 l 2 21, l j2 = (a 2j d 11 l 21 l j1 /d 22 )(j = 3,..., n) d 33 = a 33 d 11 l 2 31 d 22 l 2 32, l j3 = (a 3j d 11 l 31 l j1 d 2 l 32 l j2 )/d 33 (j = 4,..., n) d[1]=a[1][1]; for(k=2;k<=n;k++){ s=0;t=0; for(i=1;i<=k-1;i++){ LDL T C for(j=1;j<=i-1;j++){ s+=l[k][j]*l[i][j]*d[j]; l[k][i]=(a[k][i] - s)/d[i]; for(i=1;i<=k-1;i++){ t+=l[k][i]*l[k][i]*d[i]; d[k]=a[k][k]-t ; 34

5.5 A LDL T LU LDL T 35

6 36

6.1 C LU LDL T. 6.2,,. CPU memory OS Celeron 2.0GHz 512MB Windows XP cygwin+gcc 6.1: 6.3 LU 1. LDL T A 2. A a 11 2,3,...,n 37

A = 1 2 3 n 2 2 3 n 3 3 3 n...... n n n n n n A 3. A x x i = i, (i = 1,..., n) 4. Ax = b b 5. A LU 6. Ly = b y 7. Ux = y x 8. 5 7 LDL T 1. A, x, b LU 2. A LDL T 3. Ly = b y 4. DL T x = y x 5. 2 4 38

6.4 C LU LDL T 39

7 40

7.1 C LU LDL T 7.2 100 1000 0.006 "LU.txt" "LDLT.txt" 0.005 0.004 time 0.003 0.002 0.001 0 100 200 300 400 500 600 700 800 900 1000 dimension 7.1: LU LDL T 100,LDL T 0.000008 LU 0.000015 1000 LU n 3 /3 + O(n 2 ) LDL T n 3 /6 + O(n 2 ) 41

7.3 C LU LDL T 42

8 43

8.1 8.2 2, Gauss 3,Gauss LU 4,LDL T Cholesky. 5,LDL T 6, 7, 8.3 A LDL T LU 44

8.4 LU LDL T LU n 3 /3 + O(n 2 ) LDL T n 3 /6 + O(n 2 ) LDL T LU 45

46

,,..,,. 47

[1] :,, 2003, [2] :, 1999 [3] :,, 2002, [4] G.H.Golub and C.F.Van Loan: Matrix Computations,Johns Hopkins,1996 48

Appendix A 49

///////////////LDL^T_Factorization/////////////// #include <stdio.h> #include <stdlib.h> #include <time.h> #define N 1000 /* */ int main(void); int main(void){ int i,j,k; double s; double t; double sum; static double a[n+1][n+1]; static double b[n+1]; static double w[n+1][n+1]; static double x[n+1]; static double y[n+1]; static double l[n+1][n+1]; static double d[n+1]; static double ld[n+1][n+1]; unsigned long t1,t2 ; /*a_{11 a[1][1] */ /*A */ for(i=1;i<=n;i++){ for(j=i;j<=n;j++){ a[i][j]=j; for(j=i-1;j>=1;j--){ a[i][j]=a[j][i]; /* X */ 50

for(i=1;i<=n;i++){ x[i]=i; /*Ax=b B */ for(i=1;i<=n;i++){ for(k=1,sum=0;k<=n;k++){ sum+=a[i][k]*x[k]; b[i] = sum; /* */ for(i=1;i<=n;i++){ x[i]=0; t1=clock(); /* LDL^T */ /*L */ for(i=1;i<=n;i++){ for(j=1;j<=n;j++){ if(i==j){ else{ l[i][j]=1; l[i][j]=0; d[1]=a[1][1]; for(k=2;k<=n;k++){ for(i=1;i<=k-1;i++){ for(j=1,s=0;j<=i-1;j++){ 51

s+=l[k][j]*l[i][j]*d[j]; l[k][i]=(a[k][i] - s)/d[i]; for(i=1,t=0;i<=k-1;i++){ t+=l[k][i]*l[k][i]*d[i]; d[k]=a[k][k]-t ; /* LDL^T */ for(i=1;i<=n;i++){ y[i]=b[i]; for(j=1;j<i;j++){ y[i]-=l[i][j]*y[j]; /* DL^T*x=y */ /*DL^T */ for(i=1;i<=n;i++){ for(j=1;j<=n;j++){ for(i=n;i>=1;i--){ w[i][j]=d[i]*l[j][i]; x[i]=y[i]; for(j=n;j>i;j--){ x[i]-=w[i][j]*x[j]; x[i]/=w[i][i]; t2=clock(); /*l[j][i] ^ */ 52

printf(" ; %f[ms]\n", (t2-t1)/1000000.0 ); return(0); ///////////////LU_Factorization/////////////// #include <stdio.h> #include <stdlib.h> #include <time.h> #define N 1000 int main(void); int main(void){ int i,j,k; double s,t,sum; static double a[n+1][n+1]; static double b[n+1]; static double u[n+1][n+1]; static double w[n+1][n+1]; static double x[n+1]; static double y[n+1]; static double l[n+1][n+1]; static double d[n+1]; static double ld[n+1][n+1]; unsigned long t1,t2 ; /*a_{11 a[1][1] */ /*A */ for(i=1;i<=n;i++){ for(j=i;j<=n;j++){ a[i][j]=j; for(j=i-1;j>=1;j--){ a[i][j]=a[j][i]; 53

/* X */ for(i=1;i<=n;i++){ x[i]=i; /*Ax=b B */ for(i=1;i<=n;i++){ for(k=1,sum=0;k<=n;k++){ sum+=a[i][k]*x[k]; b[i] = sum; /* */ for(i=1;i<=n;i++){ x[i]=0; t1=clock(); /* LU */ for( i = 1; i <= N; i++ ){ for( j = 1; j <= N; j++ ){ if( i == j ){ l[i][j] = 1; else{ l[i][j] = 0; u[i][j] = a[i][j]; for( i = 1; i <= N; i++ ){ for( j = i+1; j <= N; j++ ){ l[j][i] = u[j][i] / u[i][i]; for( k = i; k <= N; k++ ){ u[j][k] = u[j][k] - u[i][k] * l[j][i]; 54

for(i=1;i<=n;i++){ y[i]=b[i]; for(j=1;j<i;j++){ y[i]-=l[i][j]*y[j]; for(i=n;i>=1;i--){ x[i]=y[i]; for(j=n;j>i;j--){ x[i]-=u[i][j]*x[j]; x[i]/=u[i][i]; t2=clock(); printf(" ; %f[ms]\n", (t2-t1)/1000000.0 ); return(0); 55