ILIAS  release_5-3 Revision v5.3.23-19-g915713cf615
SingularValueDecomposition.php
Go to the documentation of this file.
1<?php
20class SingularValueDecomposition {
21
26 private $U = array();
27
32 private $V = array();
33
38 private $s = array();
39
44 private $m;
45
50 private $n;
51
52
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
442 public function getU() {
443 return new Matrix($this->U, $this->m, min($this->m + 1, $this->n));
444 }
445
446
453 public function getV() {
454 return new Matrix($this->V);
455 }
456
457
464 public function getSingularValues() {
465 return $this->s;
466 }
467
468
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
492 public function norm2() {
493 return $this->s[0];
494 }
495
496
503 public function cond() {
504 return $this->s[0] / $this->s[min($this->m, $this->n) - 1];
505 }
506
507
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
hypo
hypo($a, $b)
Definition: Maths.php:14
php
An exception for terminatinating execution or to throw for unit testing.
SingularValueDecomposition\getU
getU()
Return the left singular vectors.
Definition: SingularValueDecomposition.php:442
SingularValueDecomposition\cond
cond()
Two norm condition number.
Definition: SingularValueDecomposition.php:503
SingularValueDecomposition\rank
rank()
Effective numerical matrix rank.
Definition: SingularValueDecomposition.php:514
SingularValueDecomposition\__construct
__construct($Arg)
Construct the singular value decomposition.
Definition: SingularValueDecomposition.php:61
SingularValueDecomposition\getSingularValues
getSingularValues()
Return the one-dimensional array of singular values.
Definition: SingularValueDecomposition.php:464
SingularValueDecomposition\getS
getS()
Return the diagonal matrix of singular values.
Definition: SingularValueDecomposition.php:475
SingularValueDecomposition\getV
getV()
Return the right singular vectors.
Definition: SingularValueDecomposition.php:453
$i
$i
Definition: disco.tpl.php:19
$r
$r
Definition: example_031.php:79
$eps
$eps
Definition: metadata.php:61