diff options
Diffstat (limited to 'admin/survey/excel/PHPExcel/Shared/JAMA/SingularValueDecomposition.php')
| -rw-r--r-- | admin/survey/excel/PHPExcel/Shared/JAMA/SingularValueDecomposition.php | 526 |
1 files changed, 526 insertions, 0 deletions
diff --git a/admin/survey/excel/PHPExcel/Shared/JAMA/SingularValueDecomposition.php b/admin/survey/excel/PHPExcel/Shared/JAMA/SingularValueDecomposition.php new file mode 100644 index 0000000..ecfb0cd --- /dev/null +++ b/admin/survey/excel/PHPExcel/Shared/JAMA/SingularValueDecomposition.php @@ -0,0 +1,526 @@ +<?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
|