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.

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
hypo
hypo($a, $b)
Definition: Maths.php:14
$i
$i
Definition: disco.tpl.php:19
$eps
$eps
Definition: metadata.php:61

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

+ Here is the call graph for this function:

Member Function Documentation

◆ cond()

SingularValueDecomposition::cond ( )

Two norm condition number.

@access public

Returns
max(S)/min(S)

Definition at line 503 of file SingularValueDecomposition.php.

503 {
504 return $this->s[0] / $this->s[min($this->m, $this->n) - 1];
505 }

◆ getS()

SingularValueDecomposition::getS ( )

Return the diagonal matrix of singular values.

@access public

Returns
S

Definition at line 475 of file SingularValueDecomposition.php.

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 }

References $i, and $n.

◆ getSingularValues()

SingularValueDecomposition::getSingularValues ( )

Return the one-dimensional array of singular values.

@access public

Returns
diagonal of S.

Definition at line 464 of file SingularValueDecomposition.php.

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

References $s.

◆ getU()

SingularValueDecomposition::getU ( )

Return the left singular vectors.

@access public

Returns
U

Definition at line 442 of file SingularValueDecomposition.php.

442 {
443 return new Matrix($this->U, $this->m, min($this->m + 1, $this->n));
444 }

◆ getV()

SingularValueDecomposition::getV ( )

Return the right singular vectors.

@access public

Returns
V

Definition at line 453 of file SingularValueDecomposition.php.

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

◆ norm2()

SingularValueDecomposition::norm2 ( )

Two norm.

@access public

Returns
max(S)

Definition at line 492 of file SingularValueDecomposition.php.

492 {
493 return $this->s[0];
494 }

◆ rank()

SingularValueDecomposition::rank ( )

Effective numerical matrix rank.

@access public

Returns
Number of nonnegligible singular values.

Definition at line 514 of file SingularValueDecomposition.php.

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 }
$r
$r
Definition: example_031.php:79

References $eps, $i, and $r.

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: