001/** 002 * Copyright (c) 2011, The University of Southampton and the individual contributors. 003 * All rights reserved. 004 * 005 * Redistribution and use in source and binary forms, with or without modification, 006 * are permitted provided that the following conditions are met: 007 * 008 * * Redistributions of source code must retain the above copyright notice, 009 * this list of conditions and the following disclaimer. 010 * 011 * * Redistributions in binary form must reproduce the above copyright notice, 012 * this list of conditions and the following disclaimer in the documentation 013 * and/or other materials provided with the distribution. 014 * 015 * * Neither the name of the University of Southampton nor the names of its 016 * contributors may be used to endorse or promote products derived from this 017 * software without specific prior written permission. 018 * 019 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND 020 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 021 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 022 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR 023 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES 024 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 025 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON 026 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 027 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 028 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 029 */ 030package org.openimaj.math.matrix; 031 032import java.io.PrintWriter; 033import java.io.StringWriter; 034import java.util.Random; 035 036import Jama.EigenvalueDecomposition; 037import Jama.Matrix; 038import ch.akuhn.matrix.SparseMatrix; 039import no.uib.cipr.matrix.DenseMatrix; 040import no.uib.cipr.matrix.NotConvergedException; 041import no.uib.cipr.matrix.SVD; 042 043/** 044 * Miscellaneous matrix operations. 045 * 046 * @author Sina Samangooei (ss@ecs.soton.ac.uk) 047 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk) 048 * 049 */ 050public class MatrixUtils { 051 private MatrixUtils() { 052 } 053 054 /** 055 * Are any values NaN or Inf? 056 * 057 * @param matrix 058 * matrix to test 059 * @return true if any elements are NaN or Inf; false otherwise 060 */ 061 public static boolean anyNaNorInf(Matrix matrix) { 062 for (final double[] arrLine : matrix.getArray()) { 063 for (final double d : arrLine) { 064 if (Double.isNaN(d) || Double.isInfinite(d)) 065 return true; 066 } 067 } 068 return false; 069 } 070 071 /** 072 * Get the maximum absolute value of the diagonal. 073 * 074 * @param matrix 075 * the matrix 076 * @return the maximum absolute value 077 */ 078 public static double maxAbsDiag(Matrix matrix) { 079 double max = -1; 080 081 for (int i = 0; i < matrix.getColumnDimension(); i++) { 082 final double curr = Math.abs(matrix.get(i, i)); 083 if (max < curr) { 084 max = curr; 085 } 086 } 087 return max; 088 } 089 090 /** 091 * Get the minimum absolute value of the diagonal. 092 * 093 * @param matrix 094 * the matrix 095 * @return the minimum absolute value 096 */ 097 public static double minAbsDiag(Matrix matrix) { 098 double min = Double.MAX_VALUE; 099 100 for (int i = 0; i < matrix.getColumnDimension(); i++) { 101 final double curr = Math.abs(matrix.get(i, i)); 102 if (min > curr) { 103 min = curr; 104 } 105 } 106 return min; 107 } 108 109 /** 110 * Compute the principle square root, X, of the matrix A such that A=X*X 111 * 112 * @param matrix 113 * the matrix 114 * @return the sqrt of the matrix 115 */ 116 public static Matrix sqrt(Matrix matrix) { 117 // A = V*D*V' 118 final EigenvalueDecomposition evd = matrix.eig(); 119 final Matrix v = evd.getV(); 120 final Matrix d = evd.getD(); 121 122 // sqrt of cells of D and store in-place 123 for (int r = 0; r < d.getRowDimension(); r++) 124 for (int c = 0; c < d.getColumnDimension(); c++) 125 d.set(r, c, Math.sqrt(d.get(r, c))); 126 127 // Y = V*D/V 128 // Y = V'.solve(V*D)' 129 final Matrix a = v.inverse(); 130 final Matrix b = v.times(d).inverse(); 131 return a.solve(b).inverse(); 132 } 133 134 /** 135 * Computes the Moore-Penrose pseudoinverse. This is just a convenience wrapper 136 * around {@link PseudoInverse#pseudoInverse(Matrix)}. 137 * 138 * @param matrix 139 * the matrix to invert. 140 * @return the inverted matrix 141 * @see PseudoInverse#pseudoInverse(Matrix) 142 */ 143 public static Matrix pseudoInverse(Matrix matrix) { 144 return PseudoInverse.pseudoInverse(matrix); 145 } 146 147 /** 148 * Compute the inverse square root, X, of the symmetric matrix A; A^-(1/2) 149 * 150 * @param matrix 151 * the symmetric matrix 152 * @return the inverse sqrt of the matrix 153 */ 154 public static Matrix invSqrtSym(Matrix matrix) { 155 // A = V*D*V' 156 final EigenvalueDecomposition evd = matrix.eig(); 157 final Matrix v = evd.getV(); 158 final Matrix d = evd.getD(); 159 160 // sqrt of cells of D and store in-place 161 for (int r = 0; r < d.getRowDimension(); r++) { 162 for (int c = 0; c < d.getColumnDimension(); c++) { 163 if (d.get(r, c) > 0) 164 d.set(r, c, 1 / Math.sqrt(d.get(r, c))); 165 else 166 d.set(r, c, 0); 167 } 168 } 169 170 return v.times(d).times(v.transpose()); 171 } 172 173 /** 174 * Return a copy of the input matrix with all elements set to their absolute 175 * value. 176 * 177 * @param mat 178 * the matrix. 179 * @return the absolute matrix. 180 */ 181 public static Matrix abs(Matrix mat) { 182 final Matrix copy = mat.copy(); 183 for (int i = 0; i < mat.getRowDimension(); i++) { 184 for (int j = 0; j < mat.getColumnDimension(); j++) { 185 copy.set(i, j, Math.abs(mat.get(i, j))); 186 } 187 } 188 return copy; 189 } 190 191 /** 192 * Check if two matrices are equal 193 * 194 * @param m1 195 * first matrix 196 * @param m2 197 * second matrix 198 * @param eps 199 * epsilon for checking values 200 * @return true if matrices have same size and all elements are equal within 201 * eps; false otherwise 202 */ 203 public static boolean equals(Matrix m1, Matrix m2, double eps) { 204 final double[][] a1 = m1.getArray(); 205 final double[][] a2 = m2.getArray(); 206 207 if (a1.length != a2.length || a1[0].length != a2[0].length) 208 return false; 209 210 for (int r = 0; r < a1.length; r++) 211 for (int c = 0; c < a1[r].length; c++) 212 if (Math.abs(a1[r][c] - a2[r][c]) > eps) 213 return false; 214 215 return true; 216 } 217 218 /** 219 * Return a copy of the matrix with all the values raised to a power. 220 * 221 * @param mat 222 * the matrix. 223 * @param exp 224 * the power. 225 * @return a matrix. 226 */ 227 public static Matrix pow(Matrix mat, double exp) { 228 final Matrix copy = mat.copy(); 229 for (int i = 0; i < mat.getRowDimension(); i++) { 230 for (int j = 0; j < mat.getColumnDimension(); j++) { 231 copy.set(i, j, Math.pow(mat.get(i, j), exp)); 232 } 233 } 234 return copy; 235 } 236 237 /** 238 * Generate a {@link String} representation of a matrix. 239 * 240 * @param mat 241 * the matrix 242 * @return a string representation 243 */ 244 public static String toString(Matrix mat) { 245 final StringWriter matWriter = new StringWriter(); 246 mat.print(new PrintWriter(matWriter), 5, 5); 247 return matWriter.getBuffer().toString(); 248 } 249 250 /** 251 * Compute the sum of all elements of the matrix. 252 * 253 * @param mat 254 * the matrix. 255 * @return the sum. 256 */ 257 public static double sum(Matrix mat) { 258 double sum = 0; 259 for (int i = 0; i < mat.getRowDimension(); i++) { 260 for (int j = 0; j < mat.getColumnDimension(); j++) { 261 sum += mat.get(i, j); 262 } 263 } 264 return sum; 265 } 266 267 /** 268 * Zero the matrix 269 * 270 * @param m 271 * the matrix 272 */ 273 public static void zero(Matrix m) { 274 m.timesEquals(0); 275 } 276 277 /** 278 * Compute the real Eigen decomposition of a symmetric 2x2 matrix. Warning: 279 * Doesn't check the size or whether the input is symmetric. 280 * 281 * @param m 282 * the matrix 283 * @return the Eigen vectors and values. 284 */ 285 public static EigenValueVectorPair symmetricEig2x2(Matrix m) { 286 final double a = m.get(0, 0); 287 final double b = m.get(0, 1); 288 final double c = b; 289 final double d = m.get(1, 1); 290 291 final double trace = a + d; 292 final double det = a * d - b * c; 293 294 final Matrix val = new Matrix(2, 2); 295 double sqrtInner = (trace * trace / 4) - det; 296 // FIXME: make it deal with imaginary numbers. 297 if (sqrtInner < 0) { 298 final EigenvalueDecomposition e = m.eig(); 299 return new EigenValueVectorPair(e.getD(), e.getV()); 300 } 301 302 sqrtInner = Math.sqrt(sqrtInner); 303 double firstEig = trace / 2 + sqrtInner; 304 double secondEig = trace / 2 - sqrtInner; 305 if (firstEig > secondEig) { 306 final double tmp = firstEig; 307 firstEig = secondEig; 308 secondEig = tmp; 309 } 310 311 val.set(0, 0, firstEig); 312 val.set(1, 1, secondEig); 313 314 final Matrix vec = new Matrix(2, 2); 315 316 final double v1 = firstEig - a; 317 final double v2 = secondEig - a; 318 final double norm1 = Math.sqrt(v1 * v1 + b * b); 319 final double norm2 = Math.sqrt(b * b + v2 * v2); 320 vec.set(0, 0, b / norm1); 321 vec.set(0, 1, b / norm2); 322 vec.set(1, 0, v1 / norm1); 323 vec.set(1, 1, v2 / norm2); 324 325 // To deal with rounding error 326 vec.set(1, 0, vec.get(0, 1)); 327 328 final EigenValueVectorPair ret = new EigenValueVectorPair(val, vec); 329 return ret; 330 } 331 332 /** 333 * An eigen decomposition that uses a deterministic method if the matrix is 2x2. 334 * 335 * This function returns values as in {@link EigenvalueDecomposition} i.e. the 336 * largest eigen value is held in the [m.rows - 1,m.cols-1] (i.e. [1,1]) 337 * location 338 * 339 * @param m 340 * @return the decomposition 341 */ 342 public static EigenValueVectorPair eig2x2(Matrix m) { 343 if (m.getColumnDimension() != 2 || m.getRowDimension() != 2) { 344 final EigenvalueDecomposition e = m.eig(); 345 return new EigenValueVectorPair(e.getD(), e.getV()); 346 } 347 /** 348 * A = 1 B = a + d C = ad - bc 349 * 350 * x = ( - B (+/-) sqrt(B^2 - 4AC) )/ (2A) 351 */ 352 final double a = m.get(0, 0); 353 final double b = m.get(0, 1); 354 final double c = m.get(1, 0); 355 final double d = m.get(1, 1); 356 357 final double trace = a + d; 358 final double det = a * d - b * c; 359 360 final Matrix val = new Matrix(2, 2); 361 double sqrtInner = (trace * trace / 4) - det; 362 // FIXME: make it deal with imaginary numbers. 363 if (sqrtInner < 0) { 364 final EigenvalueDecomposition e = m.eig(); 365 return new EigenValueVectorPair(e.getD(), e.getV()); 366 } 367 368 sqrtInner = Math.sqrt(sqrtInner); 369 double firstEig = trace / 2 + sqrtInner; 370 double secondEig = trace / 2 - sqrtInner; 371 if (firstEig > secondEig) { 372 final double tmp = firstEig; 373 firstEig = secondEig; 374 secondEig = tmp; 375 } 376 377 val.set(0, 0, firstEig); 378 val.set(1, 1, secondEig); 379 380 final Matrix vec = new Matrix(2, 2); 381 if (b == 0 && c == 0) { 382 vec.set(0, 0, 1); 383 vec.set(1, 1, 1); 384 } else { 385 if (c != 0) { 386 final double v1 = firstEig - d; 387 final double v2 = secondEig - d; 388 final double norm1 = Math.sqrt(v1 * v1 + c * c); 389 final double norm2 = Math.sqrt(c * c + v2 * v2); 390 vec.set(0, 0, v1 / norm1); 391 vec.set(0, 1, v2 / norm2); 392 vec.set(1, 0, c / norm1); 393 vec.set(1, 1, c / norm2); 394 } else if (b != 0) { 395 final double v1 = firstEig - a; 396 final double v2 = secondEig - a; 397 final double norm1 = Math.sqrt(v1 * v1 + b * b); 398 final double norm2 = Math.sqrt(b * b + v2 * v2); 399 vec.set(0, 0, b / norm1); 400 vec.set(0, 1, b / norm2); 401 vec.set(1, 0, v1 / norm1); 402 vec.set(1, 1, v2 / norm2); 403 } 404 } 405 406 final EigenValueVectorPair ret = new EigenValueVectorPair(val, vec); 407 return ret; 408 } 409 410 /** 411 * Construct a matrix from a 2D float array of data. 412 * 413 * @param data 414 * the data. 415 * @return the matrix. 416 */ 417 public static Matrix matrixFromFloat(float[][] data) { 418 final Matrix out = new Matrix(data.length, data[0].length); 419 for (int i = 0; i < data.length; i++) { 420 for (int j = 0; j < data[i].length; j++) { 421 out.set(i, j, data[i][j]); 422 } 423 } 424 return out; 425 } 426 427 /** 428 * Reduce the rank a matrix by estimating a the best (in a least-squares sense) 429 * approximation using the thin SVD. 430 * 431 * @param m 432 * the matrix to reduce. 433 * @param rank 434 * the desired rank. 435 * @return the rank-reduced matrix. 436 */ 437 public static Matrix reduceRank(Matrix m, int rank) { 438 if (rank > Math.min(m.getColumnDimension(), m.getRowDimension())) { 439 return m; 440 } 441 442 final no.uib.cipr.matrix.DenseMatrix mjtA = new no.uib.cipr.matrix.DenseMatrix( 443 m.getArray()); 444 no.uib.cipr.matrix.SVD svd; 445 try { 446 svd = no.uib.cipr.matrix.SVD.factorize(mjtA); 447 } catch (final NotConvergedException e) { 448 throw new RuntimeException(e); 449 } 450 451 final DenseMatrix U = svd.getU(); 452 final DenseMatrix Vt = svd.getVt(); 453 final double[] svector = svd.getS(); 454 final DenseMatrix S = new DenseMatrix(U.numColumns(), Vt.numRows()); 455 for (int i = 0; i < rank; i++) 456 S.set(i, i, svector[i]); 457 458 final DenseMatrix C = new DenseMatrix(U.numRows(), S.numColumns()); 459 final DenseMatrix out = new DenseMatrix(C.numRows(), Vt.numColumns()); 460 U.mult(S, C); 461 C.mult(Vt, out); 462 463 final Matrix outFinal = convert(out); 464 return outFinal; 465 } 466 467 /** 468 * Convert a {@link DenseMatrix} to a {@link Matrix}. 469 * 470 * @param mjt 471 * {@link DenseMatrix} to convert 472 * @return converted matrix. 473 */ 474 public static Matrix convert(DenseMatrix mjt) { 475 return convert(mjt, mjt.numRows(), mjt.numColumns()); 476 } 477 478 /** 479 * Convert a {@link DenseMatrix} to a {@link Matrix}. 480 * 481 * @param mjt 482 * {@link DenseMatrix} to convert 483 * @param nrows 484 * number of rows to copy 485 * @param ncols 486 * number of columns to copy 487 * @return converted matrix. 488 */ 489 public static Matrix convert(DenseMatrix mjt, int nrows, int ncols) { 490 final double[][] d = new double[nrows][ncols]; 491 492 final double[] mjtd = mjt.getData(); 493 for (int r = 0; r < nrows; r++) { 494 for (int c = 0; c < ncols; c++) { 495 d[r][c] = mjtd[r + c * d.length]; 496 } 497 } 498 return new Matrix(d); 499 } 500 501 /** 502 * Create a copy of a matrix with the columns in reverse order. 503 * 504 * @param m 505 * the input matrix 506 * @return a copy with the column order reversed 507 */ 508 public static Matrix reverseColumns(Matrix m) { 509 return reverseColumnsInplace(m.copy()); 510 } 511 512 /** 513 * Reverse the column order of the input matrix inplace. 514 * 515 * @param m 516 * the input matrix 517 * @return the input matrix 518 */ 519 public static Matrix reverseColumnsInplace(Matrix m) { 520 final double[][] data = m.getArray(); 521 final int rows = data.length; 522 final int cols = data[0].length; 523 final int halfCols = cols / 2; 524 525 for (int r = 0; r < rows; r++) { 526 for (int c = 0; c < halfCols; c++) { 527 final double tmp = data[r][c]; 528 data[r][c] = data[r][cols - c - 1]; 529 data[r][cols - c - 1] = tmp; 530 } 531 } 532 533 return m; 534 } 535 536 /** 537 * Create a copy of a matrix with the rows in reverse order. 538 * 539 * @param m 540 * the input matrix 541 * @return a copy with the row order reversed 542 */ 543 public static Matrix reverseRows(Matrix m) { 544 return reverseRowsInplace(m.copy()); 545 } 546 547 /** 548 * Reverse the row order of the input matrix inplace. 549 * 550 * @param m 551 * the input matrix 552 * @return the input matrix 553 */ 554 public static Matrix reverseRowsInplace(Matrix m) { 555 final double[][] data = m.getArray(); 556 final int rows = data.length; 557 final int halfRows = rows / 2; 558 559 for (int r = 0; r < halfRows; r++) { 560 final double[] tmp = data[r]; 561 data[r] = data[rows - r - 1]; 562 data[rows - r - 1] = tmp; 563 } 564 565 return m; 566 } 567 568 /** 569 * Create a diagonal matrix 570 * 571 * @param s 572 * length diagonal numbers 573 * @return new Matrix(s.length,s.length) s.t. diagonal element i,i = s[i] 574 */ 575 public static Matrix diag(double[] s) { 576 final Matrix r = new Matrix(s.length, s.length); 577 for (int i = 0; i < s.length; i++) { 578 r.set(i, i, s[i]); 579 } 580 return r; 581 } 582 583 /** 584 * Set the values of the elements in a single column to a constant value. 585 * 586 * @param m 587 * the matrix 588 * @param c 589 * the column 590 * @param v 591 * the constant value 592 * @return the matrix 593 */ 594 public static Matrix setColumn(Matrix m, int c, double v) { 595 final double[][] data = m.getArray(); 596 final int rows = m.getRowDimension(); 597 598 for (int r = 0; r < rows; r++) 599 data[r][c] = v; 600 601 return m; 602 } 603 604 /** 605 * Set the values of the elements in a single column to a constant value. 606 * 607 * @param m 608 * the matrix 609 * @param r 610 * the row 611 * @param v 612 * the constant value 613 * @return the matrix 614 */ 615 public static Matrix setRow(Matrix m, int r, double v) { 616 final double[][] data = m.getArray(); 617 final int cols = m.getColumnDimension(); 618 619 for (int c = 0; c < cols; c++) 620 data[r][c] = v; 621 622 return m; 623 } 624 625 /** 626 * Fill a matrix with a constant value. 627 * 628 * @param m 629 * the matrix 630 * @param v 631 * the constant value 632 * @return the matrix 633 */ 634 public static Matrix fill(Matrix m, double v) { 635 final double[][] data = m.getArray(); 636 637 final int rows = m.getRowDimension(); 638 final int cols = m.getColumnDimension(); 639 640 for (int r = 0; r < rows; r++) 641 for (int c = 0; c < cols; c++) 642 data[r][c] = v; 643 644 return m; 645 } 646 647 /** 648 * Subtract a constant from all values 649 * 650 * @param m 651 * the matrix 652 * @param v 653 * the constant value 654 * @return the matrix 655 */ 656 public static Matrix minus(Matrix m, double v) { 657 final double[][] data = m.getArray(); 658 659 final int rows = m.getRowDimension(); 660 final int cols = m.getColumnDimension(); 661 662 for (int r = 0; r < rows; r++) 663 for (int c = 0; c < cols; c++) 664 data[r][c] -= v; 665 666 return m; 667 } 668 669 /** 670 * Add a constant to all values 671 * 672 * @param m 673 * the matrix 674 * @param v 675 * the constant value 676 * @return the matrix 677 */ 678 public static Matrix plus(Matrix m, double v) { 679 final double[][] data = m.getArray(); 680 681 final int rows = m.getRowDimension(); 682 final int cols = m.getColumnDimension(); 683 684 for (int r = 0; r < rows; r++) 685 for (int c = 0; c < cols; c++) 686 data[r][c] += v; 687 688 return m; 689 } 690 691 /** 692 * Get a reshaped copy of the input matrix 693 * 694 * @param m 695 * the matrix to reshape 696 * @param newRows 697 * the new number of rows 698 * @return new matrix 699 */ 700 public static Matrix reshape(Matrix m, int newRows) { 701 final int oldCols = m.getColumnDimension(); 702 final int length = oldCols * m.getRowDimension(); 703 final int newCols = length / newRows; 704 final Matrix mat = new Matrix(newRows, newCols); 705 706 final double[][] m1v = m.getArray(); 707 final double[][] m2v = mat.getArray(); 708 709 int r1 = 0, r2 = 0, c1 = 0, c2 = 0; 710 for (int i = 0; i < length; i++) { 711 m2v[r2][c2] = m1v[r1][c1]; 712 713 c1++; 714 if (c1 >= oldCols) { 715 c1 = 0; 716 r1++; 717 } 718 719 c2++; 720 if (c2 >= newCols) { 721 c2 = 0; 722 r2++; 723 } 724 } 725 726 return mat; 727 } 728 729 /** 730 * Get a reshaped copy of the input matrix 731 * 732 * @param m 733 * the matrix to reshape 734 * @param newRows 735 * the new number of rows 736 * @param columnMajor 737 * if true, values are drawn and placed down columns first. if false 738 * values are drawn and placed across rows first 739 * @return new matrix 740 */ 741 public static Matrix reshape(Matrix m, int newRows, boolean columnMajor) { 742 final int oldCols = m.getColumnDimension(); 743 final int oldRows = m.getRowDimension(); 744 final int length = oldCols * m.getRowDimension(); 745 final int newCols = length / newRows; 746 final Matrix mat = new Matrix(newRows, newCols); 747 748 final double[][] m1v = m.getArray(); 749 final double[][] m2v = mat.getArray(); 750 751 int r1 = 0, r2 = 0, c1 = 0, c2 = 0; 752 if (!columnMajor) { 753 for (int i = 0; i < length; i++) { 754 m2v[r2][c2] = m1v[r1][c1]; 755 756 c1++; 757 if (c1 >= oldCols) { 758 c1 = 0; 759 r1++; 760 } 761 762 c2++; 763 if (c2 >= newCols) { 764 c2 = 0; 765 r2++; 766 } 767 } 768 } else { 769 for (int i = 0; i < length; i++) { 770 m2v[r2][c2] = m1v[r1][c1]; 771 772 r1++; 773 if (r1 >= oldRows) { 774 r1 = 0; 775 c1++; 776 } 777 778 r2++; 779 if (r2 >= newRows) { 780 r2 = 0; 781 c2++; 782 } 783 } 784 } 785 786 return mat; 787 } 788 789 /** 790 * Compute the sum of values in a single column 791 * 792 * @param m 793 * the matrix 794 * @param col 795 * the column 796 * @return the sum of values in column col 797 */ 798 public static double sumColumn(Matrix m, int col) { 799 final double[][] data = m.getArray(); 800 final int rows = m.getRowDimension(); 801 802 double sum = 0; 803 for (int r = 0; r < rows; r++) 804 sum += data[r][col]; 805 806 return sum; 807 } 808 809 /** 810 * Compute the sum of values in a single row 811 * 812 * @param m 813 * the matrix 814 * @param row 815 * the row 816 * @return the sum of values in row row 817 */ 818 public static double sumRow(Matrix m, int row) { 819 final double[][] data = m.getArray(); 820 final int cols = m.getColumnDimension(); 821 822 double sum = 0; 823 for (int c = 0; c < cols; c++) 824 sum += data[row][c]; 825 826 return sum; 827 } 828 829 /** 830 * Increment values in a single column by a constant 831 * 832 * @param m 833 * the matrix 834 * @param col 835 * the column 836 * @param value 837 * the constant 838 * @return the matrix 839 */ 840 public static Matrix incrColumn(Matrix m, int col, double value) { 841 final double[][] data = m.getArray(); 842 final int rows = m.getRowDimension(); 843 844 for (int r = 0; r < rows; r++) 845 data[r][col] += value; 846 847 return m; 848 } 849 850 /** 851 * Increment values in a single column by a constant 852 * 853 * @param m 854 * the matrix 855 * @param row 856 * the row 857 * @param value 858 * the constant 859 * @return the sum of values in row row 860 */ 861 public static Matrix incrRow(Matrix m, int row, double value) { 862 final double[][] data = m.getArray(); 863 final int cols = m.getColumnDimension(); 864 865 for (int c = 0; c < cols; c++) 866 data[row][c] += value; 867 868 return m; 869 } 870 871 /** 872 * round (using {@link Math#round(double)} each value of the matrix 873 * 874 * @param times 875 * @return same matrix as handed in 876 */ 877 public static Matrix round(Matrix times) { 878 final double[][] data = times.getArray(); 879 for (int i = 0; i < data.length; i++) { 880 for (int j = 0; j < data[i].length; j++) { 881 data[i][j] = Math.round(data[i][j]); 882 } 883 } 884 return times; 885 } 886 887 /** 888 * min(A,B) returns an array the same size as A and B with the smallest elements 889 * taken from A or B. The dimensions of A and B must match 890 * 891 * @param A 892 * @param B 893 * @return new Matrix filled with min from A and B 894 */ 895 public static Matrix min(Matrix A, Matrix B) { 896 final double[][] dataA = A.getArray(); 897 final double[][] dataB = B.getArray(); 898 final Matrix ret = A.copy(); 899 final double[][] dataRet = ret.getArray(); 900 for (int i = 0; i < dataA.length; i++) { 901 for (int j = 0; j < dataB[i].length; j++) { 902 dataRet[i][j] = Math.min(dataA[i][j], dataB[i][j]); 903 } 904 } 905 return ret; 906 } 907 908 /** 909 * d to the power of each value in range. as with {@link #range} range is: - a 910 * single number (a) (0:1:a) - two numbers (a,b) (a:1:b) - three numbers (a,b,c) 911 * (a:b:c) 912 * 913 * any other amount of range results in a {@link RuntimeException} 914 * 915 * @param d 916 * @param range 917 * @return d to the power of each value in range 918 */ 919 public static Matrix rangePow(double d, double... range) { 920 double start, end, delta; 921 if (range.length == 1) { 922 start = 0; 923 end = range[0]; 924 delta = 1; 925 } else if (range.length == 2) { 926 start = range[0]; 927 end = range[1]; 928 delta = 1; 929 } else if (range.length == 3) { 930 start = range[0]; 931 end = range[2]; 932 delta = range[1]; 933 } else { 934 throw new RuntimeException("Invalid range options selected"); 935 } 936 final int l = (int) ((end - start + 1) / delta); 937 final double[][] out = new double[1][l]; 938 for (int i = 0; i < l; i++) { 939 out[0][i] = Math.pow(d, start + (i * delta)); 940 } 941 return new Matrix(out); 942 } 943 944 /** 945 * range is: - a single number (a) (0:1:a) - two numbers (a,b) (a:1:b) - three 946 * numbers (a,b,c) (a:b:c) 947 * 948 * @param range 949 * @return the range defined 950 */ 951 public static Matrix range(double... range) { 952 double start, end, delta; 953 if (range.length == 1) { 954 start = 0; 955 end = range[0]; 956 delta = 1; 957 } else if (range.length == 2) { 958 start = range[0]; 959 end = range[1]; 960 delta = 1; 961 } else if (range.length == 3) { 962 start = range[0]; 963 end = range[2]; 964 delta = range[1]; 965 } else { 966 throw new RuntimeException("Invalid range options selected"); 967 } 968 final int l = (int) Math.floor((end - start) / delta) + 1; 969 final double[][] out = new double[1][l]; 970 for (int i = 0; i < l; i++) { 971 out[0][i] = start + (i * delta); 972 } 973 return new Matrix(out); 974 } 975 976 /** 977 * range is: - a single number (a) (0:1:a) - two numbers (a,b) (a:1:b) - three 978 * numbers (a,b,c) (a:b:c) 979 * 980 * @param range 981 * @return the range defined 982 */ 983 public static Matrix range(int... range) { 984 int start, end, delta; 985 if (range.length == 1) { 986 start = 0; 987 end = range[0]; 988 delta = 1; 989 } else if (range.length == 2) { 990 start = range[0]; 991 end = range[1]; 992 delta = 1; 993 } else if (range.length == 3) { 994 start = range[0]; 995 end = range[2]; 996 delta = range[1]; 997 } else { 998 throw new RuntimeException("Invalid range options selected"); 999 } 1000 final int l = (int) Math.floor((end - start) / delta) + 1; 1001 final double[][] out = new double[1][l]; 1002 for (int i = 0; i < l; i++) { 1003 out[0][i] = start + (i * delta); 1004 } 1005 return new Matrix(out); 1006 } 1007 1008 /** 1009 * Given two row vectors, construct the power set of rowvector combinations 1010 * 1011 * @param A 1012 * @param B 1013 * @return a new matrix of size A.cols * B.cols 1014 */ 1015 public static Matrix ntuples(Matrix A, Matrix B) { 1016 final double[][] Adata = A.getArray(); 1017 final double[][] Bdata = B.getArray(); 1018 1019 final double[][] out = new double[2][Adata[0].length * Bdata[0].length]; 1020 int i = 0; 1021 for (final double a : Adata[0]) { 1022 for (final double b : Bdata[0]) { 1023 out[0][i] = a; 1024 out[1][i] = b; 1025 i++; 1026 } 1027 } 1028 return new Matrix(out); 1029 } 1030 1031 /** 1032 * Given a matrix, repeat the matrix over i rows and j columns 1033 * 1034 * @param x 1035 * @param i 1036 * @param j 1037 * @return repeated matrix 1038 */ 1039 public static Matrix repmat(Matrix x, int i, int j) { 1040 final double[][] xdata = x.getArray(); 1041 final double[][] newmat = new double[xdata.length * i][xdata[0].length * j]; 1042 for (int k = 0; k < newmat.length; k += xdata.length) { 1043 for (int l = 0; l < newmat[0].length; l += xdata[0].length) { 1044 int rowcopyindex = 0; 1045 for (final double[] ds : xdata) { 1046 System.arraycopy(ds, 0, newmat[k + rowcopyindex], l, xdata[0].length); 1047 rowcopyindex += 1; 1048 } 1049 } 1050 } 1051 return new Matrix(newmat); 1052 } 1053 1054 /** 1055 * horizontally stack all the matrices provided. i.e. ret = [x1 x2 x3 x4 ... xn] 1056 * 1057 * @param x 1058 * @return horizontally stacked 1059 */ 1060 public static Matrix hstack(Matrix... x) { 1061 final int height = x[0].getRowDimension(); 1062 int width = 0; 1063 for (final Matrix matrix : x) { 1064 width += matrix.getColumnDimension(); 1065 } 1066 final double[][] newmat = new double[height][width]; 1067 int colindex = 0; 1068 for (final Matrix matrix : x) { 1069 final double[][] matdata = matrix.getArray(); 1070 final int w = matrix.getColumnDimension(); 1071 for (int i = 0; i < height; i++) { 1072 System.arraycopy(matdata[i], 0, newmat[i], colindex, w); 1073 } 1074 colindex += w; 1075 } 1076 return new Matrix(newmat); 1077 } 1078 1079 /** 1080 * Add the rows to the mat at rowIndex. Assumes MANY things with no checks: 1081 * rows.rows == rowIndex.length mat.cols == rows.cols rowIndex.length < 1082 * mat.rows for x in rowIndex: x < mat.rows && x >= 0 etc. 1083 * 1084 * @param mat 1085 * @param rows 1086 * @param rowIndex 1087 * @return the input matrix 1088 */ 1089 public static Matrix plusEqualsRow(Matrix mat, Matrix rows, int[] rowIndex) { 1090 final double[][] matdata = mat.getArray(); 1091 final double[][] rowdata = rows.getArray(); 1092 int i = 0; 1093 for (final int row : rowIndex) { 1094 for (int j = 0; j < rowdata[i].length; j++) { 1095 matdata[row][j] += rowdata[i][j]; 1096 } 1097 i++; 1098 } 1099 return mat; 1100 } 1101 1102 /** 1103 * @param x 1104 * @param val 1105 * @return a new matrix for x < val 1106 */ 1107 public static Matrix lessThan(Matrix x, double val) { 1108 final Matrix retMat = x.copy(); 1109 final double[][] data = x.getArray(); 1110 final double[][] retdata = retMat.getArray(); 1111 for (int i = 0; i < data.length; i++) { 1112 for (int j = 0; j < data[i].length; j++) { 1113 retdata[i][j] = data[i][j] < val ? 1 : 0; 1114 } 1115 } 1116 return retMat; 1117 } 1118 1119 /** 1120 * @param x 1121 * @param val 1122 * @return a new matrix for x > val 1123 */ 1124 public static Matrix greaterThan(Matrix x, double val) { 1125 final Matrix retMat = x.copy(); 1126 final double[][] data = x.getArray(); 1127 final double[][] retdata = retMat.getArray(); 1128 for (int i = 0; i < data.length; i++) { 1129 for (int j = 0; j < data[i].length; j++) { 1130 retdata[i][j] = data[i][j] > val ? 1 : 0; 1131 } 1132 } 1133 return retMat; 1134 } 1135 1136 /** 1137 * @param x 1138 * @return a new matrix for x1 && x2 && ... && xn where 1139 * && means "!=0" 1140 */ 1141 public static Matrix and(Matrix... x) { 1142 final Matrix retMat = MatrixUtils.ones(x[0].getRowDimension(), x[0].getColumnDimension()); 1143 final double[][] retdata = retMat.getArray(); 1144 1145 for (final Matrix matrix : x) { 1146 final double[][] data = matrix.getArray(); 1147 for (int i = 0; i < data.length; i++) { 1148 for (int j = 0; j < data[i].length; j++) { 1149 retdata[i][j] = (retdata[i][j] != 0 && data[i][j] != 0) ? 1 : 0; 1150 } 1151 } 1152 } 1153 return retMat; 1154 } 1155 1156 /** 1157 * @param rowDimension 1158 * @param columnDimension 1159 * @return matrix of dimensions filled with ones 1160 */ 1161 public static Matrix ones(int rowDimension, int columnDimension) { 1162 final Matrix ret = new Matrix(rowDimension, columnDimension); 1163 return plus(ret, 1); 1164 } 1165 1166 /** 1167 * @param x 1168 * @return logical-and each column of x 1169 */ 1170 public static Matrix all(Matrix x) { 1171 final int cols = x.getColumnDimension(); 1172 final int rows = x.getRowDimension(); 1173 final Matrix ret = new Matrix(1, cols); 1174 final double[][] retdata = ret.getArray(); 1175 final double[][] data = x.getArray(); 1176 for (int i = 0; i < cols; i++) { 1177 boolean cool = true; 1178 for (int j = 0; j < rows; j++) { 1179 cool = data[j][i] != 0 && cool; 1180 if (!cool) 1181 break; 1182 } 1183 retdata[0][i] = cool ? 1 : 0; 1184 } 1185 return ret; 1186 } 1187 1188 /** 1189 * @param vals 1190 * @return given vals, return the array indexes where vals != 0 1191 */ 1192 public static int[] valsToIndex(double[] vals) { 1193 int nindex = 0; 1194 for (final double d : vals) { 1195 nindex += d != 0 ? 1 : 0; 1196 } 1197 1198 final int[] indexes = new int[nindex]; 1199 nindex = 0; 1200 int i = 0; 1201 for (final double d : vals) { 1202 if (d != 0) { 1203 indexes[i] = nindex; 1204 i++; 1205 } 1206 nindex++; 1207 } 1208 return indexes; 1209 } 1210 1211 /** 1212 * for every value in x greater than val set toset 1213 * 1214 * @param x 1215 * @param val 1216 * @param toset 1217 * @return same matrix handed in 1218 */ 1219 public static Matrix greaterThanSet(Matrix x, int val, int toset) { 1220 final double[][] data = x.getArray(); 1221 for (int i = 0; i < data.length; i++) { 1222 for (int j = 0; j < data[i].length; j++) { 1223 data[i][j] = data[i][j] > val ? toset : data[i][j]; 1224 } 1225 } 1226 return x; 1227 } 1228 1229 /** 1230 * for every value in x less than val set toset 1231 * 1232 * @param x 1233 * @param val 1234 * @param toset 1235 * @return same matrix handed in 1236 */ 1237 public static Matrix lessThanSet(Matrix x, int val, int toset) { 1238 final double[][] data = x.getArray(); 1239 for (int i = 0; i < data.length; i++) { 1240 for (int j = 0; j < data[i].length; j++) { 1241 data[i][j] = data[i][j] < val ? toset : data[i][j]; 1242 } 1243 } 1244 return x; 1245 } 1246 1247 /** 1248 * Subtract the given row vector from every row of the given matrix, returning 1249 * the result in a new matrix. 1250 * 1251 * @param in 1252 * the matrix 1253 * @param row 1254 * the row vector 1255 * @return the resultant matrix 1256 */ 1257 public static Matrix minusRow(Matrix in, double[] row) { 1258 final Matrix out = in.copy(); 1259 final double[][] outData = out.getArray(); 1260 final int rows = out.getRowDimension(); 1261 final int cols = out.getColumnDimension(); 1262 1263 for (int r = 0; r < rows; r++) { 1264 for (int c = 0; c < cols; c++) { 1265 outData[r][c] -= row[c]; 1266 } 1267 } 1268 1269 return out; 1270 } 1271 1272 /** 1273 * Subtract the given col vector (held as a Matrix) from every col of the given 1274 * matrix, returning the result in a new matrix. 1275 * 1276 * @param in 1277 * the matrix 1278 * @param col 1279 * the col Matrix (Only the first column is used) 1280 * @return the resultant matrix 1281 */ 1282 public static Matrix minusCol(Matrix in, Matrix col) { 1283 final Matrix out = in.copy(); 1284 final double[][] outData = out.getArray(); 1285 final int rows = out.getRowDimension(); 1286 final int cols = out.getColumnDimension(); 1287 final double[][] colArr = col.getArray(); 1288 1289 for (int r = 0; r < rows; r++) { 1290 for (int c = 0; c < cols; c++) { 1291 outData[r][c] -= colArr[r][0]; 1292 } 1293 } 1294 1295 return out; 1296 } 1297 1298 /** 1299 * Add a matrix to another inline. 1300 * 1301 * @param result 1302 * the matrix to add to 1303 * @param add 1304 * the matrix to add 1305 * @return the result matrix 1306 */ 1307 public static Matrix plusEquals(Matrix result, Matrix add) { 1308 final int rows = result.getRowDimension(); 1309 final int cols = result.getColumnDimension(); 1310 1311 final double[][] resultData = result.getArray(); 1312 final double[][] addData = add.getArray(); 1313 1314 for (int r = 0; r < rows; r++) { 1315 for (int c = 0; c < cols; c++) { 1316 resultData[r][c] += addData[r][c]; 1317 } 1318 } 1319 1320 return result; 1321 } 1322 1323 /** 1324 * Multiply a matrix by a constant inplace, returning the matrix. 1325 * 1326 * @param m 1327 * the matrix 1328 * @param val 1329 * the value to multiply by 1330 * @return the matrix 1331 */ 1332 public static Matrix times(Matrix m, double val) { 1333 final double[][] data = m.getArray(); 1334 1335 final int rows = m.getRowDimension(); 1336 final int cols = m.getColumnDimension(); 1337 1338 for (int r = 0; r < rows; r++) 1339 for (int c = 0; c < cols; c++) 1340 data[r][c] *= val; 1341 1342 return m; 1343 } 1344 1345 /** 1346 * Convert an mtj matrix into a 2d double array 1347 * 1348 * @param mat 1349 * @return a double array 1350 */ 1351 public static double[][] mtjToDoubleArray(no.uib.cipr.matrix.DenseMatrix mat) { 1352 final double[][] out = new double[mat.numRows()][mat.numColumns()]; 1353 final double[] data = mat.getData(); 1354 for (int r = 0; r < out.length; r++) { 1355 final double[] outr = out[r]; 1356 for (int c = 0; c < out[0].length; c++) { 1357 outr[c] = data[r + c * out.length]; 1358 } 1359 } 1360 return out; 1361 } 1362 1363 /** 1364 * Compute the sum of values in all rows 1365 * 1366 * @param m 1367 * the matrix 1368 * @return the sum of values across all cols in all rows 1369 */ 1370 public static Matrix sumRows(Matrix m) { 1371 final double[][] data = m.getArray(); 1372 final int cols = m.getColumnDimension(); 1373 final int rows = m.getRowDimension(); 1374 1375 final Matrix sum = new Matrix(rows, 1); 1376 final double[][] sumArr = sum.getArray(); 1377 for (int c = 0; c < cols; c++) { 1378 for (int r = 0; r < rows; r++) { 1379 sumArr[r][0] += data[r][c]; 1380 } 1381 } 1382 1383 return sum; 1384 } 1385 1386 /** 1387 * Compute the sum of values in all cols 1388 * 1389 * @param m 1390 * the matrix 1391 * @return the sum of values across all rows in all cols 1392 */ 1393 public static Matrix sumCols(Matrix m) { 1394 final double[][] data = m.getArray(); 1395 final int cols = m.getColumnDimension(); 1396 final int rows = m.getRowDimension(); 1397 1398 final Matrix sum = new Matrix(1, cols); 1399 final double[][] sumArr = sum.getArray(); 1400 for (int c = 0; c < cols; c++) { 1401 for (int r = 0; r < rows; r++) { 1402 sumArr[0][c] += data[r][c]; 1403 } 1404 } 1405 1406 return sum; 1407 } 1408 1409 /** 1410 * Generate a matrix with Gaussian distributed randoms 1411 * 1412 * @param rows 1413 * the number of rows 1414 * @param cols 1415 * the number of columns 1416 * @return a matrix containing values drawn from a 0 mean 1.0 sdev gaussian 1417 */ 1418 public static Matrix randGaussian(int rows, int cols) { 1419 final Matrix m = new Matrix(rows, cols); 1420 final double[][] d = m.getArray(); 1421 final Random r = new Random(); 1422 for (int row = 0; row < d.length; row++) { 1423 for (int col = 0; col < d[row].length; col++) { 1424 d[row][col] = r.nextGaussian(); 1425 } 1426 } 1427 return m; 1428 } 1429 1430 /** 1431 * Compute the sparsity (i.e. ratio of non-zero elements to matrix size) of the 1432 * given matrix 1433 * 1434 * @param matrix 1435 * the matrix 1436 * @return the sparsity 1437 */ 1438 public static double sparsity(SparseMatrix matrix) { 1439 final double density = matrix.used() / ((double) matrix.rowCount() * (double) matrix.columnCount()); 1440 return 1 - density; 1441 } 1442 1443 /** 1444 * Extract the diagonal component from the given matrix 1445 * 1446 * @param cv 1447 * the matrix 1448 * @return a new matrix with the diagonal values from the input matrix and other 1449 * values set to zero. 1450 */ 1451 public static Matrix diag(Matrix cv) { 1452 final Matrix d = new Matrix(cv.getRowDimension(), cv.getColumnDimension()); 1453 1454 for (int i = 0; i < Math.min(cv.getRowDimension(), cv.getColumnDimension()); i++) 1455 d.set(i, i, cv.get(i, i)); 1456 1457 return d; 1458 } 1459 1460 /** 1461 * Extract the diagonal component from the given matrix 1462 * 1463 * @param cv 1464 * the matrix 1465 * @return a new matrix with the diagonal values from the input matrix and other 1466 * values set to zero. 1467 */ 1468 public static double[] diagVector(Matrix cv) { 1469 final double[] d = new double[Math.min(cv.getRowDimension(), cv.getColumnDimension())]; 1470 1471 for (int i = 0; i < Math.min(cv.getRowDimension(), cv.getColumnDimension()); i++) 1472 d[i] = cv.get(i, i); 1473 1474 return d; 1475 } 1476 1477 /** 1478 * Format a matrix as a single-line string suitable for using in matlab or 1479 * octave 1480 * 1481 * @param mat 1482 * the matrix to format 1483 * @return the string 1484 */ 1485 public static String toMatlabString(Matrix mat) { 1486 return "[" + toString(mat).trim().replace("\n ", ";").replace(" ", ",") + "]"; 1487 } 1488 1489 /** 1490 * Format a matrix as a single-line string suitable for using in python 1491 * 1492 * @param mat 1493 * the matrix to format 1494 * @return the string 1495 */ 1496 public static String toPythonString(Matrix mat) { 1497 return "[[" + toString(mat).trim().replace("\n ", "][").replace(" ", ",") + "]]"; 1498 } 1499 1500 /** 1501 * @param mat 1502 * @return trace of the matrix 1503 */ 1504 public static double trace(Matrix mat) { 1505 double sum = 0; 1506 for (int i = 0; i < mat.getRowDimension(); i++) { 1507 sum += mat.get(i, i); 1508 } 1509 return sum; 1510 } 1511 1512 /** 1513 * Solves the system <code>Ax = 0</code>, returning the vector x as an array. 1514 * Internally computes the least-squares solution using the SVD of 1515 * <code>A</code>. 1516 * 1517 * @param A 1518 * the matrix describing the system 1519 * @return the solution vector 1520 */ 1521 public static double[] solveHomogeneousSystem(Matrix A) { 1522 return solveHomogeneousSystem(new DenseMatrix(A.getArray())); 1523 } 1524 1525 /** 1526 * Solves the system <code>Ax = 0</code>, returning the vector x as an array. 1527 * Internally computes the least-squares solution using the SVD of 1528 * <code>A</code>. 1529 * 1530 * @param A 1531 * the matrix describing the system 1532 * @return the solution vector 1533 */ 1534 public static double[] solveHomogeneousSystem(double[][] A) { 1535 return solveHomogeneousSystem(new DenseMatrix(A)); 1536 } 1537 1538 /** 1539 * Solves the system <code>Ax = 0</code>, returning the vector x as an array. 1540 * Internally computes the least-squares solution using the SVD of 1541 * <code>A</code>. 1542 * 1543 * @param A 1544 * the matrix describing the system 1545 * @return the solution vector 1546 */ 1547 public static double[] solveHomogeneousSystem(DenseMatrix A) { 1548 try { 1549 final SVD svd = SVD.factorize(A); 1550 1551 final double[] x = new double[svd.getVt().numRows()]; 1552 final int c = svd.getVt().numColumns() - 1; 1553 1554 for (int i = 0; i < x.length; i++) { 1555 x[i] = svd.getVt().get(c, i); 1556 } 1557 1558 return x; 1559 } catch (final NotConvergedException e) { 1560 throw new RuntimeException(e); 1561 } 1562 } 1563 1564 /** 1565 * Format a matrix as a single-line string suitable for using in java 1566 * 1567 * @param mat 1568 * the matrix to format 1569 * @return the string 1570 */ 1571 public static String toJavaString(Matrix mat) { 1572 return "{" + toString(mat).trim().replaceAll("^", "{").replace("\n ", "},{").replace(" ", ",") + "}}"; 1573 } 1574 1575 /** 1576 * Construct a matrix from a row-packed (i.e. row-by-row) vector of the data. 1577 * 1578 * @param vector 1579 * the row-packed vector 1580 * @param ncols 1581 * the number of columns 1582 * @return the reconstructed matrix 1583 */ 1584 public static Matrix fromRowPacked(double[] vector, int ncols) { 1585 final int nrows = vector.length / ncols; 1586 1587 final Matrix a = new Matrix(nrows, ncols); 1588 final double[][] ad = a.getArray(); 1589 for (int r = 0, i = 0; r < nrows; r++) 1590 for (int c = 0; c < ncols; c++, i++) 1591 ad[r][c] = vector[i]; 1592 1593 return a; 1594 } 1595 1596 /** 1597 * Compute the covariance matrix of the given samples (assumed each sample is a 1598 * row). 1599 * 1600 * @param m 1601 * the samples matrix 1602 * @return the covariance matrix 1603 */ 1604 public static Matrix covariance(Matrix m) { 1605 final int N = m.getRowDimension(); 1606 return times(m.transpose().times(m), 1.0 / (N > 1 ? N - 1 : N)); 1607 } 1608 1609 /** 1610 * For each element of X, sign(X) returns 1 if the element is greater than zero, 1611 * 0 if it equals zero and -1 if it is less than zero. 1612 * 1613 * @param m 1614 * the matrix 1615 * @return the sign matrix 1616 */ 1617 public static Matrix sign(Matrix m) { 1618 final Matrix o = new Matrix(m.getRowDimension(), m.getColumnDimension()); 1619 final double[][] md = m.getArray(); 1620 final double[][] od = o.getArray(); 1621 1622 for (int r = 0; r < o.getRowDimension(); r++) { 1623 for (int c = 0; c < o.getColumnDimension(); c++) { 1624 if (md[r][c] > 0) 1625 od[r][c] = 1; 1626 if (md[r][c] < 0) 1627 od[r][c] = -1; 1628 } 1629 } 1630 1631 return o; 1632 } 1633 1634 /** 1635 * Return a copy of the input matrix where every value is the exponential of the 1636 * elements, e to the X. 1637 * 1638 * @param m 1639 * the input matrix 1640 * @return the exponential matrix 1641 */ 1642 public static Matrix exp(Matrix m) { 1643 final Matrix o = new Matrix(m.getRowDimension(), m.getColumnDimension()); 1644 final double[][] md = m.getArray(); 1645 final double[][] od = o.getArray(); 1646 1647 for (int r = 0; r < o.getRowDimension(); r++) { 1648 for (int c = 0; c < o.getColumnDimension(); c++) { 1649 od[r][c] = Math.exp(md[r][c]); 1650 } 1651 } 1652 1653 return o; 1654 } 1655 1656 /** 1657 * Return a copy of the input matrix where every value is the hyperbolic tangent 1658 * of the elements. 1659 * 1660 * @param m 1661 * the input matrix 1662 * @return the tanh matrix 1663 */ 1664 public static Matrix tanh(Matrix m) { 1665 final Matrix o = new Matrix(m.getRowDimension(), m.getColumnDimension()); 1666 final double[][] md = m.getArray(); 1667 final double[][] od = o.getArray(); 1668 1669 for (int r = 0; r < o.getRowDimension(); r++) { 1670 for (int c = 0; c < o.getColumnDimension(); c++) { 1671 od[r][c] = Math.tanh(md[r][c]); 1672 } 1673 } 1674 1675 return o; 1676 } 1677}