|
1
|
|
|
<?php |
|
2
|
|
|
|
|
3
|
|
|
namespace PhpOffice\PhpSpreadsheet\Shared\JAMA; |
|
4
|
|
|
|
|
5
|
|
|
/** |
|
6
|
|
|
* For an m-by-n matrix A with m >= n, the singular value decomposition is |
|
7
|
|
|
* an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and |
|
8
|
|
|
* an n-by-n orthogonal matrix V so that A = U*S*V'. |
|
9
|
|
|
* |
|
10
|
|
|
* The singular values, sigma[$k] = S[$k][$k], are ordered so that |
|
11
|
|
|
* sigma[0] >= sigma[1] >= ... >= sigma[n-1]. |
|
12
|
|
|
* |
|
13
|
|
|
* The singular value decompostion always exists, so the constructor will |
|
14
|
|
|
* never fail. The matrix condition number and the effective numerical |
|
15
|
|
|
* rank can be computed from this decomposition. |
|
16
|
|
|
* |
|
17
|
|
|
* @author Paul Meagher |
|
18
|
|
|
* |
|
19
|
|
|
* @version 1.1 |
|
20
|
|
|
*/ |
|
21
|
|
|
class SingularValueDecomposition |
|
22
|
|
|
{ |
|
23
|
|
|
/** |
|
24
|
|
|
* Internal storage of U. |
|
25
|
|
|
* |
|
26
|
|
|
* @var array |
|
27
|
|
|
*/ |
|
28
|
|
|
private $U = []; |
|
29
|
|
|
|
|
30
|
|
|
/** |
|
31
|
|
|
* Internal storage of V. |
|
32
|
|
|
* |
|
33
|
|
|
* @var array |
|
34
|
|
|
*/ |
|
35
|
|
|
private $V = []; |
|
36
|
|
|
|
|
37
|
|
|
/** |
|
38
|
|
|
* Internal storage of singular values. |
|
39
|
|
|
* |
|
40
|
|
|
* @var array |
|
41
|
|
|
*/ |
|
42
|
|
|
private $s = []; |
|
43
|
|
|
|
|
44
|
|
|
/** |
|
45
|
|
|
* Row dimension. |
|
46
|
|
|
* |
|
47
|
|
|
* @var int |
|
48
|
|
|
*/ |
|
49
|
|
|
private $m; |
|
50
|
|
|
|
|
51
|
|
|
/** |
|
52
|
|
|
* Column dimension. |
|
53
|
|
|
* |
|
54
|
|
|
* @var int |
|
55
|
|
|
*/ |
|
56
|
|
|
private $n; |
|
57
|
|
|
|
|
58
|
|
|
/** |
|
59
|
|
|
* Construct the singular value decomposition. |
|
60
|
|
|
* |
|
61
|
|
|
* Derived from LINPACK code. |
|
62
|
|
|
* |
|
63
|
|
|
* @param mixed $Arg Rectangular matrix |
|
64
|
|
|
*/ |
|
65
|
|
|
public function __construct($Arg) |
|
66
|
|
|
{ |
|
67
|
|
|
// Initialize. |
|
68
|
|
|
$A = $Arg->getArrayCopy(); |
|
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 and ($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) and ($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) and ($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 = pow(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 and ($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 and ($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
|
|
|
|
|
442
|
|
|
/** |
|
443
|
|
|
* Return the left singular vectors. |
|
444
|
|
|
* |
|
445
|
|
|
* @return Matrix U |
|
446
|
|
|
*/ |
|
447
|
|
|
public function getU() |
|
448
|
|
|
{ |
|
449
|
|
|
return new Matrix($this->U, $this->m, min($this->m + 1, $this->n)); |
|
450
|
|
|
} |
|
451
|
|
|
|
|
452
|
|
|
/** |
|
453
|
|
|
* Return the right singular vectors. |
|
454
|
|
|
* |
|
455
|
|
|
* @return Matrix V |
|
456
|
|
|
*/ |
|
457
|
|
|
public function getV() |
|
458
|
|
|
{ |
|
459
|
|
|
return new Matrix($this->V); |
|
460
|
|
|
} |
|
461
|
|
|
|
|
462
|
|
|
/** |
|
463
|
|
|
* Return the one-dimensional array of singular values. |
|
464
|
|
|
* |
|
465
|
|
|
* @return array diagonal of S |
|
466
|
|
|
*/ |
|
467
|
|
|
public function getSingularValues() |
|
468
|
|
|
{ |
|
469
|
|
|
return $this->s; |
|
470
|
|
|
} |
|
471
|
|
|
|
|
472
|
|
|
/** |
|
473
|
|
|
* Return the diagonal matrix of singular values. |
|
474
|
|
|
* |
|
475
|
|
|
* @return Matrix S |
|
476
|
|
|
*/ |
|
477
|
|
|
public function getS() |
|
478
|
|
|
{ |
|
479
|
|
|
for ($i = 0; $i < $this->n; ++$i) { |
|
480
|
|
|
for ($j = 0; $j < $this->n; ++$j) { |
|
481
|
|
|
$S[$i][$j] = 0.0; |
|
482
|
|
|
} |
|
483
|
|
|
$S[$i][$i] = $this->s[$i]; |
|
484
|
|
|
} |
|
485
|
|
|
|
|
486
|
|
|
return new Matrix($S); |
|
|
|
|
|
|
487
|
|
|
} |
|
488
|
|
|
|
|
489
|
|
|
/** |
|
490
|
|
|
* Two norm. |
|
491
|
|
|
* |
|
492
|
|
|
* @return float max(S) |
|
493
|
|
|
*/ |
|
494
|
|
|
public function norm2() |
|
495
|
|
|
{ |
|
496
|
|
|
return $this->s[0]; |
|
497
|
|
|
} |
|
498
|
|
|
|
|
499
|
|
|
/** |
|
500
|
|
|
* Two norm condition number. |
|
501
|
|
|
* |
|
502
|
|
|
* @return float max(S)/min(S) |
|
503
|
|
|
*/ |
|
504
|
|
|
public function cond() |
|
505
|
|
|
{ |
|
506
|
|
|
return $this->s[0] / $this->s[min($this->m, $this->n) - 1]; |
|
507
|
|
|
} |
|
508
|
|
|
|
|
509
|
|
|
/** |
|
510
|
|
|
* Effective numerical matrix rank. |
|
511
|
|
|
* |
|
512
|
|
|
* @return int Number of nonnegligible singular values |
|
513
|
|
|
*/ |
|
514
|
|
|
public function rank() |
|
515
|
|
|
{ |
|
516
|
|
|
$eps = pow(2.0, -52.0); |
|
517
|
|
|
$tol = max($this->m, $this->n) * $this->s[0] * $eps; |
|
518
|
|
|
$r = 0; |
|
519
|
|
|
$iMax = count($this->s); |
|
520
|
|
|
for ($i = 0; $i < $iMax; ++$i) { |
|
521
|
|
|
if ($this->s[$i] > $tol) { |
|
522
|
|
|
++$r; |
|
523
|
|
|
} |
|
524
|
|
|
} |
|
525
|
|
|
|
|
526
|
|
|
return $r; |
|
527
|
|
|
} |
|
528
|
|
|
} |
|
529
|
|
|
|
PHPOffice /
PhpSpreadsheet
Loading history...