Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

SingularValueDecomposition.php @ 1

Last change on this file since 1 was 1, checked in by dungnv, 11 years ago

File size: 12.7 KB

Line
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 = USV'.
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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: sourcecode/application/libraries/PHPExcel/Shared/JAMA/SingularValueDecomposition.php @ 1

Download in other formats: