root / tmp / org.txm.analec.rcp / src / JamaPlus / SingularValueDecomposition.java0 @ 1726
History  View  Annotate  Download (15.8 kB)
1 
package JamaPlus; 

2 
import JamaPlus.util.*; 
3  
4 
/** Singular Value Decomposition. 
5 
<P> 
6 
For an mbyn matrix A with m >= n, the singular value decomposition is 
7 
an mbyn orthogonal matrix U, an nbyn diagonal matrix S, and 
8 
an nbyn orthogonal matrix V so that A = U*S*V'. 
9 
<P> 
10 
The singular values, sigma[k] = S[k][k], are ordered so that 
11 
sigma[0] >= sigma[1] >= ... >= sigma[n1]. 
12 
<P> 
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  
18 
public class SingularValueDecomposition implements java.io.Serializable { 
19  
20 
/*  
21 
Class variables 
22 
*  */ 
23  
24 
/** Arrays for internal storage of U and V. 
25 
@serial internal storage of U. 
26 
@serial internal storage of V. 
27 
*/ 
28 
private double[][] U, V; 
29  
30 
/** Array for internal storage of singular values. 
31 
@serial internal storage of singular values. 
32 
*/ 
33 
private double[] s; 
34  
35 
/** Row and column dimensions. 
36 
@serial row dimension. 
37 
@serial column dimension. 
38 
*/ 
39 
private int m, n; 
40  
41 
/*  
42 
Constructor 
43 
*  */ 
44  
45 
/** Construct the singular value decomposition 
46 
@param A Rectangular matrix 
47 
@return Structure to access U, S and V. 
48 
*/ 
49  
50 
public SingularValueDecomposition (Matrix Arg) { 
51  
52 
// Derived from LINPACK code. 
53 
// Initialize. 
54 
double[][] A = Arg.getArrayCopy(); 
55 
m = Arg.getRowDimension(); 
56 
n = Arg.getColumnDimension(); 
57  
58 
/* Apparently the failing cases are only a proper subset of (m<n), 
59 
so let's not throw error. Correct fix to come later? 
60 
if (m<n) { 
61 
throw new IllegalArgumentException("Jama SVD only works for m >= n"); } 
62 
*/ 
63 
int nu = Math.min(m,n); 
64 
s = new double [Math.min(m+1,n)]; 
65 
U = new double [m][nu]; 
66 
V = new double [n][n]; 
67 
double[] e = new double [n]; 
68 
double[] work = new double [m]; 
69 
boolean wantu = true; 
70 
boolean wantv = true; 
71  
72 
// Reduce A to bidiagonal form, storing the diagonal elements 
73 
// in s and the superdiagonal elements in e. 
74  
75 
int nct = Math.min(m1,n); 
76 
int nrt = Math.max(0,Math.min(n2,m)); 
77 
for (int k = 0; k < Math.max(nct,nrt); k++) { 
78 
if (k < nct) { 
79  
80 
// Compute the transformation for the kth column and 
81 
// place the kth diagonal in s[k]. 
82 
// Compute 2norm of kth column without under/overflow. 
83 
s[k] = 0; 
84 
for (int i = k; i < m; i++) { 
85 
s[k] = Math.hypot(s[k],A[i][k]); 
86 
} 
87 
if (s[k] != 0.0) { 
88 
if (A[k][k] < 0.0) { 
89 
s[k] = s[k]; 
90 
} 
91 
for (int i = k; i < m; i++) { 
92 
A[i][k] /= s[k]; 
93 
} 
94 
A[k][k] += 1.0; 
95 
} 
96 
s[k] = s[k]; 
97 
} 
98 
for (int j = k+1; j < n; j++) { 
99 
if ((k < nct) & (s[k] != 0.0)) { 
100  
101 
// Apply the transformation. 
102  
103 
double t = 0; 
104 
for (int i = k; i < m; i++) { 
105 
t += A[i][k]*A[i][j]; 
106 
} 
107 
t = t/A[k][k]; 
108 
for (int i = k; i < m; i++) { 
109 
A[i][j] += t*A[i][k]; 
110 
} 
111 
} 
112  
113 
// Place the kth row of A into e for the 
114 
// subsequent calculation of the row transformation. 
115  
116 
e[j] = A[k][j]; 
117 
} 
118 
if (wantu & (k < nct)) { 
119  
120 
// Place the transformation in U for subsequent back 
121 
// multiplication. 
122  
123 
for (int i = k; i < m; i++) { 
124 
U[i][k] = A[i][k]; 
125 
} 
126 
} 
127 
if (k < nrt) { 
128  
129 
// Compute the kth row transformation and place the 
130 
// kth superdiagonal in e[k]. 
131 
// Compute 2norm without under/overflow. 
132 
e[k] = 0; 
133 
for (int i = k+1; i < n; i++) { 
134 
e[k] = Math.hypot(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 (int i = k+1; i < n; i++) { 
141 
e[i] /= e[k]; 
142 
} 
143 
e[k+1] += 1.0; 
144 
} 
145 
e[k] = e[k]; 
146 
if ((k+1 < m) & (e[k] != 0.0)) { 
147  
148 
// Apply the transformation. 
149  
150 
for (int i = k+1; i < m; i++) { 
151 
work[i] = 0.0; 
152 
} 
153 
for (int j = k+1; j < n; j++) { 
154 
for (int i = k+1; i < m; i++) { 
155 
work[i] += e[j]*A[i][j]; 
156 
} 
157 
} 
158 
for (int j = k+1; j < n; j++) { 
159 
double t = e[j]/e[k+1]; 
160 
for (int i = k+1; i < m; i++) { 
161 
A[i][j] += t*work[i]; 
162 
} 
163 
} 
164 
} 
165 
if (wantv) { 
166  
167 
// Place the transformation in V for subsequent 
168 
// back multiplication. 
169  
170 
for (int i = k+1; i < n; i++) { 
171 
V[i][k] = e[i]; 
172 
} 
173 
} 
174 
} 
175 
} 
176  
177 
// Set up the final bidiagonal matrix or order p. 
178  
179 
int p = Math.min(n,m+1); 
180 
if (nct < n) { 
181 
s[nct] = A[nct][nct]; 
182 
} 
183 
if (m < p) { 
184 
s[p1] = 0.0; 
185 
} 
186 
if (nrt+1 < p) { 
187 
e[nrt] = A[nrt][p1]; 
188 
} 
189 
e[p1] = 0.0; 
190  
191 
// If required, generate U. 
192  
193 
if (wantu) { 
194 
for (int j = nct; j < nu; j++) { 
195 
for (int i = 0; i < m; i++) { 
196 
U[i][j] = 0.0; 
197 
} 
198 
U[j][j] = 1.0; 
199 
} 
200 
for (int k = nct1; k >= 0; k) { 
201 
if (s[k] != 0.0) { 
202 
for (int j = k+1; j < nu; j++) { 
203 
double t = 0; 
204 
for (int i = k; i < m; i++) { 
205 
t += U[i][k]*U[i][j]; 
206 
} 
207 
t = t/U[k][k]; 
208 
for (int i = k; i < m; i++) { 
209 
U[i][j] += t*U[i][k]; 
210 
} 
211 
} 
212 
for (int i = k; i < m; i++ ) { 
213 
U[i][k] = U[i][k]; 
214 
} 
215 
U[k][k] = 1.0 + U[k][k]; 
216 
for (int i = 0; i < k1; i++) { 
217 
U[i][k] = 0.0; 
218 
} 
219 
} else { 
220 
for (int i = 0; i < m; i++) { 
221 
U[i][k] = 0.0; 
222 
} 
223 
U[k][k] = 1.0; 
224 
} 
225 
} 
226 
} 
227  
228 
// If required, generate V. 
229  
230 
if (wantv) { 
231 
for (int k = n1; k >= 0; k) { 
232 
if ((k < nrt) & (e[k] != 0.0)) { 
233 
for (int j = k+1; j < nu; j++) { 
234 
double t = 0; 
235 
for (int i = k+1; i < n; i++) { 
236 
t += V[i][k]*V[i][j]; 
237 
} 
238 
t = t/V[k+1][k]; 
239 
for (int i = k+1; i < n; i++) { 
240 
V[i][j] += t*V[i][k]; 
241 
} 
242 
} 
243 
} 
244 
for (int i = 0; i < n; i++) { 
245 
V[i][k] = 0.0; 
246 
} 
247 
V[k][k] = 1.0; 
248 
} 
249 
} 
250  
251 
// Main iteration loop for the singular values. 
252  
253 
int pp = p1; 
254 
int iter = 0; 
255 
double eps = Math.pow(2.0,52.0); 
256 
double tiny = Math.pow(2.0,966.0); 
257 
while (p > 0) { 
258 
int k,kase; 
259  
260 
// Here is where a test for too many iterations would go. 
261  
262 
// This section of the program inspects for 
263 
// negligible elements in the s and e arrays. On 
264 
// completion the variables kase and k are set as follows. 
265  
266 
// kase = 1 if s(p) and e[k1] are negligible and k<p 
267 
// kase = 2 if s(k) is negligible and k<p 
268 
// kase = 3 if e[k1] is negligible, k<p, and 
269 
// s(k), ..., s(p) are not negligible (qr step). 
270 
// kase = 4 if e(p1) is negligible (convergence). 
271  
272 
for (k = p2; k >= 1; k) { 
273 
if (k == 1) { 
274 
break; 
275 
} 
276 
if (Math.abs(e[k]) <= 
277 
tiny + eps*(Math.abs(s[k]) + Math.abs(s[k+1]))) { 
278 
e[k] = 0.0; 
279 
break; 
280 
} 
281 
} 
282 
if (k == p2) { 
283 
kase = 4; 
284 
} else { 
285 
int ks; 
286 
for (ks = p1; ks >= k; ks) { 
287 
if (ks == k) { 
288 
break; 
289 
} 
290 
double t = (ks != p ? Math.abs(e[ks]) : 0.) + 
291 
(ks != k+1 ? Math.abs(e[ks1]) : 0.); 
292 
if (Math.abs(s[ks]) <= tiny + eps*t) { 
293 
s[ks] = 0.0; 
294 
break; 
295 
} 
296 
} 
297 
if (ks == k) { 
298 
kase = 3; 
299 
} else if (ks == p1) { 
300 
kase = 1; 
301 
} else { 
302 
kase = 2; 
303 
k = ks; 
304 
} 
305 
} 
306 
k++; 
307  
308 
// Perform the task indicated by kase. 
309  
310 
switch (kase) { 
311  
312 
// Deflate negligible s(p). 
313  
314 
case 1: { 
315 
double f = e[p2]; 
316 
e[p2] = 0.0; 
317 
for (int j = p2; j >= k; j) { 
318 
double t = Math.hypot(s[j],f); 
319 
double cs = s[j]/t; 
320 
double sn = f/t; 
321 
s[j] = t; 
322 
if (j != k) { 
323 
f = sn*e[j1]; 
324 
e[j1] = cs*e[j1]; 
325 
} 
326 
if (wantv) { 
327 
for (int i = 0; i < n; i++) { 
328 
t = cs*V[i][j] + sn*V[i][p1]; 
329 
V[i][p1] = sn*V[i][j] + cs*V[i][p1]; 
330 
V[i][j] = t; 
331 
} 
332 
} 
333 
} 
334 
} 
335 
break; 
336  
337 
// Split at negligible s(k). 
338  
339 
case 2: { 
340 
double f = e[k1]; 
341 
e[k1] = 0.0; 
342 
for (int j = k; j < p; j++) { 
343 
double t = Math.hypot(s[j],f); 
344 
double cs = s[j]/t; 
345 
double sn = f/t; 
346 
s[j] = t; 
347 
f = sn*e[j]; 
348 
e[j] = cs*e[j]; 
349 
if (wantu) { 
350 
for (int i = 0; i < m; i++) { 
351 
t = cs*U[i][j] + sn*U[i][k1]; 
352 
U[i][k1] = sn*U[i][j] + cs*U[i][k1]; 
353 
U[i][j] = t; 
354 
} 
355 
} 
356 
} 
357 
} 
358 
break; 
359  
360 
// Perform one qr step. 
361  
362 
case 3: { 
363  
364 
// Calculate the shift. 
365 

366 
double scale = Math.max(Math.max(Math.max(Math.max( 
367 
Math.abs(s[p1]),Math.abs(s[p2])),Math.abs(e[p2])), 
368 
Math.abs(s[k])),Math.abs(e[k])); 
369 
double sp = s[p1]/scale; 
370 
double spm1 = s[p2]/scale; 
371 
double epm1 = e[p2]/scale; 
372 
double sk = s[k]/scale; 
373 
double ek = e[k]/scale; 
374 
double b = ((spm1 + sp)*(spm1  sp) + epm1*epm1)/2.0; 
375 
double c = (sp*epm1)*(sp*epm1); 
376 
double shift = 0.0; 
377 
if ((b != 0.0)  (c != 0.0)) { 
378 
shift = Math.sqrt(b*b + c); 
379 
if (b < 0.0) { 
380 
shift = shift; 
381 
} 
382 
shift = c/(b + shift); 
383 
} 
384 
double f = (sk + sp)*(sk  sp) + shift; 
385 
double g = sk*ek; 
386 

387 
// Chase zeros. 
388 

389 
for (int j = k; j < p1; j++) { 
390 
double t = Math.hypot(f,g); 
391 
double cs = f/t; 
392 
double sn = g/t; 
393 
if (j != k) { 
394 
e[j1] = t; 
395 
} 
396 
f = cs*s[j] + sn*e[j]; 
397 
e[j] = cs*e[j]  sn*s[j]; 
398 
g = sn*s[j+1]; 
399 
s[j+1] = cs*s[j+1]; 
400 
if (wantv) { 
401 
for (int i = 0; i < n; i++) { 
402 
t = cs*V[i][j] + sn*V[i][j+1]; 
403 
V[i][j+1] = sn*V[i][j] + cs*V[i][j+1]; 
404 
V[i][j] = t; 
405 
} 
406 
} 
407 
t = Math.hypot(f,g); 
408 
cs = f/t; 
409 
sn = g/t; 
410 
s[j] = t; 
411 
f = cs*e[j] + sn*s[j+1]; 
412 
s[j+1] = sn*e[j] + cs*s[j+1]; 
413 
g = sn*e[j+1]; 
414 
e[j+1] = cs*e[j+1]; 
415 
if (wantu && (j < m1)) { 
416 
for (int i = 0; i < m; i++) { 
417 
t = cs*U[i][j] + sn*U[i][j+1]; 
418 
U[i][j+1] = sn*U[i][j] + cs*U[i][j+1]; 
419 
U[i][j] = t; 
420 
} 
421 
} 
422 
} 
423 
e[p2] = f; 
424 
iter = iter + 1; 
425 
} 
426 
break; 
427  
428 
// Convergence. 
429  
430 
case 4: { 
431  
432 
// Make the singular values positive. 
433 

434 
if (s[k] <= 0.0) { 
435 
s[k] = (s[k] < 0.0 ? s[k] : 0.0); 
436 
if (wantv) { 
437 
for (int i = 0; i <= pp; i++) { 
438 
V[i][k] = V[i][k]; 
439 
} 
440 
} 
441 
} 
442 

443 
// Order the singular values. 
444 

445 
while (k < pp) { 
446 
if (s[k] >= s[k+1]) { 
447 
break; 
448 
} 
449 
double t = s[k]; 
450 
s[k] = s[k+1]; 
451 
s[k+1] = t; 
452 
if (wantv && (k < n1)) { 
453 
for (int i = 0; i < n; i++) { 
454 
t = V[i][k+1]; V[i][k+1] = V[i][k]; V[i][k] = t; 
455 
} 
456 
} 
457 
if (wantu && (k < m1)) { 
458 
for (int i = 0; i < m; i++) { 
459 
t = U[i][k+1]; U[i][k+1] = U[i][k]; U[i][k] = t; 
460 
} 
461 
} 
462 
k++; 
463 
} 
464 
iter = 0; 
465 
p; 
466 
} 
467 
break; 
468 
} 
469 
} 
470 
} 
471  
472 
/*  
473 
Public Methods 
474 
*  */ 
475  
476 
/** Return the left singular vectors 
477 
@return U 
478 
*/ 
479  
480 
public Matrix getU () { 
481 
return new Matrix(U,m,Math.min(m+1,n)); 
482 
} 
483  
484 
/** Return the right singular vectors 
485 
@return V 
486 
*/ 
487  
488 
public Matrix getV () { 
489 
return new Matrix(V,n,n); 
490 
} 
491  
492 
/** Return the onedimensional array of singular values 
493 
@return diagonal of S. 
494 
*/ 
495  
496 
public double[] getSingularValues () { 
497 
return s; 
498 
} 
499  
500 
/** Return the diagonal matrix of singular values 
501 
@return S 
502 
*/ 
503  
504 
public Matrix getS () { 
505 
Matrix X = new Matrix(n,n); 
506 
double[][] S = X.getArray(); 
507 
for (int i = 0; i < n; i++) { 
508 
for (int j = 0; j < n; j++) { 
509 
S[i][j] = 0.0; 
510 
} 
511 
S[i][i] = this.s[i]; 
512 
} 
513 
return X; 
514 
} 
515  
516 
/** Two norm 
517 
@return max(S) 
518 
*/ 
519  
520 
public double norm2 () { 
521 
return s[0]; 
522 
} 
523  
524 
/** Two norm condition number 
525 
@return max(S)/min(S) 
526 
*/ 
527  
528 
public double cond () { 
529 
return s[0]/s[Math.min(m,n)1]; 
530 
} 
531  
532 
/** Effective numerical matrix rank 
533 
@return Number of nonnegligible singular values. 
534 
*/ 
535  
536 
public int rank () { 
537 
double eps = Math.pow(2.0,52.0); 
538 
double tol = Math.max(m,n)*s[0]*eps; 
539 
int r = 0; 
540 
for (int i = 0; i < s.length; i++) { 
541 
if (s[i] > tol) { 
542 
r++; 
543 
} 
544 
} 
545 
return r; 
546 
} 
547 
} 