123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528 |
- <?php
- class SingularValueDecomposition
- {
-
- private $U = array();
-
- private $V = array();
-
- private $s = array();
-
- private $m;
-
- private $n;
-
- public function __construct($Arg)
- {
-
- $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));
-
-
- for ($k = 0; $k < max($nct, $nrt); ++$k) {
- if ($k < $nct) {
-
-
-
- $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)) {
-
- $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];
- }
-
-
- $e[$j] = $A[$k][$j];
- }
- }
- if ($wantu and ($k < $nct)) {
-
-
- for ($i = $k; $i < $this->m; ++$i) {
- $this->U[$i][$k] = $A[$i][$k];
- }
- }
- if ($k < $nrt) {
-
-
-
- $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)) {
-
- 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) {
-
-
- for ($i = $k + 1; $i < $this->n; ++$i) {
- $this->V[$i][$k] = $e[$i];
- }
- }
- }
- }
-
- $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 ($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 ($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;
- }
- }
-
- $pp = $p - 1;
- $iter = 0;
- $eps = pow(2.0, -52.0);
- while ($p > 0) {
-
-
-
-
-
-
-
-
-
- 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;
- } elseif ($ks == $p-1) {
- $kase = 1;
- } else {
- $kase = 2;
- $k = $ks;
- }
- }
- ++$k;
-
- switch ($kase) {
-
- 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;
-
- 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;
-
- case 3:
-
- $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;
-
- 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;
-
- case 4:
-
- 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];
- }
- }
- }
-
- 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;
- }
- }
- }
-
- public function getU()
- {
- return new Matrix($this->U, $this->m, min($this->m + 1, $this->n));
- }
-
- public function getV()
- {
- return new Matrix($this->V);
- }
-
- public function getSingularValues()
- {
- return $this->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);
- }
-
- public function norm2()
- {
- return $this->s[0];
- }
-
- public function cond()
- {
- return $this->s[0] / $this->s[min($this->m, $this->n) - 1];
- }
-
- 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;
- }
- }
|