ILIAS  release_5-3 Revision v5.3.23-19-g915713cf615
Public Member Functions | Private Attributes
SingularValueDecomposition Class Reference
+ Collaboration diagram for SingularValueDecomposition:

Public Member Functions

 __construct ($Arg)
 Construct the singular value decomposition. More...
 
 getU ()
 Return the left singular vectors. More...
 
 getV ()
 Return the right singular vectors. More...
 
 getSingularValues ()
 Return the one-dimensional array of singular values. More...
 
 getS ()
 Return the diagonal matrix of singular values. More...
 
 norm2 ()
 Two norm. More...
 
 cond ()
 Two norm condition number. More...
 
 rank ()
 Effective numerical matrix rank. More...
 

Private Attributes

 $U = array()
 
 $V = array()
 
 $s = array()
 
 $m
 
 $n
 

Detailed Description

Definition at line 20 of file SingularValueDecomposition.php.

Constructor & Destructor Documentation

◆ __construct()

SingularValueDecomposition::__construct (   $Arg)

Construct the singular value decomposition.

Derived from LINPACK code.

Parameters
$ARectangular matrix
Returns
Structure to access U, S and V.

Definition at line 61 of file SingularValueDecomposition.php.

References $eps, $i, $m, $n, $t, array, hypo(), n, and s.

61  {
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
n
if(! $in) print Initializing normalization quick check tables n
Definition: UtfNormalGenerate.php:37
$eps
$eps
Definition: metadata.php:61
hypo
hypo($a, $b)
Definition: Maths.php:14
array
Create styles array
The data for the language used.
Definition: 40duplicateStyle.php:19
s
Add data(end) s
Definition: 40duplicateStyle.php:39
$i
$i
Definition: disco.tpl.php:19
+ Here is the call graph for this function:

Member Function Documentation

◆ cond()

SingularValueDecomposition::cond ( )

Two norm condition number.

public

Returns
max(S)/min(S)

Definition at line 503 of file SingularValueDecomposition.php.

References n, and s.

503  {
504  return $this->s[0] / $this->s[min($this->m, $this->n) - 1];
505  }
n
if(! $in) print Initializing normalization quick check tables n
Definition: UtfNormalGenerate.php:37
s
Add data(end) s
Definition: 40duplicateStyle.php:39

◆ getS()

SingularValueDecomposition::getS ( )

Return the diagonal matrix of singular values.

public

Returns
S

Definition at line 475 of file SingularValueDecomposition.php.

References $i, $n, and s.

475  {
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  }
s
Add data(end) s
Definition: 40duplicateStyle.php:39
$i
$i
Definition: disco.tpl.php:19

◆ getSingularValues()

SingularValueDecomposition::getSingularValues ( )

Return the one-dimensional array of singular values.

public

Returns
diagonal of S.

Definition at line 464 of file SingularValueDecomposition.php.

References $s.

464  {
465  return $this->s;
466  }

◆ getU()

SingularValueDecomposition::getU ( )

Return the left singular vectors.

public

Returns
U

Definition at line 442 of file SingularValueDecomposition.php.

References n.

442  {
443  return new Matrix($this->U, $this->m, min($this->m + 1, $this->n));
444  }
n
if(! $in) print Initializing normalization quick check tables n
Definition: UtfNormalGenerate.php:37

◆ getV()

SingularValueDecomposition::getV ( )

Return the right singular vectors.

public

Returns
V

Definition at line 453 of file SingularValueDecomposition.php.

453  {
454  return new Matrix($this->V);
455  }

◆ norm2()

SingularValueDecomposition::norm2 ( )

Two norm.

public

Returns
max(S)

Definition at line 492 of file SingularValueDecomposition.php.

References s.

492  {
493  return $this->s[0];
494  }
s
Add data(end) s
Definition: 40duplicateStyle.php:39

◆ rank()

SingularValueDecomposition::rank ( )

Effective numerical matrix rank.

public

Returns
Number of nonnegligible singular values.

Definition at line 514 of file SingularValueDecomposition.php.

References $eps, $i, $r, n, and s.

514  {
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  }
n
if(! $in) print Initializing normalization quick check tables n
Definition: UtfNormalGenerate.php:37
$r
$r
Definition: example_031.php:79
$eps
$eps
Definition: metadata.php:61
s
Add data(end) s
Definition: 40duplicateStyle.php:39
$i
$i
Definition: disco.tpl.php:19

Field Documentation

◆ $m

SingularValueDecomposition::$m
private

Definition at line 44 of file SingularValueDecomposition.php.

Referenced by __construct().

◆ $n

SingularValueDecomposition::$n
private

Definition at line 50 of file SingularValueDecomposition.php.

Referenced by __construct(), and getS().

◆ $s

SingularValueDecomposition::$s = array()
private

Definition at line 38 of file SingularValueDecomposition.php.

Referenced by getSingularValues().

◆ $U

SingularValueDecomposition::$U = array()
private

Definition at line 26 of file SingularValueDecomposition.php.

◆ $V

SingularValueDecomposition::$V = array()
private

Definition at line 32 of file SingularValueDecomposition.php.

The documentation for this class was generated from the following file: