Subversion Repositories eFlore/Applications.cel

Rev

Rev 2388 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2388 jpm 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