Rev 2388 | Blame | Compare with Previous | Last modification | View Log | RSS feed
<?php/*** @package JAMA** For an m-by-n matrix A with m >= n, the singular value decomposition is* an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and* an n-by-n orthogonal matrix V so that A = U*S*V'.** The singular values, sigma[$k] = S[$k][$k], are ordered so that* sigma[0] >= sigma[1] >= ... >= sigma[n-1].** The singular value decompostion always exists, so the constructor will* never fail. The matrix condition number and the effective numerical* rank can be computed from this decomposition.** @author Paul Meagher* @license PHP v3.0* @version 1.1*/class SingularValueDecomposition {/*** Internal storage of U.* @var array*/private $U = array();/*** Internal storage of V.* @var array*/private $V = array();/*** Internal storage of singular values.* @var array*/private $s = array();/*** Row dimension.* @var int*/private $m;/*** Column dimension.* @var int*/private $n;/*** Construct the singular value decomposition** Derived from LINPACK code.** @param $A Rectangular matrix* @return Structure to access U, S and V.*/public function __construct($Arg) {// Initialize.$A = $Arg->getArrayCopy();$this->m = $Arg->getRowDimension();$this->n = $Arg->getColumnDimension();$nu = min($this->m, $this->n);$e = array();$work = array();$wantu = true;$wantv = true;$nct = min($this->m - 1, $this->n);$nrt = max(0, min($this->n - 2, $this->m));// Reduce A to bidiagonal form, storing the diagonal elements// in s and the super-diagonal elements in e.for ($k = 0; $k < max($nct,$nrt); ++$k) {if ($k < $nct) {// Compute the transformation for the k-th column and// place the k-th diagonal in s[$k].// Compute 2-norm of k-th column without under/overflow.$this->s[$k] = 0;for ($i = $k; $i < $this->m; ++$i) {$this->s[$k] = hypo($this->s[$k], $A[$i][$k]);}if ($this->s[$k] != 0.0) {if ($A[$k][$k] < 0.0) {$this->s[$k] = -$this->s[$k];}for ($i = $k; $i < $this->m; ++$i) {$A[$i][$k] /= $this->s[$k];}$A[$k][$k] += 1.0;}$this->s[$k] = -$this->s[$k];}for ($j = $k + 1; $j < $this->n; ++$j) {if (($k < $nct) & ($this->s[$k] != 0.0)) {// Apply the transformation.$t = 0;for ($i = $k; $i < $this->m; ++$i) {$t += $A[$i][$k] * $A[$i][$j];}$t = -$t / $A[$k][$k];for ($i = $k; $i < $this->m; ++$i) {$A[$i][$j] += $t * $A[$i][$k];}// Place the k-th row of A into e for the// subsequent calculation of the row transformation.$e[$j] = $A[$k][$j];}}if ($wantu AND ($k < $nct)) {// Place the transformation in U for subsequent back// multiplication.for ($i = $k; $i < $this->m; ++$i) {$this->U[$i][$k] = $A[$i][$k];}}if ($k < $nrt) {// Compute the k-th row transformation and place the// k-th super-diagonal in e[$k].// Compute 2-norm without under/overflow.$e[$k] = 0;for ($i = $k + 1; $i < $this->n; ++$i) {$e[$k] = hypo($e[$k], $e[$i]);}if ($e[$k] != 0.0) {if ($e[$k+1] < 0.0) {$e[$k] = -$e[$k];}for ($i = $k + 1; $i < $this->n; ++$i) {$e[$i] /= $e[$k];}$e[$k+1] += 1.0;}$e[$k] = -$e[$k];if (($k+1 < $this->m) AND ($e[$k] != 0.0)) {// Apply the transformation.for ($i = $k+1; $i < $this->m; ++$i) {$work[$i] = 0.0;}for ($j = $k+1; $j < $this->n; ++$j) {for ($i = $k+1; $i < $this->m; ++$i) {$work[$i] += $e[$j] * $A[$i][$j];}}for ($j = $k + 1; $j < $this->n; ++$j) {$t = -$e[$j] / $e[$k+1];for ($i = $k + 1; $i < $this->m; ++$i) {$A[$i][$j] += $t * $work[$i];}}}if ($wantv) {// Place the transformation in V for subsequent// back multiplication.for ($i = $k + 1; $i < $this->n; ++$i) {$this->V[$i][$k] = $e[$i];}}}}// Set up the final bidiagonal matrix or order p.$p = min($this->n, $this->m + 1);if ($nct < $this->n) {$this->s[$nct] = $A[$nct][$nct];}if ($this->m < $p) {$this->s[$p-1] = 0.0;}if ($nrt + 1 < $p) {$e[$nrt] = $A[$nrt][$p-1];}$e[$p-1] = 0.0;// If required, generate U.if ($wantu) {for ($j = $nct; $j < $nu; ++$j) {for ($i = 0; $i < $this->m; ++$i) {$this->U[$i][$j] = 0.0;}$this->U[$j][$j] = 1.0;}for ($k = $nct - 1; $k >= 0; --$k) {if ($this->s[$k] != 0.0) {for ($j = $k + 1; $j < $nu; ++$j) {$t = 0;for ($i = $k; $i < $this->m; ++$i) {$t += $this->U[$i][$k] * $this->U[$i][$j];}$t = -$t / $this->U[$k][$k];for ($i = $k; $i < $this->m; ++$i) {$this->U[$i][$j] += $t * $this->U[$i][$k];}}for ($i = $k; $i < $this->m; ++$i ) {$this->U[$i][$k] = -$this->U[$i][$k];}$this->U[$k][$k] = 1.0 + $this->U[$k][$k];for ($i = 0; $i < $k - 1; ++$i) {$this->U[$i][$k] = 0.0;}} else {for ($i = 0; $i < $this->m; ++$i) {$this->U[$i][$k] = 0.0;}$this->U[$k][$k] = 1.0;}}}// If required, generate V.if ($wantv) {for ($k = $this->n - 1; $k >= 0; --$k) {if (($k < $nrt) AND ($e[$k] != 0.0)) {for ($j = $k + 1; $j < $nu; ++$j) {$t = 0;for ($i = $k + 1; $i < $this->n; ++$i) {$t += $this->V[$i][$k]* $this->V[$i][$j];}$t = -$t / $this->V[$k+1][$k];for ($i = $k + 1; $i < $this->n; ++$i) {$this->V[$i][$j] += $t * $this->V[$i][$k];}}}for ($i = 0; $i < $this->n; ++$i) {$this->V[$i][$k] = 0.0;}$this->V[$k][$k] = 1.0;}}// Main iteration loop for the singular values.$pp = $p - 1;$iter = 0;$eps = pow(2.0, -52.0);while ($p > 0) {// Here is where a test for too many iterations would go.// This section of the program inspects for negligible// elements in the s and e arrays. On completion the// variables kase and k are set as follows:// kase = 1 if s(p) and e[k-1] are negligible and k<p// kase = 2 if s(k) is negligible and k<p// kase = 3 if e[k-1] is negligible, k<p, and// s(k), ..., s(p) are not negligible (qr step).// kase = 4 if e(p-1) is negligible (convergence).for ($k = $p - 2; $k >= -1; --$k) {if ($k == -1) {break;}if (abs($e[$k]) <= $eps * (abs($this->s[$k]) + abs($this->s[$k+1]))) {$e[$k] = 0.0;break;}}if ($k == $p - 2) {$kase = 4;} else {for ($ks = $p - 1; $ks >= $k; --$ks) {if ($ks == $k) {break;}$t = ($ks != $p ? abs($e[$ks]) : 0.) + ($ks != $k + 1 ? abs($e[$ks-1]) : 0.);if (abs($this->s[$ks]) <= $eps * $t) {$this->s[$ks] = 0.0;break;}}if ($ks == $k) {$kase = 3;} else if ($ks == $p-1) {$kase = 1;} else {$kase = 2;$k = $ks;}}++$k;// Perform the task indicated by kase.switch ($kase) {// Deflate negligible s(p).case 1:$f = $e[$p-2];$e[$p-2] = 0.0;for ($j = $p - 2; $j >= $k; --$j) {$t = hypo($this->s[$j],$f);$cs = $this->s[$j] / $t;$sn = $f / $t;$this->s[$j] = $t;if ($j != $k) {$f = -$sn * $e[$j-1];$e[$j-1] = $cs * $e[$j-1];}if ($wantv) {for ($i = 0; $i < $this->n; ++$i) {$t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$p-1];$this->V[$i][$p-1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$p-1];$this->V[$i][$j] = $t;}}}break;// Split at negligible s(k).case 2:$f = $e[$k-1];$e[$k-1] = 0.0;for ($j = $k; $j < $p; ++$j) {$t = hypo($this->s[$j], $f);$cs = $this->s[$j] / $t;$sn = $f / $t;$this->s[$j] = $t;$f = -$sn * $e[$j];$e[$j] = $cs * $e[$j];if ($wantu) {for ($i = 0; $i < $this->m; ++$i) {$t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$k-1];$this->U[$i][$k-1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$k-1];$this->U[$i][$j] = $t;}}}break;// Perform one qr step.case 3:// Calculate the shift.$scale = max(max(max(max(abs($this->s[$p-1]),abs($this->s[$p-2])),abs($e[$p-2])),abs($this->s[$k])), abs($e[$k]));$sp = $this->s[$p-1] / $scale;$spm1 = $this->s[$p-2] / $scale;$epm1 = $e[$p-2] / $scale;$sk = $this->s[$k] / $scale;$ek = $e[$k] / $scale;$b = (($spm1 + $sp) * ($spm1 - $sp) + $epm1 * $epm1) / 2.0;$c = ($sp * $epm1) * ($sp * $epm1);$shift = 0.0;if (($b != 0.0) || ($c != 0.0)) {$shift = sqrt($b * $b + $c);if ($b < 0.0) {$shift = -$shift;}$shift = $c / ($b + $shift);}$f = ($sk + $sp) * ($sk - $sp) + $shift;$g = $sk * $ek;// Chase zeros.for ($j = $k; $j < $p-1; ++$j) {$t = hypo($f,$g);$cs = $f/$t;$sn = $g/$t;if ($j != $k) {$e[$j-1] = $t;}$f = $cs * $this->s[$j] + $sn * $e[$j];$e[$j] = $cs * $e[$j] - $sn * $this->s[$j];$g = $sn * $this->s[$j+1];$this->s[$j+1] = $cs * $this->s[$j+1];if ($wantv) {for ($i = 0; $i < $this->n; ++$i) {$t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$j+1];$this->V[$i][$j+1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$j+1];$this->V[$i][$j] = $t;}}$t = hypo($f,$g);$cs = $f/$t;$sn = $g/$t;$this->s[$j] = $t;$f = $cs * $e[$j] + $sn * $this->s[$j+1];$this->s[$j+1] = -$sn * $e[$j] + $cs * $this->s[$j+1];$g = $sn * $e[$j+1];$e[$j+1] = $cs * $e[$j+1];if ($wantu && ($j < $this->m - 1)) {for ($i = 0; $i < $this->m; ++$i) {$t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$j+1];$this->U[$i][$j+1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$j+1];$this->U[$i][$j] = $t;}}}$e[$p-2] = $f;$iter = $iter + 1;break;// Convergence.case 4:// Make the singular values positive.if ($this->s[$k] <= 0.0) {$this->s[$k] = ($this->s[$k] < 0.0 ? -$this->s[$k] : 0.0);if ($wantv) {for ($i = 0; $i <= $pp; ++$i) {$this->V[$i][$k] = -$this->V[$i][$k];}}}// Order the singular values.while ($k < $pp) {if ($this->s[$k] >= $this->s[$k+1]) {break;}$t = $this->s[$k];$this->s[$k] = $this->s[$k+1];$this->s[$k+1] = $t;if ($wantv AND ($k < $this->n - 1)) {for ($i = 0; $i < $this->n; ++$i) {$t = $this->V[$i][$k+1];$this->V[$i][$k+1] = $this->V[$i][$k];$this->V[$i][$k] = $t;}}if ($wantu AND ($k < $this->m-1)) {for ($i = 0; $i < $this->m; ++$i) {$t = $this->U[$i][$k+1];$this->U[$i][$k+1] = $this->U[$i][$k];$this->U[$i][$k] = $t;}}++$k;}$iter = 0;--$p;break;} // end switch} // end while} // end constructor/*** Return the left singular vectors** @access public* @return U*/public function getU() {return new Matrix($this->U, $this->m, min($this->m + 1, $this->n));}/*** Return the right singular vectors** @access public* @return V*/public function getV() {return new Matrix($this->V);}/*** Return the one-dimensional array of singular values** @access public* @return diagonal of S.*/public function getSingularValues() {return $this->s;}/*** Return the diagonal matrix of singular values** @access public* @return S*/public function getS() {for ($i = 0; $i < $this->n; ++$i) {for ($j = 0; $j < $this->n; ++$j) {$S[$i][$j] = 0.0;}$S[$i][$i] = $this->s[$i];}return new Matrix($S);}/*** Two norm** @access public* @return max(S)*/public function norm2() {return $this->s[0];}/*** Two norm condition number** @access public* @return max(S)/min(S)*/public function cond() {return $this->s[0] / $this->s[min($this->m, $this->n) - 1];}/*** Effective numerical matrix rank** @access public* @return Number of nonnegligible singular values.*/public function rank() {$eps = pow(2.0, -52.0);$tol = max($this->m, $this->n) * $this->s[0] * $eps;$r = 0;for ($i = 0; $i < count($this->s); ++$i) {if ($this->s[$i] > $tol) {++$r;}}return $r;}} // class SingularValueDecomposition