ILIAS  Release_4_0_x_branch Revision 61816
 All Data Structures Namespaces Files Functions Variables Groups Pages
SingularValueDecomposition.php
Go to the documentation of this file.
1 <?php
20 class 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