ILIAS  release_5-4 Revision v5.4.26-12-gabc799a52e6
SingularValueDecomposition.php
Go to the documentation of this file.
1<?php
2
4
22{
28 private $U = [];
29
35 private $V = [];
36
42 private $s = [];
43
49 private $m;
50
56 private $n;
57
65 public function __construct($Arg)
66 {
67 // Initialize.
68 $A = $Arg->getArray();
69 $this->m = $Arg->getRowDimension();
70 $this->n = $Arg->getColumnDimension();
71 $nu = min($this->m, $this->n);
72 $e = [];
73 $work = [];
74 $wantu = true;
75 $wantv = true;
76 $nct = min($this->m - 1, $this->n);
77 $nrt = max(0, min($this->n - 2, $this->m));
78
79 // Reduce A to bidiagonal form, storing the diagonal elements
80 // in s and the super-diagonal elements in e.
81 $kMax = max($nct, $nrt);
82 for ($k = 0; $k < $kMax; ++$k) {
83 if ($k < $nct) {
84 // Compute the transformation for the k-th column and
85 // place the k-th diagonal in s[$k].
86 // Compute 2-norm of k-th column without under/overflow.
87 $this->s[$k] = 0;
88 for ($i = $k; $i < $this->m; ++$i) {
89 $this->s[$k] = hypo($this->s[$k], $A[$i][$k]);
90 }
91 if ($this->s[$k] != 0.0) {
92 if ($A[$k][$k] < 0.0) {
93 $this->s[$k] = -$this->s[$k];
94 }
95 for ($i = $k; $i < $this->m; ++$i) {
96 $A[$i][$k] /= $this->s[$k];
97 }
98 $A[$k][$k] += 1.0;
99 }
100 $this->s[$k] = -$this->s[$k];
101 }
102
103 for ($j = $k + 1; $j < $this->n; ++$j) {
104 if (($k < $nct) & ($this->s[$k] != 0.0)) {
105 // Apply the transformation.
106 $t = 0;
107 for ($i = $k; $i < $this->m; ++$i) {
108 $t += $A[$i][$k] * $A[$i][$j];
109 }
110 $t = -$t / $A[$k][$k];
111 for ($i = $k; $i < $this->m; ++$i) {
112 $A[$i][$j] += $t * $A[$i][$k];
113 }
114 // Place the k-th row of A into e for the
115 // subsequent calculation of the row transformation.
116 $e[$j] = $A[$k][$j];
117 }
118 }
119
120 if ($wantu && ($k < $nct)) {
121 // Place the transformation in U for subsequent back
122 // multiplication.
123 for ($i = $k; $i < $this->m; ++$i) {
124 $this->U[$i][$k] = $A[$i][$k];
125 }
126 }
127
128 if ($k < $nrt) {
129 // Compute the k-th row transformation and place the
130 // k-th super-diagonal in e[$k].
131 // Compute 2-norm without under/overflow.
132 $e[$k] = 0;
133 for ($i = $k + 1; $i < $this->n; ++$i) {
134 $e[$k] = hypo($e[$k], $e[$i]);
135 }
136 if ($e[$k] != 0.0) {
137 if ($e[$k + 1] < 0.0) {
138 $e[$k] = -$e[$k];
139 }
140 for ($i = $k + 1; $i < $this->n; ++$i) {
141 $e[$i] /= $e[$k];
142 }
143 $e[$k + 1] += 1.0;
144 }
145 $e[$k] = -$e[$k];
146 if (($k + 1 < $this->m) && ($e[$k] != 0.0)) {
147 // Apply the transformation.
148 for ($i = $k + 1; $i < $this->m; ++$i) {
149 $work[$i] = 0.0;
150 }
151 for ($j = $k + 1; $j < $this->n; ++$j) {
152 for ($i = $k + 1; $i < $this->m; ++$i) {
153 $work[$i] += $e[$j] * $A[$i][$j];
154 }
155 }
156 for ($j = $k + 1; $j < $this->n; ++$j) {
157 $t = -$e[$j] / $e[$k + 1];
158 for ($i = $k + 1; $i < $this->m; ++$i) {
159 $A[$i][$j] += $t * $work[$i];
160 }
161 }
162 }
163 if ($wantv) {
164 // Place the transformation in V for subsequent
165 // back multiplication.
166 for ($i = $k + 1; $i < $this->n; ++$i) {
167 $this->V[$i][$k] = $e[$i];
168 }
169 }
170 }
171 }
172
173 // Set up the final bidiagonal matrix or order p.
174 $p = min($this->n, $this->m + 1);
175 if ($nct < $this->n) {
176 $this->s[$nct] = $A[$nct][$nct];
177 }
178 if ($this->m < $p) {
179 $this->s[$p - 1] = 0.0;
180 }
181 if ($nrt + 1 < $p) {
182 $e[$nrt] = $A[$nrt][$p - 1];
183 }
184 $e[$p - 1] = 0.0;
185 // If required, generate U.
186 if ($wantu) {
187 for ($j = $nct; $j < $nu; ++$j) {
188 for ($i = 0; $i < $this->m; ++$i) {
189 $this->U[$i][$j] = 0.0;
190 }
191 $this->U[$j][$j] = 1.0;
192 }
193 for ($k = $nct - 1; $k >= 0; --$k) {
194 if ($this->s[$k] != 0.0) {
195 for ($j = $k + 1; $j < $nu; ++$j) {
196 $t = 0;
197 for ($i = $k; $i < $this->m; ++$i) {
198 $t += $this->U[$i][$k] * $this->U[$i][$j];
199 }
200 $t = -$t / $this->U[$k][$k];
201 for ($i = $k; $i < $this->m; ++$i) {
202 $this->U[$i][$j] += $t * $this->U[$i][$k];
203 }
204 }
205 for ($i = $k; $i < $this->m; ++$i) {
206 $this->U[$i][$k] = -$this->U[$i][$k];
207 }
208 $this->U[$k][$k] = 1.0 + $this->U[$k][$k];
209 for ($i = 0; $i < $k - 1; ++$i) {
210 $this->U[$i][$k] = 0.0;
211 }
212 } else {
213 for ($i = 0; $i < $this->m; ++$i) {
214 $this->U[$i][$k] = 0.0;
215 }
216 $this->U[$k][$k] = 1.0;
217 }
218 }
219 }
220
221 // If required, generate V.
222 if ($wantv) {
223 for ($k = $this->n - 1; $k >= 0; --$k) {
224 if (($k < $nrt) && ($e[$k] != 0.0)) {
225 for ($j = $k + 1; $j < $nu; ++$j) {
226 $t = 0;
227 for ($i = $k + 1; $i < $this->n; ++$i) {
228 $t += $this->V[$i][$k] * $this->V[$i][$j];
229 }
230 $t = -$t / $this->V[$k + 1][$k];
231 for ($i = $k + 1; $i < $this->n; ++$i) {
232 $this->V[$i][$j] += $t * $this->V[$i][$k];
233 }
234 }
235 }
236 for ($i = 0; $i < $this->n; ++$i) {
237 $this->V[$i][$k] = 0.0;
238 }
239 $this->V[$k][$k] = 1.0;
240 }
241 }
242
243 // Main iteration loop for the singular values.
244 $pp = $p - 1;
245 $iter = 0;
246 $eps = 2.0 ** (-52.0);
247
248 while ($p > 0) {
249 // Here is where a test for too many iterations would go.
250 // This section of the program inspects for negligible
251 // elements in the s and e arrays. On completion the
252 // variables kase and k are set as follows:
253 // kase = 1 if s(p) and e[k-1] are negligible and k<p
254 // kase = 2 if s(k) is negligible and k<p
255 // kase = 3 if e[k-1] is negligible, k<p, and
256 // s(k), ..., s(p) are not negligible (qr step).
257 // kase = 4 if e(p-1) is negligible (convergence).
258 for ($k = $p - 2; $k >= -1; --$k) {
259 if ($k == -1) {
260 break;
261 }
262 if (abs($e[$k]) <= $eps * (abs($this->s[$k]) + abs($this->s[$k + 1]))) {
263 $e[$k] = 0.0;
264
265 break;
266 }
267 }
268 if ($k == $p - 2) {
269 $kase = 4;
270 } else {
271 for ($ks = $p - 1; $ks >= $k; --$ks) {
272 if ($ks == $k) {
273 break;
274 }
275 $t = ($ks != $p ? abs($e[$ks]) : 0.) + ($ks != $k + 1 ? abs($e[$ks - 1]) : 0.);
276 if (abs($this->s[$ks]) <= $eps * $t) {
277 $this->s[$ks] = 0.0;
278
279 break;
280 }
281 }
282 if ($ks == $k) {
283 $kase = 3;
284 } elseif ($ks == $p - 1) {
285 $kase = 1;
286 } else {
287 $kase = 2;
288 $k = $ks;
289 }
290 }
291 ++$k;
292
293 // Perform the task indicated by kase.
294 switch ($kase) {
295 // Deflate negligible s(p).
296 case 1:
297 $f = $e[$p - 2];
298 $e[$p - 2] = 0.0;
299 for ($j = $p - 2; $j >= $k; --$j) {
300 $t = hypo($this->s[$j], $f);
301 $cs = $this->s[$j] / $t;
302 $sn = $f / $t;
303 $this->s[$j] = $t;
304 if ($j != $k) {
305 $f = -$sn * $e[$j - 1];
306 $e[$j - 1] = $cs * $e[$j - 1];
307 }
308 if ($wantv) {
309 for ($i = 0; $i < $this->n; ++$i) {
310 $t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$p - 1];
311 $this->V[$i][$p - 1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$p - 1];
312 $this->V[$i][$j] = $t;
313 }
314 }
315 }
316
317 break;
318 // Split at negligible s(k).
319 case 2:
320 $f = $e[$k - 1];
321 $e[$k - 1] = 0.0;
322 for ($j = $k; $j < $p; ++$j) {
323 $t = hypo($this->s[$j], $f);
324 $cs = $this->s[$j] / $t;
325 $sn = $f / $t;
326 $this->s[$j] = $t;
327 $f = -$sn * $e[$j];
328 $e[$j] = $cs * $e[$j];
329 if ($wantu) {
330 for ($i = 0; $i < $this->m; ++$i) {
331 $t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$k - 1];
332 $this->U[$i][$k - 1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$k - 1];
333 $this->U[$i][$j] = $t;
334 }
335 }
336 }
337
338 break;
339 // Perform one qr step.
340 case 3:
341 // Calculate the shift.
342 $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]));
343 $sp = $this->s[$p - 1] / $scale;
344 $spm1 = $this->s[$p - 2] / $scale;
345 $epm1 = $e[$p - 2] / $scale;
346 $sk = $this->s[$k] / $scale;
347 $ek = $e[$k] / $scale;
348 $b = (($spm1 + $sp) * ($spm1 - $sp) + $epm1 * $epm1) / 2.0;
349 $c = ($sp * $epm1) * ($sp * $epm1);
350 $shift = 0.0;
351 if (($b != 0.0) || ($c != 0.0)) {
352 $shift = sqrt($b * $b + $c);
353 if ($b < 0.0) {
354 $shift = -$shift;
355 }
356 $shift = $c / ($b + $shift);
357 }
358 $f = ($sk + $sp) * ($sk - $sp) + $shift;
359 $g = $sk * $ek;
360 // Chase zeros.
361 for ($j = $k; $j < $p - 1; ++$j) {
362 $t = hypo($f, $g);
363 $cs = $f / $t;
364 $sn = $g / $t;
365 if ($j != $k) {
366 $e[$j - 1] = $t;
367 }
368 $f = $cs * $this->s[$j] + $sn * $e[$j];
369 $e[$j] = $cs * $e[$j] - $sn * $this->s[$j];
370 $g = $sn * $this->s[$j + 1];
371 $this->s[$j + 1] = $cs * $this->s[$j + 1];
372 if ($wantv) {
373 for ($i = 0; $i < $this->n; ++$i) {
374 $t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$j + 1];
375 $this->V[$i][$j + 1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$j + 1];
376 $this->V[$i][$j] = $t;
377 }
378 }
379 $t = hypo($f, $g);
380 $cs = $f / $t;
381 $sn = $g / $t;
382 $this->s[$j] = $t;
383 $f = $cs * $e[$j] + $sn * $this->s[$j + 1];
384 $this->s[$j + 1] = -$sn * $e[$j] + $cs * $this->s[$j + 1];
385 $g = $sn * $e[$j + 1];
386 $e[$j + 1] = $cs * $e[$j + 1];
387 if ($wantu && ($j < $this->m - 1)) {
388 for ($i = 0; $i < $this->m; ++$i) {
389 $t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$j + 1];
390 $this->U[$i][$j + 1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$j + 1];
391 $this->U[$i][$j] = $t;
392 }
393 }
394 }
395 $e[$p - 2] = $f;
396 $iter = $iter + 1;
397
398 break;
399 // Convergence.
400 case 4:
401 // Make the singular values positive.
402 if ($this->s[$k] <= 0.0) {
403 $this->s[$k] = ($this->s[$k] < 0.0 ? -$this->s[$k] : 0.0);
404 if ($wantv) {
405 for ($i = 0; $i <= $pp; ++$i) {
406 $this->V[$i][$k] = -$this->V[$i][$k];
407 }
408 }
409 }
410 // Order the singular values.
411 while ($k < $pp) {
412 if ($this->s[$k] >= $this->s[$k + 1]) {
413 break;
414 }
415 $t = $this->s[$k];
416 $this->s[$k] = $this->s[$k + 1];
417 $this->s[$k + 1] = $t;
418 if ($wantv && ($k < $this->n - 1)) {
419 for ($i = 0; $i < $this->n; ++$i) {
420 $t = $this->V[$i][$k + 1];
421 $this->V[$i][$k + 1] = $this->V[$i][$k];
422 $this->V[$i][$k] = $t;
423 }
424 }
425 if ($wantu && ($k < $this->m - 1)) {
426 for ($i = 0; $i < $this->m; ++$i) {
427 $t = $this->U[$i][$k + 1];
428 $this->U[$i][$k + 1] = $this->U[$i][$k];
429 $this->U[$i][$k] = $t;
430 }
431 }
432 ++$k;
433 }
434 $iter = 0;
435 --$p;
436
437 break;
438 } // end switch
439 } // end while
440 }
441
447 public function getU()
448 {
449 return new Matrix($this->U, $this->m, min($this->m + 1, $this->n));
450 }
451
457 public function getV()
458 {
459 return new Matrix($this->V);
460 }
461
467 public function getSingularValues()
468 {
469 return $this->s;
470 }
471
477 public function getS()
478 {
479 $S = [];
480 for ($i = 0; $i < $this->n; ++$i) {
481 for ($j = 0; $j < $this->n; ++$j) {
482 $S[$i][$j] = 0.0;
483 }
484 $S[$i][$i] = $this->s[$i];
485 }
486
487 return new Matrix($S);
488 }
489
495 public function norm2()
496 {
497 return $this->s[0];
498 }
499
505 public function cond()
506 {
507 return $this->s[0] / $this->s[min($this->m, $this->n) - 1];
508 }
509
515 public function rank()
516 {
517 $eps = 2.0 ** (-52.0);
518 $tol = max($this->m, $this->n) * $this->s[0] * $eps;
519 $r = 0;
520 $iMax = count($this->s);
521 for ($i = 0; $i < $iMax; ++$i) {
522 if ($this->s[$i] > $tol) {
523 ++$r;
524 }
525 }
526
527 return $r;
528 }
529}
hypo($a, $b)
Pythagorean Theorem:.
Definition: Maths.php:18
An exception for terminatinating execution or to throw for unit testing.
For an m-by-n matrix A with m >= n, the singular value decomposition is an m-by-n orthogonal matrix U...
__construct($Arg)
Construct the singular value decomposition.
getSingularValues()
Return the one-dimensional array of singular values.
$i
Definition: disco.tpl.php:19
$r
Definition: example_031.php:79
$eps
Definition: metadata.php:61
Class for the creating "special" Matrices.
Definition: Builder.php:11