001 /* 002 * Copyright (C) 2008-2010 by Holger Arndt 003 * 004 * This file is part of the Universal Java Matrix Package (UJMP). 005 * See the NOTICE file distributed with this work for additional 006 * information regarding copyright ownership and licensing. 007 * 008 * UJMP is free software; you can redistribute it and/or modify 009 * it under the terms of the GNU Lesser General Public License as 010 * published by the Free Software Foundation; either version 2 011 * of the License, or (at your option) any later version. 012 * 013 * UJMP is distributed in the hope that it will be useful, 014 * but WITHOUT ANY WARRANTY; without even the implied warranty of 015 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 016 * GNU Lesser General Public License for more details. 017 * 018 * You should have received a copy of the GNU Lesser General Public 019 * License along with UJMP; if not, write to the 020 * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, 021 * Boston, MA 02110-1301 USA 022 */ 023 024 package org.ujmp.core.doublematrix.calculation.general.decomposition; 025 026 import org.ujmp.core.Matrix; 027 import org.ujmp.core.MatrixFactory; 028 import org.ujmp.core.util.DecompositionOps; 029 import org.ujmp.core.util.MathUtil; 030 import org.ujmp.core.util.UJMPSettings; 031 032 /** 033 * Singular Value Decomposition. 034 * <P> 035 * For an m-by-n matrix A, the singular value decomposition is an m-by-(m or n) 036 * orthogonal matrix U, a (m or n)-by-n diagonal matrix S, and an n-by-n 037 * orthogonal matrix V so that A = U*S*V'. 038 * <P> 039 * The singular values, sigma[k] = S[k][k], are ordered so that sigma[0] >= 040 * sigma[1] >= ... >= sigma[n-1]. 041 * <P> 042 * The singular value decompostion always exists, so the constructor will never 043 * fail. The matrix condition number and the effective numerical rank can be 044 * computed from this decomposition. 045 */ 046 047 public interface SVD<T> { 048 049 public T[] calc(T source); 050 051 public static int THRESHOLD = 100; 052 053 public static final SVD<Matrix> MATRIX = new SVD<Matrix>() { 054 055 public final Matrix[] calc(Matrix source) { 056 if (UJMPSettings.getNumberOfThreads() == 1) { 057 if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) { 058 return MATRIXLARGESINGLETHREADED.calc(source); 059 } else { 060 return MATRIXSMALLSINGLETHREADED.calc(source); 061 } 062 } else { 063 if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) { 064 return MATRIXLARGEMULTITHREADED.calc(source); 065 } else { 066 return MATRIXSMALLMULTITHREADED.calc(source); 067 } 068 } 069 } 070 }; 071 072 public static final SVD<Matrix> INSTANCE = MATRIX; 073 074 public static final SVD<Matrix> UJMP = new SVD<Matrix>() { 075 076 public final Matrix[] calc(Matrix source) { 077 SVDMatrix svd = new SVDMatrix(source); 078 return new Matrix[] { svd.getU(), svd.getS(), svd.getV() }; 079 } 080 }; 081 082 public static final SVD<Matrix> MATRIXSMALLSINGLETHREADED = UJMP; 083 084 public static final SVD<Matrix> MATRIXLARGESINGLETHREADED = new SVD<Matrix>() { 085 086 public final Matrix[] calc(Matrix source) { 087 SVD<Matrix> svd = null; 088 if (UJMPSettings.isUseMTJ()) { 089 svd = DecompositionOps.SVD_MTJ; 090 } 091 if (svd == null && UJMPSettings.isUseOjalgo()) { 092 svd = DecompositionOps.SVD_OJALGO; 093 } 094 if (svd == null) { 095 svd = UJMP; 096 } 097 return svd.calc(source); 098 } 099 }; 100 101 public static final SVD<Matrix> MATRIXSMALLMULTITHREADED = UJMP; 102 103 public static final SVD<Matrix> MATRIXLARGEMULTITHREADED = new SVD<Matrix>() { 104 105 public final Matrix[] calc(Matrix source) { 106 SVD<Matrix> svd = null; 107 if (UJMPSettings.isUseOjalgo()) { 108 svd = DecompositionOps.SVD_OJALGO; 109 } 110 if (svd == null && UJMPSettings.isUseMTJ()) { 111 svd = DecompositionOps.SVD_MTJ; 112 } 113 if (svd == null) { 114 svd = UJMP; 115 } 116 return svd.calc(source); 117 } 118 }; 119 120 final class SVDMatrix { 121 private static final double EPSILON = Math.pow(2.0, -52.0); 122 123 private static final double TINY = Math.pow(2.0, -966.0); 124 125 /** 126 * Arrays for internal storage of U and V. 127 * 128 * @serial internal storage of U. 129 * @serial internal storage of V. 130 */ 131 private final double[][] U, V; 132 133 /** 134 * Array for internal storage of singular values. 135 * 136 * @serial internal storage of singular values. 137 */ 138 private final double[] s; 139 140 /** 141 * Row and column dimensions. 142 * 143 * @serial row dimension. 144 * @serial column dimension. 145 * @serial U column dimension. 146 */ 147 private final int m, n, ncu; 148 149 /** 150 * Column specification of matrix U 151 * 152 * @serial U column dimension toggle 153 */ 154 155 private final boolean thin; 156 157 /* 158 * ------------------------ Old Constructor ------------------------ 159 */ 160 /** 161 * Construct the singular value decomposition 162 * 163 * @param Arg 164 * Rectangular matrix 165 * @return Structure to access U, S and V. 166 */ 167 168 public SVDMatrix(Matrix Arg) { 169 this(Arg, true, true, true); 170 } 171 172 /* 173 * ------------------------ Constructor ------------------------ 174 */ 175 176 /** 177 * Construct the singular value decomposition 178 * 179 * @param Arg 180 * Rectangular matrix 181 * @param thin 182 * If true U is economy sized 183 * @param wantu 184 * If true generate the U matrix 185 * @param wantv 186 * If true generate the V matrix 187 * @return Structure to access U, S and V. 188 */ 189 190 public SVDMatrix(Matrix Arg, boolean thin, boolean wantu, boolean wantv) { 191 192 // Derived from LINPACK code. 193 // Initialize. 194 final double[][] A = Arg.toDoubleArray(); 195 m = (int) Arg.getRowCount(); 196 n = (int) Arg.getColumnCount(); 197 this.thin = thin; 198 199 ncu = thin ? Math.min(m, n) : m; 200 s = new double[Math.min(m + 1, n)]; 201 U = new double[m][ncu]; 202 V = new double[n][n]; 203 final double[] e = new double[n]; 204 final double[] work = new double[m]; 205 206 // Reduce A to bidiagonal form, storing the diagonal elements 207 // in s and the super-diagonal elements in e. 208 209 final int nct = Math.min(m - 1, n); 210 final int nrt = Math.max(0, Math.min(n - 2, m)); 211 final int lu = Math.max(nct, nrt); 212 for (int k = 0; k < lu; k++) { 213 if (k < nct) { 214 215 // Compute the transformation for the k-th column and 216 // place the k-th diagonal in s[k]. 217 // Compute 2-norm of k-th column without under/overflow. 218 s[k] = 0; 219 for (int i = k; i < m; i++) { 220 s[k] = MathUtil.hypot(s[k], A[i][k]); 221 } 222 if (s[k] != 0.0) { 223 if (A[k][k] < 0.0) { 224 s[k] = -s[k]; 225 } 226 for (int i = k; i < m; i++) { 227 A[i][k] /= s[k]; 228 } 229 A[k][k] += 1.0; 230 } 231 s[k] = -s[k]; 232 } 233 for (int j = k + 1; j < n; j++) { 234 if ((k < nct) & (s[k] != 0.0)) { 235 236 // Apply the transformation. 237 238 double t = 0; 239 for (int i = k; i < m; i++) { 240 t += A[i][k] * A[i][j]; 241 } 242 t = -t / A[k][k]; 243 for (int i = k; i < m; i++) { 244 A[i][j] += t * A[i][k]; 245 } 246 } 247 248 // Place the k-th row of A into e for the 249 // subsequent calculation of the row transformation. 250 251 e[j] = A[k][j]; 252 } 253 if (wantu & (k < nct)) { 254 255 // Place the transformation in U for subsequent back 256 // multiplication. 257 258 for (int i = k; i < m; i++) { 259 U[i][k] = A[i][k]; 260 } 261 } 262 if (k < nrt) { 263 264 // Compute the k-th row transformation and place the 265 // k-th super-diagonal in e[k]. 266 // Compute 2-norm without under/overflow. 267 e[k] = 0; 268 for (int i = k + 1; i < n; i++) { 269 e[k] = MathUtil.hypot(e[k], e[i]); 270 } 271 if (e[k] != 0.0) { 272 if (e[k + 1] < 0.0) { 273 e[k] = -e[k]; 274 } 275 for (int i = k + 1; i < n; i++) { 276 e[i] /= e[k]; 277 } 278 e[k + 1] += 1.0; 279 } 280 e[k] = -e[k]; 281 if ((k + 1 < m) & (e[k] != 0.0)) { 282 283 // Apply the transformation. 284 285 for (int i = k + 1; i < m; i++) { 286 work[i] = 0.0; 287 } 288 for (int j = k + 1; j < n; j++) { 289 for (int i = k + 1; i < m; i++) { 290 work[i] += e[j] * A[i][j]; 291 } 292 } 293 for (int j = k + 1; j < n; j++) { 294 double t = -e[j] / e[k + 1]; 295 for (int i = k + 1; i < m; i++) { 296 A[i][j] += t * work[i]; 297 } 298 } 299 } 300 if (wantv) { 301 302 // Place the transformation in V for subsequent 303 // back multiplication. 304 305 for (int i = k + 1; i < n; i++) { 306 V[i][k] = e[i]; 307 } 308 } 309 } 310 } 311 312 // Set up the final bidiagonal matrix or order p. 313 int p = Math.min(n, m + 1); 314 if (nct < n) { 315 s[nct] = A[nct][nct]; 316 } 317 if (m < p) { 318 s[p - 1] = 0.0; 319 } 320 if (nrt + 1 < p) { 321 e[nrt] = A[nrt][p - 1]; 322 } 323 e[p - 1] = 0.0; 324 325 // If required, generate U. 326 if (wantu) { 327 for (int j = nct; j < ncu; j++) { 328 for (int i = 0; i < m; i++) { 329 U[i][j] = 0.0; 330 } 331 U[j][j] = 1.0; 332 } 333 for (int k = nct - 1; k >= 0; k--) { 334 if (s[k] != 0.0) { 335 for (int j = k + 1; j < ncu; j++) { 336 double t = 0; 337 for (int i = k; i < m; i++) { 338 t += U[i][k] * U[i][j]; 339 } 340 t = -t / U[k][k]; 341 for (int i = k; i < m; i++) { 342 U[i][j] += t * U[i][k]; 343 } 344 } 345 for (int i = k; i < m; i++) { 346 U[i][k] = -U[i][k]; 347 } 348 U[k][k] += 1.0; 349 for (int i = 0; i < k - 1; i++) { 350 U[i][k] = 0.0; 351 } 352 } else { 353 for (int i = 0; i < m; i++) { 354 U[i][k] = 0.0; 355 } 356 U[k][k] = 1.0; 357 } 358 } 359 } 360 361 // If required, generate V. 362 if (wantv) { 363 for (int k = n - 1; k >= 0; k--) { 364 if ((k < nrt) & (e[k] != 0.0)) { 365 for (int j = k + 1; j < n; j++) { 366 double t = 0; 367 for (int i = k + 1; i < n; i++) { 368 t += V[i][k] * V[i][j]; 369 } 370 t = -t / V[k + 1][k]; 371 for (int i = k + 1; i < n; i++) { 372 V[i][j] += t * V[i][k]; 373 } 374 } 375 } 376 for (int i = 0; i < n; i++) { 377 V[i][k] = 0.0; 378 } 379 V[k][k] = 1.0; 380 } 381 } 382 383 // Main iteration loop for the singular values. 384 385 final int pp = p - 1; 386 int iter = 0; 387 388 while (p > 0) { 389 int k, kase; 390 391 // Here is where a test for too many iterations would go. 392 393 // This section of the program inspects for 394 // negligible elements in the s and e arrays. On 395 // completion the variables kase and k are set as follows. 396 397 // kase = 1 if s(p) and e[k-1] are negligible and k<p 398 // kase = 2 if s(k) is negligible and k<p 399 // kase = 3 if e[k-1] is negligible, k<p, and 400 // s(k), ..., s(p) are not negligible (qr step). 401 // kase = 4 if e(p-1) is negligible (convergence). 402 403 for (k = p - 2; k >= -1; k--) { 404 if (k == -1) { 405 break; 406 } 407 if (Math.abs(e[k]) <= TINY + EPSILON * (Math.abs(s[k]) + Math.abs(s[k + 1]))) { 408 e[k] = 0.0; 409 break; 410 } 411 } 412 if (k == p - 2) { 413 kase = 4; 414 } else { 415 int ks; 416 for (ks = p - 1; ks >= k; ks--) { 417 if (ks == k) { 418 break; 419 } 420 double t = (ks != p ? Math.abs(e[ks]) : 0.) 421 + (ks != k + 1 ? Math.abs(e[ks - 1]) : 0.); 422 if (Math.abs(s[ks]) <= TINY + EPSILON * t) { 423 s[ks] = 0.0; 424 break; 425 } 426 } 427 if (ks == k) { 428 kase = 3; 429 } else if (ks == p - 1) { 430 kase = 1; 431 } else { 432 kase = 2; 433 k = ks; 434 } 435 } 436 k++; 437 438 // Perform the task indicated by kase. 439 440 switch (kase) { 441 442 // Deflate negligible s(p). 443 444 case 1: { 445 double f = e[p - 2]; 446 e[p - 2] = 0.0; 447 for (int j = p - 2; j >= k; j--) { 448 double t = MathUtil.hypot(s[j], f); 449 double cs = s[j] / t; 450 double sn = f / t; 451 s[j] = t; 452 if (j != k) { 453 f = -sn * e[j - 1]; 454 e[j - 1] = cs * e[j - 1]; 455 } 456 if (wantv) { 457 for (int i = 0; i < n; i++) { 458 t = cs * V[i][j] + sn * V[i][p - 1]; 459 V[i][p - 1] = -sn * V[i][j] + cs * V[i][p - 1]; 460 V[i][j] = t; 461 } 462 } 463 } 464 } 465 break; 466 467 // Split at negligible s(k). 468 469 case 2: { 470 double f = e[k - 1]; 471 e[k - 1] = 0.0; 472 for (int j = k; j < p; j++) { 473 double t = MathUtil.hypot(s[j], f); 474 double cs = s[j] / t; 475 double sn = f / t; 476 s[j] = t; 477 f = -sn * e[j]; 478 e[j] = cs * e[j]; 479 if (wantu) { 480 for (int i = 0; i < m; i++) { 481 t = cs * U[i][j] + sn * U[i][k - 1]; 482 U[i][k - 1] = -sn * U[i][j] + cs * U[i][k - 1]; 483 U[i][j] = t; 484 } 485 } 486 } 487 } 488 break; 489 490 // Perform one qr step. 491 492 case 3: { 493 494 // Calculate the shift. 495 496 final double scale = Math.max(Math.max(Math.max(Math.max(Math.abs(s[p - 1]), 497 Math.abs(s[p - 2])), Math.abs(e[p - 2])), Math.abs(s[k])), Math 498 .abs(e[k])); 499 final double sp = s[p - 1] / scale; 500 final double spm1 = s[p - 2] / scale; 501 final double epm1 = e[p - 2] / scale; 502 final double sk = s[k] / scale; 503 final double ek = e[k] / scale; 504 final double b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / 2.0; 505 final double c = (sp * epm1) * (sp * epm1); 506 double shift = 0.0; 507 if ((b != 0.0) | (c != 0.0)) { 508 shift = Math.sqrt(b * b + c); 509 if (b < 0.0) { 510 shift = -shift; 511 } 512 shift = c / (b + shift); 513 } 514 double f = (sk + sp) * (sk - sp) + shift; 515 double g = sk * ek; 516 517 // Chase zeros. 518 519 for (int j = k; j < p - 1; j++) { 520 double t = MathUtil.hypot(f, g); 521 double cs = f / t; 522 double sn = g / t; 523 if (j != k) { 524 e[j - 1] = t; 525 } 526 f = cs * s[j] + sn * e[j]; 527 e[j] = cs * e[j] - sn * s[j]; 528 g = sn * s[j + 1]; 529 s[j + 1] = cs * s[j + 1]; 530 if (wantv) { 531 for (int i = 0; i < n; i++) { 532 t = cs * V[i][j] + sn * V[i][j + 1]; 533 V[i][j + 1] = -sn * V[i][j] + cs * V[i][j + 1]; 534 V[i][j] = t; 535 } 536 } 537 t = MathUtil.hypot(f, g); 538 cs = f / t; 539 sn = g / t; 540 s[j] = t; 541 f = cs * e[j] + sn * s[j + 1]; 542 s[j + 1] = -sn * e[j] + cs * s[j + 1]; 543 g = sn * e[j + 1]; 544 e[j + 1] = cs * e[j + 1]; 545 if (wantu && (j < m - 1)) { 546 for (int i = 0; i < m; i++) { 547 t = cs * U[i][j] + sn * U[i][j + 1]; 548 U[i][j + 1] = -sn * U[i][j] + cs * U[i][j + 1]; 549 U[i][j] = t; 550 } 551 } 552 } 553 e[p - 2] = f; 554 iter++; 555 } 556 break; 557 558 // Convergence. 559 560 case 4: { 561 562 // Make the singular values positive. 563 564 if (s[k] <= 0.0) { 565 s[k] = (s[k] < 0.0 ? -s[k] : 0.0); 566 if (wantv) { 567 for (int i = 0; i < n; i++) { 568 V[i][k] = -V[i][k]; 569 } 570 } 571 } 572 573 // Order the singular values. 574 575 while (k < pp) { 576 if (s[k] >= s[k + 1]) { 577 break; 578 } 579 double t = s[k]; 580 s[k] = s[k + 1]; 581 s[k + 1] = t; 582 if (wantv && (k < n - 1)) { 583 for (int i = 0; i < n; i++) { 584 t = V[i][k + 1]; 585 V[i][k + 1] = V[i][k]; 586 V[i][k] = t; 587 } 588 } 589 if (wantu && (k < m - 1)) { 590 for (int i = 0; i < m; i++) { 591 t = U[i][k + 1]; 592 U[i][k + 1] = U[i][k]; 593 U[i][k] = t; 594 } 595 } 596 k++; 597 } 598 iter = 0; 599 p--; 600 } 601 break; 602 } 603 } 604 } 605 606 /* 607 * ------------------------ Public Methods ------------------------ 608 */ 609 610 /** 611 * Return the left singular vectors 612 * 613 * @return U 614 */ 615 616 public final Matrix getU() { 617 final double[][] x = new double[m][m >= n ? (thin ? Math.min(m + 1, n) : ncu) : ncu]; 618 619 for (int r = 0; r < m; r++) { 620 for (int c = x[0].length; --c >= 0;) { 621 x[r][c] = U[r][c]; 622 } 623 } 624 625 return MatrixFactory.linkToArray(x); 626 } 627 628 /** 629 * Return the right singular vectors 630 * 631 * @return V 632 */ 633 634 public final Matrix getV() { 635 return V == null ? null : MatrixFactory.linkToArray(V); 636 } 637 638 /** 639 * Return the one-dimensional array of singular values 640 * 641 * @return diagonal of S. 642 */ 643 644 public final double[] getSingularValues() { 645 return s; 646 } 647 648 /** 649 * Return the diagonal matrix of singular values 650 * 651 * @return S 652 */ 653 654 public final Matrix getS() { 655 final double[][] X = new double[m >= n ? (thin ? n : ncu) : ncu][n]; 656 for (int i = Math.min(m, n); --i >= 0;) 657 X[i][i] = s[i]; 658 return MatrixFactory.linkToArray(X); 659 } 660 661 /** 662 * Return the diagonal matrix of the reciprocals of the singular values 663 * 664 * @return S+ 665 */ 666 667 public final Matrix getreciprocalS() { 668 final double[][] X = new double[n][m >= n ? (thin ? n : ncu) : ncu]; 669 for (int i = Math.min(m, n) - 1; i >= 0; i--) 670 X[i][i] = s[i] == 0.0 ? 0.0 : 1.0 / s[i]; 671 return MatrixFactory.linkToArray(X); 672 } 673 674 /** 675 * Return the Moore-Penrose (generalized) inverse Slightly modified 676 * version of Kim van der Linde's code 677 * 678 * @param omit 679 * if true tolerance based omitting of negligible singular 680 * values 681 * @return A+ 682 */ 683 684 public final Matrix inverse(boolean omit) { 685 final double[][] inverse = new double[n][m]; 686 if (rank() > 0) { 687 final double[] reciprocalS = new double[s.length]; 688 if (omit) { 689 double tol = Math.max(m, n) * s[0] * EPSILON; 690 for (int i = s.length - 1; i >= 0; i--) 691 reciprocalS[i] = Math.abs(s[i]) < tol ? 0.0 : 1.0 / s[i]; 692 } else 693 for (int i = s.length - 1; i >= 0; i--) 694 reciprocalS[i] = s[i] == 0.0 ? 0.0 : 1.0 / s[i]; 695 int min = Math.min(n, ncu); 696 for (int i = n - 1; i >= 0; i--) 697 for (int j = m - 1; j >= 0; j--) 698 for (int k = min - 1; k >= 0; k--) 699 inverse[i][j] += V[i][k] * reciprocalS[k] * U[j][k]; 700 } 701 return MatrixFactory.linkToArray(inverse); 702 } 703 704 /** 705 * Two norm 706 * 707 * @return max(S) 708 */ 709 710 public final double norm2() { 711 return s[0]; 712 } 713 714 /** 715 * Two norm condition number 716 * 717 * @return max(S)/min(S) 718 */ 719 720 public final double cond() { 721 return s[0] / s[Math.min(m, n) - 1]; 722 } 723 724 /** 725 * Effective numerical matrix rank 726 * 727 * @return Number of nonnegligible singular values. 728 */ 729 730 public final int rank() { 731 final double tol = Math.max(m, n) * s[0] * EPSILON; 732 int r = 0; 733 for (int i = 0; i < s.length; i++) { 734 if (s[i] > tol) { 735 r++; 736 } 737 } 738 return r; 739 } 740 741 } 742 }