ILIAS  release_5-2 Revision v5.2.25-18-g3f80b82851
SingularValueDecomposition.php
Go to the documentation of this file.
1 <?php
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
__construct($Arg)
Construct the singular value decomposition.
getU()
Return the left singular vectors.
getV()
Return the right singular vectors.
if(! $in) print Initializing normalization quick check tables n
$r
Definition: example_031.php:79
hypo($a, $b)
Definition: Maths.php:14
Create styles array
The data for the language used.
getS()
Return the diagonal matrix of singular values.
cond()
Two norm condition number.
rank()
Effective numerical matrix rank.
getSingularValues()
Return the one-dimensional array of singular values.