| [ Index ] |
PHP Cross Reference of moodle-2.8 |
[Summary view] [Print] [Text view]
1 <?php 2 /** 3 * @package JAMA 4 * 5 * For an m-by-n matrix A with m >= n, the singular value decomposition is 6 * an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and 7 * an n-by-n orthogonal matrix V so that A = U*S*V'. 8 * 9 * The singular values, sigma[$k] = S[$k][$k], are ordered so that 10 * sigma[0] >= sigma[1] >= ... >= sigma[n-1]. 11 * 12 * The singular value decompostion always exists, so the constructor will 13 * never fail. The matrix condition number and the effective numerical 14 * rank can be computed from this decomposition. 15 * 16 * @author Paul Meagher 17 * @license PHP v3.0 18 * @version 1.1 19 */ 20 class SingularValueDecomposition { 21 22 /** 23 * Internal storage of U. 24 * @var array 25 */ 26 private $U = array(); 27 28 /** 29 * Internal storage of V. 30 * @var array 31 */ 32 private $V = array(); 33 34 /** 35 * Internal storage of singular values. 36 * @var array 37 */ 38 private $s = array(); 39 40 /** 41 * Row dimension. 42 * @var int 43 */ 44 private $m; 45 46 /** 47 * Column dimension. 48 * @var int 49 */ 50 private $n; 51 52 53 /** 54 * Construct the singular value decomposition 55 * 56 * Derived from LINPACK code. 57 * 58 * @param $A Rectangular matrix 59 * @return Structure to access U, S and V. 60 */ 61 public function __construct($Arg) { 62 63 // Initialize. 64 $A = $Arg->getArrayCopy(); 65 $this->m = $Arg->getRowDimension(); 66 $this->n = $Arg->getColumnDimension(); 67 $nu = min($this->m, $this->n); 68 $e = array(); 69 $work = array(); 70 $wantu = true; 71 $wantv = true; 72 $nct = min($this->m - 1, $this->n); 73 $nrt = max(0, min($this->n - 2, $this->m)); 74 75 // Reduce A to bidiagonal form, storing the diagonal elements 76 // in s and the super-diagonal elements in e. 77 for ($k = 0; $k < max($nct,$nrt); ++$k) { 78 79 if ($k < $nct) { 80 // Compute the transformation for the k-th column and 81 // place the k-th diagonal in s[$k]. 82 // Compute 2-norm of k-th column without under/overflow. 83 $this->s[$k] = 0; 84 for ($i = $k; $i < $this->m; ++$i) { 85 $this->s[$k] = hypo($this->s[$k], $A[$i][$k]); 86 } 87 if ($this->s[$k] != 0.0) { 88 if ($A[$k][$k] < 0.0) { 89 $this->s[$k] = -$this->s[$k]; 90 } 91 for ($i = $k; $i < $this->m; ++$i) { 92 $A[$i][$k] /= $this->s[$k]; 93 } 94 $A[$k][$k] += 1.0; 95 } 96 $this->s[$k] = -$this->s[$k]; 97 } 98 99 for ($j = $k + 1; $j < $this->n; ++$j) { 100 if (($k < $nct) & ($this->s[$k] != 0.0)) { 101 // Apply the transformation. 102 $t = 0; 103 for ($i = $k; $i < $this->m; ++$i) { 104 $t += $A[$i][$k] * $A[$i][$j]; 105 } 106 $t = -$t / $A[$k][$k]; 107 for ($i = $k; $i < $this->m; ++$i) { 108 $A[$i][$j] += $t * $A[$i][$k]; 109 } 110 // Place the k-th row of A into e for the 111 // subsequent calculation of the row transformation. 112 $e[$j] = $A[$k][$j]; 113 } 114 } 115 116 if ($wantu AND ($k < $nct)) { 117 // Place the transformation in U for subsequent back 118 // multiplication. 119 for ($i = $k; $i < $this->m; ++$i) { 120 $this->U[$i][$k] = $A[$i][$k]; 121 } 122 } 123 124 if ($k < $nrt) { 125 // Compute the k-th row transformation and place the 126 // k-th super-diagonal in e[$k]. 127 // Compute 2-norm without under/overflow. 128 $e[$k] = 0; 129 for ($i = $k + 1; $i < $this->n; ++$i) { 130 $e[$k] = hypo($e[$k], $e[$i]); 131 } 132 if ($e[$k] != 0.0) { 133 if ($e[$k+1] < 0.0) { 134 $e[$k] = -$e[$k]; 135 } 136 for ($i = $k + 1; $i < $this->n; ++$i) { 137 $e[$i] /= $e[$k]; 138 } 139 $e[$k+1] += 1.0; 140 } 141 $e[$k] = -$e[$k]; 142 if (($k+1 < $this->m) AND ($e[$k] != 0.0)) { 143 // Apply the transformation. 144 for ($i = $k+1; $i < $this->m; ++$i) { 145 $work[$i] = 0.0; 146 } 147 for ($j = $k+1; $j < $this->n; ++$j) { 148 for ($i = $k+1; $i < $this->m; ++$i) { 149 $work[$i] += $e[$j] * $A[$i][$j]; 150 } 151 } 152 for ($j = $k + 1; $j < $this->n; ++$j) { 153 $t = -$e[$j] / $e[$k+1]; 154 for ($i = $k + 1; $i < $this->m; ++$i) { 155 $A[$i][$j] += $t * $work[$i]; 156 } 157 } 158 } 159 if ($wantv) { 160 // Place the transformation in V for subsequent 161 // back multiplication. 162 for ($i = $k + 1; $i < $this->n; ++$i) { 163 $this->V[$i][$k] = $e[$i]; 164 } 165 } 166 } 167 } 168 169 // Set up the final bidiagonal matrix or order p. 170 $p = min($this->n, $this->m + 1); 171 if ($nct < $this->n) { 172 $this->s[$nct] = $A[$nct][$nct]; 173 } 174 if ($this->m < $p) { 175 $this->s[$p-1] = 0.0; 176 } 177 if ($nrt + 1 < $p) { 178 $e[$nrt] = $A[$nrt][$p-1]; 179 } 180 $e[$p-1] = 0.0; 181 // If required, generate U. 182 if ($wantu) { 183 for ($j = $nct; $j < $nu; ++$j) { 184 for ($i = 0; $i < $this->m; ++$i) { 185 $this->U[$i][$j] = 0.0; 186 } 187 $this->U[$j][$j] = 1.0; 188 } 189 for ($k = $nct - 1; $k >= 0; --$k) { 190 if ($this->s[$k] != 0.0) { 191 for ($j = $k + 1; $j < $nu; ++$j) { 192 $t = 0; 193 for ($i = $k; $i < $this->m; ++$i) { 194 $t += $this->U[$i][$k] * $this->U[$i][$j]; 195 } 196 $t = -$t / $this->U[$k][$k]; 197 for ($i = $k; $i < $this->m; ++$i) { 198 $this->U[$i][$j] += $t * $this->U[$i][$k]; 199 } 200 } 201 for ($i = $k; $i < $this->m; ++$i ) { 202 $this->U[$i][$k] = -$this->U[$i][$k]; 203 } 204 $this->U[$k][$k] = 1.0 + $this->U[$k][$k]; 205 for ($i = 0; $i < $k - 1; ++$i) { 206 $this->U[$i][$k] = 0.0; 207 } 208 } else { 209 for ($i = 0; $i < $this->m; ++$i) { 210 $this->U[$i][$k] = 0.0; 211 } 212 $this->U[$k][$k] = 1.0; 213 } 214 } 215 } 216 217 // If required, generate V. 218 if ($wantv) { 219 for ($k = $this->n - 1; $k >= 0; --$k) { 220 if (($k < $nrt) AND ($e[$k] != 0.0)) { 221 for ($j = $k + 1; $j < $nu; ++$j) { 222 $t = 0; 223 for ($i = $k + 1; $i < $this->n; ++$i) { 224 $t += $this->V[$i][$k]* $this->V[$i][$j]; 225 } 226 $t = -$t / $this->V[$k+1][$k]; 227 for ($i = $k + 1; $i < $this->n; ++$i) { 228 $this->V[$i][$j] += $t * $this->V[$i][$k]; 229 } 230 } 231 } 232 for ($i = 0; $i < $this->n; ++$i) { 233 $this->V[$i][$k] = 0.0; 234 } 235 $this->V[$k][$k] = 1.0; 236 } 237 } 238 239 // Main iteration loop for the singular values. 240 $pp = $p - 1; 241 $iter = 0; 242 $eps = pow(2.0, -52.0); 243 244 while ($p > 0) { 245 // Here is where a test for too many iterations would go. 246 // This section of the program inspects for negligible 247 // elements in the s and e arrays. On completion the 248 // variables kase and k are set as follows: 249 // kase = 1 if s(p) and e[k-1] are negligible and k<p 250 // kase = 2 if s(k) is negligible and k<p 251 // kase = 3 if e[k-1] is negligible, k<p, and 252 // s(k), ..., s(p) are not negligible (qr step). 253 // kase = 4 if e(p-1) is negligible (convergence). 254 for ($k = $p - 2; $k >= -1; --$k) { 255 if ($k == -1) { 256 break; 257 } 258 if (abs($e[$k]) <= $eps * (abs($this->s[$k]) + abs($this->s[$k+1]))) { 259 $e[$k] = 0.0; 260 break; 261 } 262 } 263 if ($k == $p - 2) { 264 $kase = 4; 265 } else { 266 for ($ks = $p - 1; $ks >= $k; --$ks) { 267 if ($ks == $k) { 268 break; 269 } 270 $t = ($ks != $p ? abs($e[$ks]) : 0.) + ($ks != $k + 1 ? abs($e[$ks-1]) : 0.); 271 if (abs($this->s[$ks]) <= $eps * $t) { 272 $this->s[$ks] = 0.0; 273 break; 274 } 275 } 276 if ($ks == $k) { 277 $kase = 3; 278 } else if ($ks == $p-1) { 279 $kase = 1; 280 } else { 281 $kase = 2; 282 $k = $ks; 283 } 284 } 285 ++$k; 286 287 // Perform the task indicated by kase. 288 switch ($kase) { 289 // Deflate negligible s(p). 290 case 1: 291 $f = $e[$p-2]; 292 $e[$p-2] = 0.0; 293 for ($j = $p - 2; $j >= $k; --$j) { 294 $t = hypo($this->s[$j],$f); 295 $cs = $this->s[$j] / $t; 296 $sn = $f / $t; 297 $this->s[$j] = $t; 298 if ($j != $k) { 299 $f = -$sn * $e[$j-1]; 300 $e[$j-1] = $cs * $e[$j-1]; 301 } 302 if ($wantv) { 303 for ($i = 0; $i < $this->n; ++$i) { 304 $t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$p-1]; 305 $this->V[$i][$p-1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$p-1]; 306 $this->V[$i][$j] = $t; 307 } 308 } 309 } 310 break; 311 // Split at negligible s(k). 312 case 2: 313 $f = $e[$k-1]; 314 $e[$k-1] = 0.0; 315 for ($j = $k; $j < $p; ++$j) { 316 $t = hypo($this->s[$j], $f); 317 $cs = $this->s[$j] / $t; 318 $sn = $f / $t; 319 $this->s[$j] = $t; 320 $f = -$sn * $e[$j]; 321 $e[$j] = $cs * $e[$j]; 322 if ($wantu) { 323 for ($i = 0; $i < $this->m; ++$i) { 324 $t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$k-1]; 325 $this->U[$i][$k-1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$k-1]; 326 $this->U[$i][$j] = $t; 327 } 328 } 329 } 330 break; 331 // Perform one qr step. 332 case 3: 333 // Calculate the shift. 334 $scale = max(max(max(max( 335 abs($this->s[$p-1]),abs($this->s[$p-2])),abs($e[$p-2])), 336 abs($this->s[$k])), abs($e[$k])); 337 $sp = $this->s[$p-1] / $scale; 338 $spm1 = $this->s[$p-2] / $scale; 339 $epm1 = $e[$p-2] / $scale; 340 $sk = $this->s[$k] / $scale; 341 $ek = $e[$k] / $scale; 342 $b = (($spm1 + $sp) * ($spm1 - $sp) + $epm1 * $epm1) / 2.0; 343 $c = ($sp * $epm1) * ($sp * $epm1); 344 $shift = 0.0; 345 if (($b != 0.0) || ($c != 0.0)) { 346 $shift = sqrt($b * $b + $c); 347 if ($b < 0.0) { 348 $shift = -$shift; 349 } 350 $shift = $c / ($b + $shift); 351 } 352 $f = ($sk + $sp) * ($sk - $sp) + $shift; 353 $g = $sk * $ek; 354 // Chase zeros. 355 for ($j = $k; $j < $p-1; ++$j) { 356 $t = hypo($f,$g); 357 $cs = $f/$t; 358 $sn = $g/$t; 359 if ($j != $k) { 360 $e[$j-1] = $t; 361 } 362 $f = $cs * $this->s[$j] + $sn * $e[$j]; 363 $e[$j] = $cs * $e[$j] - $sn * $this->s[$j]; 364 $g = $sn * $this->s[$j+1]; 365 $this->s[$j+1] = $cs * $this->s[$j+1]; 366 if ($wantv) { 367 for ($i = 0; $i < $this->n; ++$i) { 368 $t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$j+1]; 369 $this->V[$i][$j+1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$j+1]; 370 $this->V[$i][$j] = $t; 371 } 372 } 373 $t = hypo($f,$g); 374 $cs = $f/$t; 375 $sn = $g/$t; 376 $this->s[$j] = $t; 377 $f = $cs * $e[$j] + $sn * $this->s[$j+1]; 378 $this->s[$j+1] = -$sn * $e[$j] + $cs * $this->s[$j+1]; 379 $g = $sn * $e[$j+1]; 380 $e[$j+1] = $cs * $e[$j+1]; 381 if ($wantu && ($j < $this->m - 1)) { 382 for ($i = 0; $i < $this->m; ++$i) { 383 $t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$j+1]; 384 $this->U[$i][$j+1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$j+1]; 385 $this->U[$i][$j] = $t; 386 } 387 } 388 } 389 $e[$p-2] = $f; 390 $iter = $iter + 1; 391 break; 392 // Convergence. 393 case 4: 394 // Make the singular values positive. 395 if ($this->s[$k] <= 0.0) { 396 $this->s[$k] = ($this->s[$k] < 0.0 ? -$this->s[$k] : 0.0); 397 if ($wantv) { 398 for ($i = 0; $i <= $pp; ++$i) { 399 $this->V[$i][$k] = -$this->V[$i][$k]; 400 } 401 } 402 } 403 // Order the singular values. 404 while ($k < $pp) { 405 if ($this->s[$k] >= $this->s[$k+1]) { 406 break; 407 } 408 $t = $this->s[$k]; 409 $this->s[$k] = $this->s[$k+1]; 410 $this->s[$k+1] = $t; 411 if ($wantv AND ($k < $this->n - 1)) { 412 for ($i = 0; $i < $this->n; ++$i) { 413 $t = $this->V[$i][$k+1]; 414 $this->V[$i][$k+1] = $this->V[$i][$k]; 415 $this->V[$i][$k] = $t; 416 } 417 } 418 if ($wantu AND ($k < $this->m-1)) { 419 for ($i = 0; $i < $this->m; ++$i) { 420 $t = $this->U[$i][$k+1]; 421 $this->U[$i][$k+1] = $this->U[$i][$k]; 422 $this->U[$i][$k] = $t; 423 } 424 } 425 ++$k; 426 } 427 $iter = 0; 428 --$p; 429 break; 430 } // end switch 431 } // end while 432 433 } // end constructor 434 435 436 /** 437 * Return the left singular vectors 438 * 439 * @access public 440 * @return U 441 */ 442 public function getU() { 443 return new Matrix($this->U, $this->m, min($this->m + 1, $this->n)); 444 } 445 446 447 /** 448 * Return the right singular vectors 449 * 450 * @access public 451 * @return V 452 */ 453 public function getV() { 454 return new Matrix($this->V); 455 } 456 457 458 /** 459 * Return the one-dimensional array of singular values 460 * 461 * @access public 462 * @return diagonal of S. 463 */ 464 public function getSingularValues() { 465 return $this->s; 466 } 467 468 469 /** 470 * Return the diagonal matrix of singular values 471 * 472 * @access public 473 * @return S 474 */ 475 public function getS() { 476 for ($i = 0; $i < $this->n; ++$i) { 477 for ($j = 0; $j < $this->n; ++$j) { 478 $S[$i][$j] = 0.0; 479 } 480 $S[$i][$i] = $this->s[$i]; 481 } 482 return new Matrix($S); 483 } 484 485 486 /** 487 * Two norm 488 * 489 * @access public 490 * @return max(S) 491 */ 492 public function norm2() { 493 return $this->s[0]; 494 } 495 496 497 /** 498 * Two norm condition number 499 * 500 * @access public 501 * @return max(S)/min(S) 502 */ 503 public function cond() { 504 return $this->s[0] / $this->s[min($this->m, $this->n) - 1]; 505 } 506 507 508 /** 509 * Effective numerical matrix rank 510 * 511 * @access public 512 * @return Number of nonnegligible singular values. 513 */ 514 public function rank() { 515 $eps = pow(2.0, -52.0); 516 $tol = max($this->m, $this->n) * $this->s[0] * $eps; 517 $r = 0; 518 for ($i = 0; $i < count($this->s); ++$i) { 519 if ($this->s[$i] > $tol) { 520 ++$r; 521 } 522 } 523 return $r; 524 } 525 526 } // class SingularValueDecomposition
title
Description
Body
title
Description
Body
title
Description
Body
title
Body
| Generated: Fri Nov 28 20:29:05 2014 | Cross-referenced by PHPXref 0.7.1 |