QR - PDF 無料ダウンロード

1 7 16 13 1 13.1 QR...................................... 2 13.1.1............................................ 2 13.1.2..................................... 3 13.1.3 QR........................................ 4 13.1.4 QR.................................... 4 13.2 QR................................. 6 13.3............................................. 7 13.3.1........................................... 8 13.3.2............................................ 8 13.3.3.................................... 8 13.3.4...................................... 9 13.3.5 QR................................ 9 13.3.6..................................... 9 13.4.................................................... 14 13 10 (x i, y i ), i = 1,..., n y i = y(x i ), i = 1,..., n y(x) y i y(x i ) y(x) ϕ i (x) y(x) = c 1 ϕ 1 (x) + c 2 ϕ(x) + + c m ϕ m (x), y i = y(x i ) = c 1 ϕ 1 (x i ) + c 2 ϕ(x i ) + + c m ϕ m (x i ), i = 1,..., n c i

13. 2 ϕ 1 (x 1 ) ϕ 2 (x 1 ) ϕ m (x 1 ) ϕ 1 (x 2 ) ϕ 2 (x 2 ) ϕ m (x 2 ).... ϕ 1 (x n ) ϕ 2 (x n ) ϕ m (x n ) c 1 c 2. c m = y 1 y 2. y n (1) n m A A = (a ij ) = (ϕ j (x i )) (2) c y c T = (c 1, c 2,, c m ) y T = (y 1, y 2,, y n ) (1) Ac = y (3) 13.1 QR 13.1.1 (3) m n c e = y Ac 2 = (y Ac) T (y Ac) (4) e = y T y c T A T y y T Ac c T A T Ac = y T y 2c T A T y + c T A T Ac e c T = 2AT y + 2A T Ac = 0 A T Ac = A T y A T A 0 c = (A T A) 1 A T y (3) A T A T Ac = A T y (5)

13. 3 13.1.2 1 2 // least squre method 3 4 #include "pch.h" 5 // #include <iostream> 6 #include <math.h> 7 8 #define N 9 9 #define M 5 10 11 12 double x[] = { 2.10000, 2.20000, 2.30000, 2.40000, 13 2.50000, 2.60000, 2.70000, 2.80000, 2.90000 }; 14 double y[] = { 0.17664, 0.27467, 0.46421, 0.83988, 15 1.54782, 2.64146, 3.78487, 4.47578, 4.63628 }; 16 17 int 18 main(int ac, char av[]) 19 { 20 double a[n][m], r[m][m]; 21 double b[m], c[m], s; 22 int n = N, m = M; 23 int i, j, k; 24 double p, q; 25 26 for (i = 0; i < n; i++) { 27 s = x[i]; 28 a[i][0] = 1; 29 a[i][1] = s; 30 s = s s; 31 a[i][2] = s; 32 s = s x[i]; 33 a[i][3] = s; 34 s = s x[i]; 35 a[i][4] = s; 36 } 37 38 for (j = 0; j < m; j++) { 39 for (k = 0; k < m; k++) { 40 s = 0; 41 for (i = 0; i < n; i++) 42 s = s + a[i][j] a[i][k]; 43 r[j][k] = s; 44 } 45 } 46 for (j = 0; j < m; j++) { // aˆt y > b 47 s = 0; 48 for (i = 0; i < n; i++) 49 s = s + a[i][j] y[i]; 50 b[j] = s; 51 } 52 for (i = 0; i < m; i++) { 53 p = r[i][i]; 54 for (j = i + 1; j < m; j++) { 55 q = r[j][i] / p; 56 for (k = i + 1; k < n; k++) 57 r[j][k] = r[j][k] q r[i][k];

13. 4 58 r[j][i] = 0; 59 b[j] = b[j] q b[i]; 60 } 61 } 62 63 for (i = m 1; i >= 0; i ) { 64 c[i] = b[i]; 65 for (j = m 1; j > i; j ) { 66 c[i] = c[i] r[i][j] c[j]; 67 } 68 c[i] = c[i] / r[i][i]; 69 } 70 printf("solution\n"); 71 for (i = 0; i < m; i++) { 72 printf("c%d: %18.15f, ", i + 1, c[i]); 73 } 74 printf("\n"); 75 } 13.1.3 QR (3) A QR A = QR Q n m (Q T Q = I) R m m Ac = QRc = y Q T Q T QRc = Rc = Q T y Rc = Q T y R 13.1.4 QR 1 // qr decomposition for non square matrix 2 3 #include "pch.h" 4 #include <iostream> 5 #include <math.h> 6 7 #define N 9 8 #define M 5 9 10 void gram schmidt(double a[n][m], int n, int m, double q[n][m]); 11 void qr decomposition(double a[n][m], int n, int m, 12 double q[n][m], double r[m][m]); 13 14 double x[] = { 2.10000, 2.20000, 2.30000, 2.40000, 15 2.50000, 2.60000, 2.70000, 2.80000, 2.90000 };

13. 5 16 double y[] = { 0.17664, 0.27467, 0.46421, 0.83988, 17 1.54782, 2.64146, 3.78487, 4.47578, 4.63628 }; 18 19 int 20 main(int ac, char av[]) 21 { 22 double a[n][m]; 23 double q[n][m], r[m][m]; 24 double b[m], c[m], s; 25 int n = N, m = M; 26 int i, j; 27 28 for (i = 0; i < n; i++) { 29 s = x[i]; 30 a[i][0] = 1; 31 a[i][1] = s; 32 s = s s; 33 a[i][2] = s; 34 s = s x[i]; 35 a[i][3] = s; 36 s = s x[i]; 37 a[i][4] = s; 38 } 39 40 qr decomposition(a, n, m, q, r); 41 for (j = 0; j < m; j++) { // QˆT y > b 42 s = 0; 43 for (i = 0; i < n; i++) 44 s = s + q[i][j] y[i]; 45 b[j] = s; 46 } 47 for (i = m 1; i >= 0; i ) { // solve Rc = b 48 c[i] = b[i]; 49 for (j = m 1; j > i; j ) { 50 c[i] = c[i] r[i][j] c[j]; 51 } 52 c[i] = c[i] / r[i][i]; 53 } 54 printf("solution\n"); 55 for (i = 0; i < m; i++) { 56 printf("c%d: %18.15f, ", i + 1, c[i]); 57 } 58 printf("\n"); 59 } 60 61 62 void 63 qr decomposition(double a[n][m], int n, int m, 64 double q[n][m], double r[m][m]) 65 { 66 int i, j, k, l; 67 double s, p; 68 69 gram schmidt(a, n, m, q); 70 for (k = 0; k < m; k++) { 71 for (j = 0; j < k; j++) { 72 s = 0; 73 for (i = 0; i < n; i++) 74 s = s + a[i][k] q[i][j]; 75 r[j][k] = s; 76 }

13. 6 77 s = 0; 78 for (i = 0; i < n; i++) { 79 p = a[i][k]; 80 for (l = 0; l < k; l++) { 81 p = p r[l][k] q[i][l]; 82 } 83 s = s + p p; 84 } 85 r[k][k] = sqrt(s); 86 for (j = k + 1; j < m; j++) 87 r[j][k] = 0; 88 } 89 } 90 91 void 92 gram schmidt(double a[n][m], int n, int m, double q[n][m]) 93 { 94 int i, j, k; 95 double s; 96 for (i = 0; i < n; i++) { 97 for (j = 0; j < m; j++) 98 q[i][j] = a[i][j]; 99 } 100 for (j = 0; j < m; j++) { 101 s = 0; 102 for (i = 0; i < n; i++) 103 s = s + q[i][j] q[i][j]; // (a[j], a[j]) 104 s = 1 / sqrt(s); 105 for (i = 0; i < n; i++) 106 q[i][j] = q[i][j] s; // a[j] 107 for (k = j + 1; k < m; k++) { 108 s = 0; 109 for (i = 0; i < n; i++) 110 s = s + q[i][k] q[i][j]; // (a[j], a[k]) 111 for (i = 0; i < n; i++) 112 q[i][k] = q[i][k] s q[i][j]; // a[i] 113 } 114 } 115 } 13.2 QR A = 1.00000 1.00000 0.00000 0.00001 0.00000 0.00000 a= 1.000000000000000, -1.000000000000000, 0.000000000000000, 0.000010000000000, 0.000000000000000, 0.000000000000000, a^ta= 1.000000000000000, -1.000000000000000, -1.000000000000000, 1.000000000100000,, B = 0.00000 0.00001 1.00000

13. 7 solution c1: 99999.991725963598583, c2: 99999.991725963598583, A T A c 1, c 2 A T A QR a= 1.000000000000000, -1.000000000000000, 0.000000000000000, 0.000010000000000, 0.000000000000000, 0.000000000000000, qr decomposition q= 1.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000, r= 1.000000000000000, -1.000000000000000, 0.000000000000000, 0.000010000000000, solution c1: 1.000000000000000, c2: 1.000000000000000, A A QR 13.3 n m A n m A = UΣV T U n m (U T U = I ) V m m Σ m m A 0 Σ σ 1, σ 2,..., σ m 0 σ 1 σ 2... σ m A n < m σ n+1,..., σ m 0 U 0 A A T A = (UΣV T ) T UΣV T = V Σ 2 V T Σ 2 V Σ 2 A T A {σ i 2, i = 1,..., m} A T A

13. 8 13.3.1 n m A A + (1) AA + A = A (2) A + AA + = A + (3) (AA + ) T = AA + (4) (A + A) T = A + A A + A (A + ) + = A (A T ) + = (A + ) T A B m (AB) + = B + A + A m A + = (A T A) 1 A T A A = UΣV T A + = V Σ + U T Σ + i σ i 0 1/σ i σ i 0 0 13.3.2 A n m Ax = b (6) k n x = A + b + (I A + A)k A + b x A A + = A 1 A + b (6) A x = A + b Ax b 2 13.3.3 1. A 2. U 3. QR 4. QR R k U V 5. V

13. 9 13.3.4 13.3.5 QR 13.3.6 1 #include <stdio.h> 2 #include <stdlib.h> 3 #include <math.h> 4 5 #define M 3 6 #define N 2 7 8 #define SIGN(a, b) ((b) >= 0.0? fabs(a) : fabs(a)) 9 #define MAX(x,y) ((x)>(y)?(x):(y)) 10 #define MIN(x,y) ((x)>(y)?(y):(x)) 11 12 double PYTHAG(double a, double b); 13 int svd(double a[m][n], int m, int n, double s[n], 14 double v[n][n]); 15 16 int 17 main(int ac, char av[]) 18 { 19 int n = N, m = M; 20 double a[m][n] = { {1, 1}, {0, 1.0e 5}, {0, 0} }; 21 double b[m] = { 0, 1.0e 5, 1 }; 22 double u[m][n], v[n][n], s[n]; 23 int i, j; 24 double x[m], t; 25 26 svd(a, m, n, s, v); 27 28 printf("a=\n"); 29 for (i = 0; i < m; i++) { 30 for (j = 0; j < n; j++) 31 printf("%10.5lf, ", a[i][j]); 32 printf("\n"); 33 } 34 printf("v=\n"); 35 for (i = 0; i < n; i++) { 36 for (j = 0; j < n; j++) 37 printf("%10.5lf, ", v[i][j]); 38 printf("\n"); 39 } 40 printf("s=\n"); 41 for (i = 0; i < n; i++) { 42 printf("%10.5lf", s[i]); 43 printf("\n"); 44 } 45 // solve Ax = B 46 for (i = 0; i < m; i++) { 47 t = 0; 48 for (j = 0; j < n; j++) 49 t = t + a[j][i] b[j]; 50 x[i] = s[i] == 0.0? 0.0 : t / s[i]; 51 x[i] = t; 52 } 53 for (i = 0; i < n; i++) {

13. 10 54 t = 0; 55 for (j = 0; j < n; j++) 56 t = t + v[j][i] x[i]; 57 printf("%18.15lf\n", t); 58 } 59 } 60 61 62 double 63 PYTHAG(double a, double b) 64 { 65 double at = fabs(a), bt = fabs(b), ct, result; 66 67 if (at > bt) { 68 ct = bt / at; 69 result = at sqrt(1.0 + ct ct); 70 } else if (bt > 0.0) { 71 ct = at / bt; 72 result = bt sqrt(1.0 + ct ct); 73 } else 74 result = 0.0; 75 return (result); 76 } 77 78 int 79 svd(double a[m][n], int m, int n, double w[n], double v[n][n]) 80 { 81 int flag, i, its, j, jj, k, l, nm; 82 double c, f, h, s, x, y, z; 83 double anorm = 0.0, g = 0.0, scale = 0.0; 84 double rv1[n]; 85 86 if (m < n) { 87 fprintf(stderr, "#rows must be > #cols \n"); 88 return (0); 89 } 90 91 // Householder reduction to bidiagonal form 92 for (i = 0; i < n; i++) { 93 // left hand reduction 94 l = i + 1; 95 rv1[i] = scale g; 96 g = s = scale = 0.0; 97 if (i < m) { 98 for (k = i; k < m; k++) 99 scale += fabs(a[k][i]); 100 if (scale) { 101 for (k = i; k < m; k++) { 102 a[k][i] = (a[k][i] / scale); 103 s += (a[k][i] a[k][i]); 104 } 105 f = a[i][i]; 106 g = SIGN(sqrt(s), f); 107 h = f g s; 108 a[i][i] = (f g); 109 if (i!= n 1) { 110 for (j = l; j < n; j++) { 111 for (s = 0.0, k = i; k < m; k++) 112 s += (a[k][i] a[k][j]); 113 f = s / h; 114 for (k = i; k < m; k++)

13. 11 115 a[k][j] += (f a[k][i]); 116 } 117 } 118 for (k = i; k < m; k++) 119 a[k][i] = (a[k][i] scale); 120 } 121 } 122 w[i] = (scale g); 123 // right hand reduction 124 g = s = scale = 0.0; 125 if (i < m && i!= n 1) { 126 for (k = l; k < n; k++) 127 scale += fabs(a[i][k]); 128 if (scale) { 129 for (k = l; k < n; k++) { 130 a[i][k] = (a[i][k] / scale); 131 s += (a[i][k] a[i][k]); 132 } 133 f = a[i][l]; 134 g = SIGN(sqrt(s), f); 135 h = f g s; 136 a[i][l] = (f g); 137 for (k = l; k < n; k++) 138 rv1[k] = a[i][k] / h; 139 if (i!= m 1) { 140 for (j = l; j < m; j++) { 141 for (s = 0.0, k = l; k < n; k++) 142 s += (a[j][k] a[i][k]); 143 for (k = l; k < n; k++) 144 a[j][k] += (s rv1[k]); 145 } 146 } 147 for (k = l; k < n; k++) 148 a[i][k] = (a[i][k] scale); 149 } 150 } 151 anorm = MAX(anorm, (fabs(w[i]) + fabs(rv1[i]))); 152 } 153 154 // accumulate the right hand transformation 155 for (i = n 1; i >= 0; i ) { 156 if (i < n 1) { 157 if (g) { 158 for (j = l; j < n; j++) 159 v[j][i] = ((a[i][j] / a[i][l]) / g); 160 for (j = l; j < n; j++) { 161 for (s = 0.0, k = l; k < n; k++) 162 s += (a[i][k] v[k][j]); 163 for (k = l; k < n; k++) 164 v[k][j] += (s v[k][i]); 165 } 166 } 167 for (j = l; j < n; j++) 168 v[i][j] = v[j][i] = 0.0; 169 } 170 v[i][i] = 1.0; 171 g = rv1[i]; 172 l = i; 173 } 174 175 // accumulate the left hand transformation

13. 12 176 for (i = n 1; i >= 0; i ) { 177 l = i + 1; 178 g = w[i]; 179 if (i < n 1) 180 for (j = l; j < n; j++) 181 a[i][j] = 0.0; 182 if (g) { 183 g = 1.0 / g; 184 if (i!= n 1) { 185 for (j = l; j < n; j++) { 186 for (s = 0.0, k = l; k < m; k++) 187 s += (a[k][i] a[k][j]); 188 f = (s / a[i][i]) g; 189 for (k = i; k < m; k++) 190 a[k][j] += (f a[k][i]); 191 } 192 } 193 for (j = i; j < m; j++) 194 a[j][i] = (a[j][i] g); 195 } else { 196 for (j = i; j < m; j++) 197 a[j][i] = 0.0; 198 } 199 ++a[i][i]; 200 } 201 202 // diagonalize the bidiagonal form 203 for (k = n 1; k >= 0; k ) { // loop over singular values 204 for (its = 0; its < 30; its++) { 205 // loop over allowed iterations 206 flag = 1; 207 for (l = k; l >= 0; l ) { // test for splitting 208 nm = l 1; 209 if (fabs(rv1[l]) + anorm == anorm) { 210 flag = 0; 211 break; 212 } 213 if (fabs(w[nm]) + anorm == anorm) 214 break; 215 } 216 if (flag) { 217 c = 0.0; 218 s = 1.0; 219 for (i = l; i <= k; i++) { 220 f = s rv1[i]; 221 if (fabs(f) + anorm!= anorm) { 222 g = w[i]; 223 h = PYTHAG(f, g); 224 w[i] = h; 225 h = 1.0 / h; 226 c = g h; 227 s = ( f h); 228 for (j = 0; j < m; j++) { 229 y = a[j][nm]; 230 z = a[j][i]; 231 a[j][nm] = (y c + z s); 232 a[j][i] = (z c y s); 233 } 234 } 235 } 236 }

13. 13 237 z = w[k]; 238 if (l == k) { // convergence 239 if (z < 0.0) { // make singular value nonnegative 240 w[k] = ( z); 241 for (j = 0; j < n; j++) 242 v[j][k] = ( v[j][k]); 243 } 244 break; 245 } 246 if (its >= 30) { 247 printf("no convergence after 30,000! iterations \n"); 248 return (0); 249 } 250 251 // shift from bottom 2 x 2 minor 252 x = w[l]; 253 nm = k 1; 254 y = w[nm]; 255 g = rv1[nm]; 256 h = rv1[k]; 257 f = ((y z) (y + z) + (g h) (g + h)) / (2.0 h y); 258 g = PYTHAG(f, 1.0); 259 f = ((x z) (x + z) + h ((y / (f + SIGN(g, f))) h)) / x; 260 261 // next QR transformation 262 c = s = 1.0; 263 for (j = l; j <= nm; j++) { 264 i = j + 1; 265 g = rv1[i]; 266 y = w[i]; 267 h = s g; 268 g = c g; 269 z = PYTHAG(f, h); 270 rv1[j] = z; 271 c = f / z; 272 s = h / z; 273 f = x c + g s; 274 g = g c x s; 275 h = y s; 276 y = y c; 277 for (jj = 0; jj < n; jj++) { 278 x = v[jj][j]; 279 z = v[jj][i]; 280 v[jj][j] = (x c + z s); 281 v[jj][i] = (z c x s); 282 } 283 z = PYTHAG(f, h); 284 w[j] = z; 285 if (z) { 286 z = 1.0 / z; 287 c = f z; 288 s = h z; 289 } 290 f = (c g) + (s y); 291 x = (c y) (s g); 292 for (jj = 0; jj < m; jj++) { 293 y = a[jj][j]; 294 z = a[jj][i]; 295 a[jj][j] = (y c + z s); 296 a[jj][i] = (z c y s); 297 }

13. 14 298 } 299 rv1[l] = 0.0; 300 rv1[k] = f; 301 w[k] = x; 302 } 303 } 304 return (1); 305 } 13.4 13.1 A = 10 10 10 10 9 10 10 1 0 0 0 1