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
527: