001/* 002 * BioJava development code 003 * 004 * This code may be freely distributed and modified under the 005 * terms of the GNU Lesser General Public Licence. This should 006 * be distributed with the code. If you do not have a copy, 007 * see: 008 * 009 * http://www.gnu.org/copyleft/lesser.html 010 * 011 * Copyright for this code is held jointly by the individual 012 * authors. These should be listed in @author doc comments. 013 * 014 * For more information on the BioJava project and its aims, 015 * or to join the biojava-l mailing list, visit the home page 016 * at: 017 * 018 * http://www.biojava.org/ 019 * 020 */ 021package org.biojava.nbio.structure.jama; 022 023 /** Singular Value Decomposition. 024 <P> 025 For an m-by-n matrix A with m >= n, the singular value decomposition is 026 an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and 027 an n-by-n orthogonal matrix V so that A = U*S*V'. 028 <P> 029 The singular values, sigma[k] = S[k][k], are ordered so that 030 sigma[0] >= sigma[1] >= ... >= sigma[n-1]. 031 <P> 032 The singular value decompostion always exists, so the constructor will 033 never fail. The matrix condition number and the effective numerical 034 rank can be computed from this decomposition. 035 */ 036 037public class SingularValueDecomposition implements java.io.Serializable { 038 039 static final long serialVersionUID = 640239472093534756l; 040 041/* ------------------------ 042 Class variables 043 * ------------------------ */ 044 045 /** Arrays for internal storage of U and V. 046 @serial internal storage of U. 047 @serial internal storage of V. 048 */ 049 private double[][] U, V; 050 051 /** Array for internal storage of singular values. 052 @serial internal storage of singular values. 053 */ 054 private double[] s; 055 056 /** Row and column dimensions. 057 @serial row dimension. 058 @serial column dimension. 059 */ 060 private int m, n; 061 062/* ------------------------ 063 Constructor 064 * ------------------------ */ 065 066 /** Construct the singular value decomposition. Provides a data structure to access U, S and V. 067 @param Arg Rectangular matrix 068 */ 069 070 public SingularValueDecomposition (Matrix Arg) { 071 072 // Derived from LINPACK code. 073 // Initialize. 074 double[][] A = Arg.getArrayCopy(); 075 m = Arg.getRowDimension(); 076 n = Arg.getColumnDimension(); 077 078 /* Apparently the failing cases are only a proper subset of (m<n), 079 so let's not throw error. Correct fix to come later? 080 if (m<n) { 081 throw new IllegalArgumentException("Jama SVD only works for m >= n"); } 082 */ 083 int nu = Math.min(m,n); 084 s = new double [Math.min(m+1,n)]; 085 U = new double [m][nu]; 086 V = new double [n][n]; 087 double[] e = new double [n]; 088 double[] work = new double [m]; 089 boolean wantu = true; 090 boolean wantv = true; 091 092 // Reduce A to bidiagonal form, storing the diagonal elements 093 // in s and the super-diagonal elements in e. 094 095 int nct = Math.min(m-1,n); 096 int nrt = Math.max(0,Math.min(n-2,m)); 097 for (int k = 0; k < Math.max(nct,nrt); k++) { 098 if (k < nct) { 099 100 // Compute the transformation for the k-th column and 101 // place the k-th diagonal in s[k]. 102 // Compute 2-norm of k-th column without under/overflow. 103 s[k] = 0; 104 for (int i = k; i < m; i++) { 105 s[k] = Maths.hypot(s[k],A[i][k]); 106 } 107 if (s[k] != 0.0) { 108 if (A[k][k] < 0.0) { 109 s[k] = -s[k]; 110 } 111 for (int i = k; i < m; i++) { 112 A[i][k] /= s[k]; 113 } 114 A[k][k] += 1.0; 115 } 116 s[k] = -s[k]; 117 } 118 for (int j = k+1; j < n; j++) { 119 if ((k < nct) && (s[k] != 0.0)) { 120 121 // Apply the transformation. 122 123 double t = 0; 124 for (int i = k; i < m; i++) { 125 t += A[i][k]*A[i][j]; 126 } 127 t = -t/A[k][k]; 128 for (int i = k; i < m; i++) { 129 A[i][j] += t*A[i][k]; 130 } 131 } 132 133 // Place the k-th row of A into e for the 134 // subsequent calculation of the row transformation. 135 136 e[j] = A[k][j]; 137 } 138 if (wantu && (k < nct)) { 139 140 // Place the transformation in U for subsequent back 141 // multiplication. 142 143 for (int i = k; i < m; i++) { 144 U[i][k] = A[i][k]; 145 } 146 } 147 if (k < nrt) { 148 149 // Compute the k-th row transformation and place the 150 // k-th super-diagonal in e[k]. 151 // Compute 2-norm without under/overflow. 152 e[k] = 0; 153 for (int i = k+1; i < n; i++) { 154 e[k] = Maths.hypot(e[k],e[i]); 155 } 156 if (e[k] != 0.0) { 157 if (e[k+1] < 0.0) { 158 e[k] = -e[k]; 159 } 160 for (int i = k+1; i < n; i++) { 161 e[i] /= e[k]; 162 } 163 e[k+1] += 1.0; 164 } 165 e[k] = -e[k]; 166 if ((k+1 < m) && (e[k] != 0.0)) { 167 168 // Apply the transformation. 169 170 for (int i = k+1; i < m; i++) { 171 work[i] = 0.0; 172 } 173 for (int j = k+1; j < n; j++) { 174 for (int i = k+1; i < m; i++) { 175 work[i] += e[j]*A[i][j]; 176 } 177 } 178 for (int j = k+1; j < n; j++) { 179 double t = -e[j]/e[k+1]; 180 for (int i = k+1; i < m; i++) { 181 A[i][j] += t*work[i]; 182 } 183 } 184 } 185 if (wantv) { 186 187 // Place the transformation in V for subsequent 188 // back multiplication. 189 190 for (int i = k+1; i < n; i++) { 191 V[i][k] = e[i]; 192 } 193 } 194 } 195 } 196 197 // Set up the final bidiagonal matrix or order p. 198 199 int p = Math.min(n,m+1); 200 if (nct < n) { 201 s[nct] = A[nct][nct]; 202 } 203 if (m < p) { 204 s[p-1] = 0.0; 205 } 206 if (nrt+1 < p) { 207 e[nrt] = A[nrt][p-1]; 208 } 209 e[p-1] = 0.0; 210 211 // If required, generate U. 212 213 if (wantu) { 214 for (int j = nct; j < nu; j++) { 215 for (int i = 0; i < m; i++) { 216 U[i][j] = 0.0; 217 } 218 U[j][j] = 1.0; 219 } 220 for (int k = nct-1; k >= 0; k--) { 221 if (s[k] != 0.0) { 222 for (int j = k+1; j < nu; j++) { 223 double t = 0; 224 for (int i = k; i < m; i++) { 225 t += U[i][k]*U[i][j]; 226 } 227 t = -t/U[k][k]; 228 for (int i = k; i < m; i++) { 229 U[i][j] += t*U[i][k]; 230 } 231 } 232 for (int i = k; i < m; i++ ) { 233 U[i][k] = -U[i][k]; 234 } 235 U[k][k] = 1.0 + U[k][k]; 236 for (int i = 0; i < k-1; i++) { 237 U[i][k] = 0.0; 238 } 239 } else { 240 for (int i = 0; i < m; i++) { 241 U[i][k] = 0.0; 242 } 243 U[k][k] = 1.0; 244 } 245 } 246 } 247 248 // If required, generate V. 249 250 if (wantv) { 251 for (int k = n-1; k >= 0; k--) { 252 if ((k < nrt) && (e[k] != 0.0)) { 253 for (int j = k+1; j < nu; j++) { 254 double t = 0; 255 for (int i = k+1; i < n; i++) { 256 t += V[i][k]*V[i][j]; 257 } 258 t = -t/V[k+1][k]; 259 for (int i = k+1; i < n; i++) { 260 V[i][j] += t*V[i][k]; 261 } 262 } 263 } 264 for (int i = 0; i < n; i++) { 265 V[i][k] = 0.0; 266 } 267 V[k][k] = 1.0; 268 } 269 } 270 271 // Main iteration loop for the singular values. 272 273 int pp = p-1; 274 int iter = 0; 275 double eps = Math.pow(2.0,-52.0); 276 double tiny = Math.pow(2.0,-966.0); 277 while (p > 0) { 278 int k,kase; 279 280 // Here is where a test for too many iterations would go. 281 282 // This section of the program inspects for 283 // negligible elements in the s and e arrays. On 284 // completion the variables kase and k are set as follows. 285 286 // kase = 1 if s(p) and e[k-1] are negligible and k<p 287 // kase = 2 if s(k) is negligible and k<p 288 // kase = 3 if e[k-1] is negligible, k<p, and 289 // s(k), ..., s(p) are not negligible (qr step). 290 // kase = 4 if e(p-1) is negligible (convergence). 291 292 for (k = p-2; k >= -1; k--) { 293 if (k == -1) { 294 break; 295 } 296 if (Math.abs(e[k]) <= 297 tiny + eps*(Math.abs(s[k]) + Math.abs(s[k+1]))) { 298 e[k] = 0.0; 299 break; 300 } 301 } 302 if (k == p-2) { 303 kase = 4; 304 } else { 305 int ks; 306 for (ks = p-1; ks >= k; ks--) { 307 if (ks == k) { 308 break; 309 } 310 double t = (ks != p ? Math.abs(e[ks]) : 0.) + 311 (ks != k+1 ? Math.abs(e[ks-1]) : 0.); 312 if (Math.abs(s[ks]) <= tiny + eps*t) { 313 s[ks] = 0.0; 314 break; 315 } 316 } 317 if (ks == k) { 318 kase = 3; 319 } else if (ks == p-1) { 320 kase = 1; 321 } else { 322 kase = 2; 323 k = ks; 324 } 325 } 326 k++; 327 328 // Perform the task indicated by kase. 329 330 switch (kase) { 331 332 // Deflate negligible s(p). 333 334 case 1: { 335 double f = e[p-2]; 336 e[p-2] = 0.0; 337 for (int j = p-2; j >= k; j--) { 338 double t = Maths.hypot(s[j],f); 339 double cs = s[j]/t; 340 double sn = f/t; 341 s[j] = t; 342 if (j != k) { 343 f = -sn*e[j-1]; 344 e[j-1] = cs*e[j-1]; 345 } 346 if (wantv) { 347 for (int i = 0; i < n; i++) { 348 t = cs*V[i][j] + sn*V[i][p-1]; 349 V[i][p-1] = -sn*V[i][j] + cs*V[i][p-1]; 350 V[i][j] = t; 351 } 352 } 353 } 354 } 355 break; 356 357 // Split at negligible s(k). 358 359 case 2: { 360 double f = e[k-1]; 361 e[k-1] = 0.0; 362 for (int j = k; j < p; j++) { 363 double t = Maths.hypot(s[j],f); 364 double cs = s[j]/t; 365 double sn = f/t; 366 s[j] = t; 367 f = -sn*e[j]; 368 e[j] = cs*e[j]; 369 if (wantu) { 370 for (int i = 0; i < m; i++) { 371 t = cs*U[i][j] + sn*U[i][k-1]; 372 U[i][k-1] = -sn*U[i][j] + cs*U[i][k-1]; 373 U[i][j] = t; 374 } 375 } 376 } 377 } 378 break; 379 380 // Perform one qr step. 381 382 case 3: { 383 384 // Calculate the shift. 385 386 double scale = Math.max(Math.max(Math.max(Math.max( 387 Math.abs(s[p-1]),Math.abs(s[p-2])),Math.abs(e[p-2])), 388 Math.abs(s[k])),Math.abs(e[k])); 389 double sp = s[p-1]/scale; 390 double spm1 = s[p-2]/scale; 391 double epm1 = e[p-2]/scale; 392 double sk = s[k]/scale; 393 double ek = e[k]/scale; 394 double b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0; 395 double c = (sp*epm1)*(sp*epm1); 396 double shift = 0.0; 397 if ((b != 0.0) || (c != 0.0)) { 398 shift = Math.sqrt(b*b + c); 399 if (b < 0.0) { 400 shift = -shift; 401 } 402 shift = c/(b + shift); 403 } 404 double f = (sk + sp)*(sk - sp) + shift; 405 double g = sk*ek; 406 407 // Chase zeros. 408 409 for (int j = k; j < p-1; j++) { 410 double t = Maths.hypot(f,g); 411 double cs = f/t; 412 double sn = g/t; 413 if (j != k) { 414 e[j-1] = t; 415 } 416 f = cs*s[j] + sn*e[j]; 417 e[j] = cs*e[j] - sn*s[j]; 418 g = sn*s[j+1]; 419 s[j+1] = cs*s[j+1]; 420 if (wantv) { 421 for (int i = 0; i < n; i++) { 422 t = cs*V[i][j] + sn*V[i][j+1]; 423 V[i][j+1] = -sn*V[i][j] + cs*V[i][j+1]; 424 V[i][j] = t; 425 } 426 } 427 t = Maths.hypot(f,g); 428 cs = f/t; 429 sn = g/t; 430 s[j] = t; 431 f = cs*e[j] + sn*s[j+1]; 432 s[j+1] = -sn*e[j] + cs*s[j+1]; 433 g = sn*e[j+1]; 434 e[j+1] = cs*e[j+1]; 435 if (wantu && (j < m-1)) { 436 for (int i = 0; i < m; i++) { 437 t = cs*U[i][j] + sn*U[i][j+1]; 438 U[i][j+1] = -sn*U[i][j] + cs*U[i][j+1]; 439 U[i][j] = t; 440 } 441 } 442 } 443 e[p-2] = f; 444 iter = iter + 1; 445 } 446 break; 447 448 // Convergence. 449 450 case 4: { 451 452 // Make the singular values positive. 453 454 if (s[k] <= 0.0) { 455 s[k] = (s[k] < 0.0 ? -s[k] : 0.0); 456 if (wantv) { 457 for (int i = 0; i <= pp; i++) { 458 V[i][k] = -V[i][k]; 459 } 460 } 461 } 462 463 // Order the singular values. 464 465 while (k < pp) { 466 if (s[k] >= s[k+1]) { 467 break; 468 } 469 double t = s[k]; 470 s[k] = s[k+1]; 471 s[k+1] = t; 472 if (wantv && (k < n-1)) { 473 for (int i = 0; i < n; i++) { 474 t = V[i][k+1]; V[i][k+1] = V[i][k]; V[i][k] = t; 475 } 476 } 477 if (wantu && (k < m-1)) { 478 for (int i = 0; i < m; i++) { 479 t = U[i][k+1]; U[i][k+1] = U[i][k]; U[i][k] = t; 480 } 481 } 482 k++; 483 } 484 iter = 0; 485 p--; 486 } 487 break; 488 } 489 } 490 } 491 492/* ------------------------ 493 Public Methods 494 * ------------------------ */ 495 496 /** Return the left singular vectors 497 @return U 498 */ 499 500 public Matrix getU () { 501 return new Matrix(U,m,Math.min(m+1,n)); 502 } 503 504 /** Return the right singular vectors 505 @return V 506 */ 507 508 public Matrix getV () { 509 return new Matrix(V,n,n); 510 } 511 512 /** Return the one-dimensional array of singular values 513 @return diagonal of S. 514 */ 515 516 public double[] getSingularValues () { 517 return s; 518 } 519 520 /** Return the diagonal matrix of singular values 521 @return S 522 */ 523 524 public Matrix getS () { 525 Matrix X = new Matrix(n,n); 526 double[][] S = X.getArray(); 527 for (int i = 0; i < n; i++) { 528 for (int j = 0; j < n; j++) { 529 S[i][j] = 0.0; 530 } 531 S[i][i] = this.s[i]; 532 } 533 return X; 534 } 535 536 /** Two norm 537 @return max(S) 538 */ 539 540 public double norm2 () { 541 return s[0]; 542 } 543 544 /** Two norm condition number 545 @return max(S)/min(S) 546 */ 547 548 public double cond () { 549 return s[0]/s[Math.min(m,n)-1]; 550 } 551 552 /** Effective numerical matrix rank 553 @return Number of nonnegligible singular values. 554 */ 555 556 public int rank () { 557 double eps = Math.pow(2.0,-52.0); 558 double tol = Math.max(m,n)*s[0]*eps; 559 int r = 0; 560 for (int i = 0; i < s.length; i++) { 561 if (s[i] > tol) { 562 r++; 563 } 564 } 565 return r; 566 } 567}