Java tutorial
/* * This software is a cooperative product of The MathWorks and the National * Institute of Standards and Technology (NIST) which has been released to the * public domain. Neither The MathWorks nor NIST assumes any responsibility * whatsoever for its use by other parties, and makes no guarantees, expressed * or implied, about its quality, reliability, or any other characteristic. */ /* * Matrix.java * Copyright (C) 1999 The Mathworks and NIST and 2005 University of Waikato, * Hamilton, New Zealand * */ package weka.core.matrix; import java.io.BufferedReader; import java.io.LineNumberReader; import java.io.PrintWriter; import java.io.Reader; import java.io.Serializable; import java.io.StreamTokenizer; import java.io.StringReader; import java.io.StringWriter; import java.io.Writer; import java.text.DecimalFormat; import java.text.DecimalFormatSymbols; import java.text.NumberFormat; import java.util.Locale; import java.util.StringTokenizer; import weka.core.RevisionHandler; import weka.core.RevisionUtils; import weka.core.Utils; /** * Jama = Java Matrix class. * <P> * The Java Matrix Class provides the fundamental operations of numerical linear * algebra. Various constructors create Matrices from two dimensional arrays of * double precision floating point numbers. Various "gets" and "sets" provide * access to submatrices and matrix elements. Several methods implement basic * matrix arithmetic, including matrix addition and multiplication, matrix * norms, and element-by-element array operations. Methods for reading and * printing matrices are also included. All the operations in this version of * the Matrix Class involve real matrices. Complex matrices may be handled in a * future version. * <P> * Five fundamental matrix decompositions, which consist of pairs or triples of * matrices, permutation vectors, and the like, produce results in five * decomposition classes. These decompositions are accessed by the Matrix class * to compute solutions of simultaneous linear equations, determinants, inverses * and other matrix functions. The five decompositions are: * <P> * <UL> * <LI>Cholesky Decomposition of symmetric, positive definite matrices. * <LI>LU Decomposition of rectangular matrices. * <LI>QR Decomposition of rectangular matrices. * <LI>Singular Value Decomposition of rectangular matrices. * <LI>Eigenvalue Decomposition of both symmetric and nonsymmetric square * matrices. * </UL> * <DL> * <DT><B>Example of use:</B></DT> * <P> * <DD>Solve a linear system A x = b and compute the residual norm, ||b - A x||. * <P> * * <PRE> * double[][] vals = { { 1., 2., 3 }, { 4., 5., 6. }, { 7., 8., 10. } }; * Matrix A = new Matrix(vals); * Matrix b = Matrix.random(3, 1); * Matrix x = A.solve(b); * Matrix r = A.times(x).minus(b); * double rnorm = r.normInf(); * </PRE> * * </DD> * </DL> * <p/> * Adapted from the <a href="http://math.nist.gov/javanumerics/jama/" * target="_blank">JAMA</a> package. Additional methods are tagged with the * <code>@author</code> tag. * * @author The Mathworks and NIST * @author Fracpete (fracpete at waikato dot ac dot nz) * @version $Revision$ */ public class Matrix implements Cloneable, Serializable, RevisionHandler { /** for serialization */ private static final long serialVersionUID = 7856794138418366180L; /** * Array for internal storage of elements. * * @serial internal array storage. */ protected double[][] A; /** * Row and column dimensions. * * @serial row dimension. * @serial column dimension. */ protected int m, n; /** * Construct an m-by-n matrix of zeros. * * @param m Number of rows. * @param n Number of colums. */ public Matrix(int m, int n) { this.m = m; this.n = n; A = new double[m][n]; } /** * Construct an m-by-n constant matrix. * * @param m Number of rows. * @param n Number of colums. * @param s Fill the matrix with this scalar value. */ public Matrix(int m, int n, double s) { this.m = m; this.n = n; A = new double[m][n]; for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { A[i][j] = s; } } } /** * Construct a matrix from a 2-D array. * * @param A Two-dimensional array of doubles. * @throws IllegalArgumentException All rows must have the same length * @see #constructWithCopy */ public Matrix(double[][] A) { m = A.length; n = A[0].length; for (int i = 0; i < m; i++) { if (A[i].length != n) { throw new IllegalArgumentException("All rows must have the same length."); } } this.A = A; } /** * Construct a matrix quickly without checking arguments. * * @param A Two-dimensional array of doubles. * @param m Number of rows. * @param n Number of colums. */ public Matrix(double[][] A, int m, int n) { this.A = A; this.m = m; this.n = n; } /** * Construct a matrix from a one-dimensional packed array * * @param vals One-dimensional array of doubles, packed by columns (ala * Fortran). * @param m Number of rows. * @throws IllegalArgumentException Array length must be a multiple of m. */ public Matrix(double vals[], int m) { this.m = m; n = (m != 0 ? vals.length / m : 0); if (m * n != vals.length) { throw new IllegalArgumentException("Array length must be a multiple of m."); } A = new double[m][n]; for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { A[i][j] = vals[i + j * m]; } } } /** * Reads a matrix from a reader. The first line in the file should contain the * number of rows and columns. Subsequent lines contain elements of the * matrix. (FracPete: taken from old weka.core.Matrix class) * * @param r the reader containing the matrix * @throws Exception if an error occurs * @see #write(Writer) */ public Matrix(Reader r) throws Exception { LineNumberReader lnr = new LineNumberReader(r); String line; int currentRow = -1; while ((line = lnr.readLine()) != null) { // Comments if (line.startsWith("%")) { continue; } StringTokenizer st = new StringTokenizer(line); // Ignore blank lines if (!st.hasMoreTokens()) { continue; } if (currentRow < 0) { int rows = Integer.parseInt(st.nextToken()); if (!st.hasMoreTokens()) { throw new Exception("Line " + lnr.getLineNumber() + ": expected number of columns"); } int cols = Integer.parseInt(st.nextToken()); A = new double[rows][cols]; m = rows; n = cols; currentRow++; continue; } else { if (currentRow == getRowDimension()) { throw new Exception("Line " + lnr.getLineNumber() + ": too many rows provided"); } for (int i = 0; i < getColumnDimension(); i++) { if (!st.hasMoreTokens()) { throw new Exception("Line " + lnr.getLineNumber() + ": too few matrix elements provided"); } set(currentRow, i, Double.valueOf(st.nextToken()).doubleValue()); } currentRow++; } } if (currentRow == -1) { throw new Exception("Line " + lnr.getLineNumber() + ": expected number of rows"); } else if (currentRow != getRowDimension()) { throw new Exception("Line " + lnr.getLineNumber() + ": too few rows provided"); } } /** * Construct a matrix from a copy of a 2-D array. * * @param A Two-dimensional array of doubles. * @throws IllegalArgumentException All rows must have the same length */ public static Matrix constructWithCopy(double[][] A) { int m = A.length; int n = A[0].length; Matrix X = new Matrix(m, n); double[][] C = X.getArray(); for (int i = 0; i < m; i++) { if (A[i].length != n) { throw new IllegalArgumentException("All rows must have the same length."); } for (int j = 0; j < n; j++) { C[i][j] = A[i][j]; } } return X; } /** * Make a deep copy of a matrix */ public Matrix copy() { Matrix X = new Matrix(m, n); double[][] C = X.getArray(); for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { C[i][j] = A[i][j]; } } return X; } /** * Clone the Matrix object. */ @Override public Object clone() { return this.copy(); } /** * Access the internal two-dimensional array. * * @return Pointer to the two-dimensional array of matrix elements. */ public double[][] getArray() { return A; } /** * Copy the internal two-dimensional array. * * @return Two-dimensional array copy of matrix elements. */ public double[][] getArrayCopy() { double[][] C = new double[m][n]; for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { C[i][j] = A[i][j]; } } return C; } /** * Make a one-dimensional column packed copy of the internal array. * * @return Matrix elements packed in a one-dimensional array by columns. */ public double[] getColumnPackedCopy() { double[] vals = new double[m * n]; for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { vals[i + j * m] = A[i][j]; } } return vals; } /** * Make a one-dimensional row packed copy of the internal array. * * @return Matrix elements packed in a one-dimensional array by rows. */ public double[] getRowPackedCopy() { double[] vals = new double[m * n]; for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { vals[i * n + j] = A[i][j]; } } return vals; } /** * Get row dimension. * * @return m, the number of rows. */ public int getRowDimension() { return m; } /** * Get column dimension. * * @return n, the number of columns. */ public int getColumnDimension() { return n; } /** * Get a single element. * * @param i Row index. * @param j Column index. * @return A(i,j) * @throws ArrayIndexOutOfBoundsException */ public double get(int i, int j) { return A[i][j]; } /** * Get a submatrix. * * @param i0 Initial row index * @param i1 Final row index * @param j0 Initial column index * @param j1 Final column index * @return A(i0:i1,j0:j1) * @throws ArrayIndexOutOfBoundsException Submatrix indices */ public Matrix getMatrix(int i0, int i1, int j0, int j1) { Matrix X = new Matrix(i1 - i0 + 1, j1 - j0 + 1); double[][] B = X.getArray(); try { for (int i = i0; i <= i1; i++) { for (int j = j0; j <= j1; j++) { B[i - i0][j - j0] = A[i][j]; } } } catch (ArrayIndexOutOfBoundsException e) { throw new ArrayIndexOutOfBoundsException("Submatrix indices"); } return X; } /** * Get a submatrix. * * @param r Array of row indices. * @param c Array of column indices. * @return A(r(:),c(:)) * @throws ArrayIndexOutOfBoundsException Submatrix indices */ public Matrix getMatrix(int[] r, int[] c) { Matrix X = new Matrix(r.length, c.length); double[][] B = X.getArray(); try { for (int i = 0; i < r.length; i++) { for (int j = 0; j < c.length; j++) { B[i][j] = A[r[i]][c[j]]; } } } catch (ArrayIndexOutOfBoundsException e) { throw new ArrayIndexOutOfBoundsException("Submatrix indices"); } return X; } /** * Get a submatrix. * * @param i0 Initial row index * @param i1 Final row index * @param c Array of column indices. * @return A(i0:i1,c(:)) * @throws ArrayIndexOutOfBoundsException Submatrix indices */ public Matrix getMatrix(int i0, int i1, int[] c) { Matrix X = new Matrix(i1 - i0 + 1, c.length); double[][] B = X.getArray(); try { for (int i = i0; i <= i1; i++) { for (int j = 0; j < c.length; j++) { B[i - i0][j] = A[i][c[j]]; } } } catch (ArrayIndexOutOfBoundsException e) { throw new ArrayIndexOutOfBoundsException("Submatrix indices"); } return X; } /** * Get a submatrix. * * @param r Array of row indices. * @param j0 Initial column index * @param j1 Final column index * @return A(r(:),j0:j1) * @throws ArrayIndexOutOfBoundsException Submatrix indices */ public Matrix getMatrix(int[] r, int j0, int j1) { Matrix X = new Matrix(r.length, j1 - j0 + 1); double[][] B = X.getArray(); try { for (int i = 0; i < r.length; i++) { for (int j = j0; j <= j1; j++) { B[i][j - j0] = A[r[i]][j]; } } } catch (ArrayIndexOutOfBoundsException e) { throw new ArrayIndexOutOfBoundsException("Submatrix indices"); } return X; } /** * Set a single element. * * @param i Row index. * @param j Column index. * @param s A(i,j). * @throws ArrayIndexOutOfBoundsException */ public void set(int i, int j, double s) { A[i][j] = s; } /** * Set a submatrix. * * @param i0 Initial row index * @param i1 Final row index * @param j0 Initial column index * @param j1 Final column index * @param X A(i0:i1,j0:j1) * @throws ArrayIndexOutOfBoundsException Submatrix indices */ public void setMatrix(int i0, int i1, int j0, int j1, Matrix X) { try { for (int i = i0; i <= i1; i++) { for (int j = j0; j <= j1; j++) { A[i][j] = X.get(i - i0, j - j0); } } } catch (ArrayIndexOutOfBoundsException e) { throw new ArrayIndexOutOfBoundsException("Submatrix indices"); } } /** * Set a submatrix. * * @param r Array of row indices. * @param c Array of column indices. * @param X A(r(:),c(:)) * @throws ArrayIndexOutOfBoundsException Submatrix indices */ public void setMatrix(int[] r, int[] c, Matrix X) { try { for (int i = 0; i < r.length; i++) { for (int j = 0; j < c.length; j++) { A[r[i]][c[j]] = X.get(i, j); } } } catch (ArrayIndexOutOfBoundsException e) { throw new ArrayIndexOutOfBoundsException("Submatrix indices"); } } /** * Set a submatrix. * * @param r Array of row indices. * @param j0 Initial column index * @param j1 Final column index * @param X A(r(:),j0:j1) * @throws ArrayIndexOutOfBoundsException Submatrix indices */ public void setMatrix(int[] r, int j0, int j1, Matrix X) { try { for (int i = 0; i < r.length; i++) { for (int j = j0; j <= j1; j++) { A[r[i]][j] = X.get(i, j - j0); } } } catch (ArrayIndexOutOfBoundsException e) { throw new ArrayIndexOutOfBoundsException("Submatrix indices"); } } /** * Set a submatrix. * * @param i0 Initial row index * @param i1 Final row index * @param c Array of column indices. * @param X A(i0:i1,c(:)) * @throws ArrayIndexOutOfBoundsException Submatrix indices */ public void setMatrix(int i0, int i1, int[] c, Matrix X) { try { for (int i = i0; i <= i1; i++) { for (int j = 0; j < c.length; j++) { A[i][c[j]] = X.get(i - i0, j); } } } catch (ArrayIndexOutOfBoundsException e) { throw new ArrayIndexOutOfBoundsException("Submatrix indices"); } } /** * Returns true if the matrix is symmetric. (FracPete: taken from old * weka.core.Matrix class) * * @return boolean true if matrix is symmetric. */ public boolean isSymmetric() { int nr = A.length, nc = A[0].length; if (nr != nc) { return false; } for (int i = 0; i < nc; i++) { for (int j = 0; j < i; j++) { if (A[i][j] != A[j][i]) { return false; } } } return true; } /** * returns whether the matrix is a square matrix or not. * * @return true if the matrix is a square matrix */ public boolean isSquare() { return (getRowDimension() == getColumnDimension()); } /** * Matrix transpose. * * @return A' */ public Matrix transpose() { Matrix X = new Matrix(n, m); double[][] C = X.getArray(); for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { C[j][i] = A[i][j]; } } return X; } /** * One norm * * @return maximum column sum. */ public double norm1() { double f = 0; for (int j = 0; j < n; j++) { double s = 0; for (int i = 0; i < m; i++) { s += Math.abs(A[i][j]); } f = Math.max(f, s); } return f; } /** * Two norm * * @return maximum singular value. */ public double norm2() { return (new SingularValueDecomposition(this).norm2()); } /** * Infinity norm * * @return maximum row sum. */ public double normInf() { double f = 0; for (int i = 0; i < m; i++) { double s = 0; for (int j = 0; j < n; j++) { s += Math.abs(A[i][j]); } f = Math.max(f, s); } return f; } /** * Frobenius norm * * @return sqrt of sum of squares of all elements. */ public double normF() { double f = 0; for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { f = Maths.hypot(f, A[i][j]); } } return f; } /** * Unary minus * * @return -A */ public Matrix uminus() { Matrix X = new Matrix(m, n); double[][] C = X.getArray(); for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { C[i][j] = -A[i][j]; } } return X; } /** * C = A + B * * @param B another matrix * @return A + B */ public Matrix plus(Matrix B) { checkMatrixDimensions(B); Matrix X = new Matrix(m, n); double[][] C = X.getArray(); for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { C[i][j] = A[i][j] + B.A[i][j]; } } return X; } /** * A = A + B * * @param B another matrix * @return A + B */ public Matrix plusEquals(Matrix B) { checkMatrixDimensions(B); for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { A[i][j] = A[i][j] + B.A[i][j]; } } return this; } /** * C = A - B * * @param B another matrix * @return A - B */ public Matrix minus(Matrix B) { checkMatrixDimensions(B); Matrix X = new Matrix(m, n); double[][] C = X.getArray(); for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { C[i][j] = A[i][j] - B.A[i][j]; } } return X; } /** * A = A - B * * @param B another matrix * @return A - B */ public Matrix minusEquals(Matrix B) { checkMatrixDimensions(B); for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { A[i][j] = A[i][j] - B.A[i][j]; } } return this; } /** * Element-by-element multiplication, C = A.*B * * @param B another matrix * @return A.*B */ public Matrix arrayTimes(Matrix B) { checkMatrixDimensions(B); Matrix X = new Matrix(m, n); double[][] C = X.getArray(); for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { C[i][j] = A[i][j] * B.A[i][j]; } } return X; } /** * Element-by-element multiplication in place, A = A.*B * * @param B another matrix * @return A.*B */ public Matrix arrayTimesEquals(Matrix B) { checkMatrixDimensions(B); for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { A[i][j] = A[i][j] * B.A[i][j]; } } return this; } /** * Element-by-element right division, C = A./B * * @param B another matrix * @return A./B */ public Matrix arrayRightDivide(Matrix B) { checkMatrixDimensions(B); Matrix X = new Matrix(m, n); double[][] C = X.getArray(); for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { C[i][j] = A[i][j] / B.A[i][j]; } } return X; } /** * Element-by-element right division in place, A = A./B * * @param B another matrix * @return A./B */ public Matrix arrayRightDivideEquals(Matrix B) { checkMatrixDimensions(B); for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { A[i][j] = A[i][j] / B.A[i][j]; } } return this; } /** * Element-by-element left division, C = A.\B * * @param B another matrix * @return A.\B */ public Matrix arrayLeftDivide(Matrix B) { checkMatrixDimensions(B); Matrix X = new Matrix(m, n); double[][] C = X.getArray(); for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { C[i][j] = B.A[i][j] / A[i][j]; } } return X; } /** * Element-by-element left division in place, A = A.\B * * @param B another matrix * @return A.\B */ public Matrix arrayLeftDivideEquals(Matrix B) { checkMatrixDimensions(B); for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { A[i][j] = B.A[i][j] / A[i][j]; } } return this; } /** * Multiply a matrix by a scalar, C = s*A * * @param s scalar * @return s*A */ public Matrix times(double s) { Matrix X = new Matrix(m, n); double[][] C = X.getArray(); for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { C[i][j] = s * A[i][j]; } } return X; } /** * Multiply a matrix by a scalar in place, A = s*A * * @param s scalar * @return replace A by s*A */ public Matrix timesEquals(double s) { for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { A[i][j] = s * A[i][j]; } } return this; } /** * Linear algebraic matrix multiplication, A * B * * @param B another matrix * @return Matrix product, A * B * @throws IllegalArgumentException Matrix inner dimensions must agree. */ public Matrix times(Matrix B) { if (B.m != n) { throw new IllegalArgumentException("Matrix inner dimensions must agree."); } Matrix X = new Matrix(m, B.n); double[][] C = X.getArray(); double[] Bcolj = new double[n]; for (int j = 0; j < B.n; j++) { for (int k = 0; k < n; k++) { Bcolj[k] = B.A[k][j]; } for (int i = 0; i < m; i++) { double[] Arowi = A[i]; double s = 0; for (int k = 0; k < n; k++) { s += Arowi[k] * Bcolj[k]; } C[i][j] = s; } } return X; } /** * LU Decomposition * * @return LUDecomposition * @see LUDecomposition */ public LUDecomposition lu() { return new LUDecomposition(this); } /** * QR Decomposition * * @return QRDecomposition * @see QRDecomposition */ public QRDecomposition qr() { return new QRDecomposition(this); } /** * Cholesky Decomposition * * @return CholeskyDecomposition * @see CholeskyDecomposition */ public CholeskyDecomposition chol() { return new CholeskyDecomposition(this); } /** * Singular Value Decomposition * * @return SingularValueDecomposition * @see SingularValueDecomposition */ public SingularValueDecomposition svd() { return new SingularValueDecomposition(this); } /** * Eigenvalue Decomposition * * @return EigenvalueDecomposition * @see EigenvalueDecomposition */ public EigenvalueDecomposition eig() { return new EigenvalueDecomposition(this); } /** * Solve A*X = B * * @param B right hand side * @return solution if A is square, least squares solution otherwise */ public Matrix solve(Matrix B) { return (m == n ? (new LUDecomposition(this)).solve(B) : (new QRDecomposition(this)).solve(B)); } /** * Solve X*A = B, which is also A'*X' = B' * * @param B right hand side * @return solution if A is square, least squares solution otherwise. */ public Matrix solveTranspose(Matrix B) { return transpose().solve(B.transpose()); } /** * Matrix inverse or pseudoinverse * * @return inverse(A) if A is square, pseudoinverse otherwise. */ public Matrix inverse() { return solve(identity(m, m)); } /** * returns the square root of the matrix, i.e., X from the equation X*X = A.<br/> * Steps in the Calculation (see <a href= * "http://www.mathworks.com/access/helpdesk/help/techdoc/ref/sqrtm.html" * target="blank"><code>sqrtm</code></a> in Matlab):<br/> * <ol> * <li>perform eigenvalue decomposition<br/> * [V,D]=eig(A)</li> * <li>take the square root of all elements in D (only the ones with positive * sign are considered for further computation)<br/> * S=sqrt(D)</li> * <li>calculate the root<br/> * X=V*S/V, which can be also written as X=(V'\(V*S)')'</li> * </ol> * <p/> * <b>Note:</b> since this method uses other high-level methods, it generates * several instances of matrices. This can be problematic with large matrices. * <p/> * Examples: * <ol> * <li> * * <pre> * X = * 5 -4 1 0 0 * -4 6 -4 1 0 * 1 -4 6 -4 1 * 0 1 -4 6 -4 * 0 0 1 -4 5 * * sqrt(X) = * 2 -1 -0 -0 -0 * -1 2 -1 0 -0 * 0 -1 2 -1 0 * -0 0 -1 2 -1 * -0 -0 -0 -1 2 * * Matrix m = new Matrix(new double[][]{{5,-4,1,0,0},{-4,6,-4,1,0},{1,-4,6,-4,1},{0,1,-4,6,-4},{0,0,1,-4,5}}); * </pre> * * </li> * <li> * * <pre> * X = * 7 10 * 15 22 * * sqrt(X) = * 1.5667 1.7408 * 2.6112 4.1779 * * Matrix m = new Matrix(new double[][]{{7, 10},{15, 22}}); * </pre> * * </li> * </ol> * * @return sqrt(A) */ public Matrix sqrt() { EigenvalueDecomposition evd; Matrix s; Matrix v; Matrix d; Matrix result; Matrix a; Matrix b; int i; int n; result = null; // eigenvalue decomp. // [V, D] = eig(A) with A = this evd = this.eig(); v = evd.getV(); d = evd.getD(); // S = sqrt of cells of D s = new Matrix(d.getRowDimension(), d.getColumnDimension()); for (i = 0; i < s.getRowDimension(); i++) { for (n = 0; n < s.getColumnDimension(); n++) { s.set(i, n, StrictMath.sqrt(d.get(i, n))); } } // to calculate: // result = V*S/V // // with X = B/A // and B/A = (A'\B')' // and V=A and V*S=B // we get // result = (V'\(V*S)')' // // A*X = B // X = A\B // which is // X = A.solve(B) // // with A=V' and B=(V*S)' // we get // X = V'.solve((V*S)') // or // result = X' // // which is in full length // result = (V'.solve((V*S)'))' a = v.inverse(); b = v.times(s).inverse(); v = null; d = null; evd = null; s = null; result = a.solve(b).inverse(); return result; } /** * Performs a (ridged) linear regression. (FracPete: taken from old * weka.core.Matrix class) * * @param y the dependent variable vector * @param ridge the ridge parameter * @return the coefficients * @throws IllegalArgumentException if not successful */ public LinearRegression regression(Matrix y, double ridge) { return new LinearRegression(this, y, ridge); } /** * Performs a weighted (ridged) linear regression. (FracPete: taken from old * weka.core.Matrix class) * * @param y the dependent variable vector * @param w the array of data point weights * @param ridge the ridge parameter * @return the coefficients * @throws IllegalArgumentException if the wrong number of weights were * provided. */ public final LinearRegression regression(Matrix y, double[] w, double ridge) { return new LinearRegression(this, y, w, ridge); } /** * Matrix determinant * * @return determinant */ public double det() { return new LUDecomposition(this).det(); } /** * Matrix rank * * @return effective numerical rank, obtained from SVD. */ public int rank() { return new SingularValueDecomposition(this).rank(); } /** * Matrix condition (2 norm) * * @return ratio of largest to smallest singular value. */ public double cond() { return new SingularValueDecomposition(this).cond(); } /** * Matrix trace. * * @return sum of the diagonal elements. */ public double trace() { double t = 0; for (int i = 0; i < Math.min(m, n); i++) { t += A[i][i]; } return t; } /** * Generate matrix with random elements * * @param m Number of rows. * @param n Number of colums. * @return An m-by-n matrix with uniformly distributed random elements. */ public static Matrix random(int m, int n) { Matrix A = new Matrix(m, n); double[][] X = A.getArray(); for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { X[i][j] = Math.random(); } } return A; } /** * Generate identity matrix * * @param m Number of rows. * @param n Number of colums. * @return An m-by-n matrix with ones on the diagonal and zeros elsewhere. */ public static Matrix identity(int m, int n) { Matrix A = new Matrix(m, n); double[][] X = A.getArray(); for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { X[i][j] = (i == j ? 1.0 : 0.0); } } return A; } /** * Print the matrix to stdout. Line the elements up in columns with a * Fortran-like 'Fw.d' style format. * * @param w Column width. * @param d Number of digits after the decimal. */ public void print(int w, int d) { print(new PrintWriter(System.out, true), w, d); } /** * Print the matrix to the output stream. Line the elements up in columns with * a Fortran-like 'Fw.d' style format. * * @param output Output stream. * @param w Column width. * @param d Number of digits after the decimal. */ public void print(PrintWriter output, int w, int d) { DecimalFormat format = new DecimalFormat(); format.setDecimalFormatSymbols(new DecimalFormatSymbols(Locale.US)); format.setMinimumIntegerDigits(1); format.setMaximumFractionDigits(d); format.setMinimumFractionDigits(d); format.setGroupingUsed(false); print(output, format, w + 2); } /** * Print the matrix to stdout. Line the elements up in columns. Use the format * object, and right justify within columns of width characters. Note that is * the matrix is to be read back in, you probably will want to use a * NumberFormat that is set to US Locale. * * @param format A Formatting object for individual elements. * @param width Field width for each column. * @see java.text.DecimalFormat#setDecimalFormatSymbols */ public void print(NumberFormat format, int width) { print(new PrintWriter(System.out, true), format, width); } // DecimalFormat is a little disappointing coming from Fortran or C's printf. // Since it doesn't pad on the left, the elements will come out different // widths. Consequently, we'll pass the desired column width in as an // argument and do the extra padding ourselves. /** * Print the matrix to the output stream. Line the elements up in columns. Use * the format object, and right justify within columns of width characters. * Note that is the matrix is to be read back in, you probably will want to * use a NumberFormat that is set to US Locale. * * @param output the output stream. * @param format A formatting object to format the matrix elements * @param width Column width. * @see java.text.DecimalFormat#setDecimalFormatSymbols */ public void print(PrintWriter output, NumberFormat format, int width) { output.println(); // start on new line. for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { String s = format.format(A[i][j]); // format the number int padding = Math.max(1, width - s.length()); // At _least_ 1 space for (int k = 0; k < padding; k++) { output.print(' '); } output.print(s); } output.println(); } output.println(); // end with blank line. } /** * Read a matrix from a stream. The format is the same the print method, so * printed matrices can be read back in (provided they were printed using US * Locale). Elements are separated by whitespace, all the elements for each * row appear on a single line, the last row is followed by a blank line. * <p/> * Note: This format differs from the one that can be read via the * Matrix(Reader) constructor! For that format, the write(Writer) method is * used (from the original weka.core.Matrix class). * * @param input the input stream. * @see #Matrix(Reader) * @see #write(Writer) */ public static Matrix read(BufferedReader input) throws java.io.IOException { StreamTokenizer tokenizer = new StreamTokenizer(input); // Although StreamTokenizer will parse numbers, it doesn't recognize // scientific notation (E or D); however, Double.valueOf does. // The strategy here is to disable StreamTokenizer's number parsing. // We'll only get whitespace delimited words, EOL's and EOF's. // These words should all be numbers, for Double.valueOf to parse. tokenizer.resetSyntax(); tokenizer.wordChars(0, 255); tokenizer.whitespaceChars(0, ' '); tokenizer.eolIsSignificant(true); java.util.Vector<Object> v = new java.util.Vector<Object>(); // Ignore initial empty lines while (tokenizer.nextToken() == StreamTokenizer.TT_EOL) { ; } if (tokenizer.ttype == StreamTokenizer.TT_EOF) { throw new java.io.IOException("Unexpected EOF on matrix read."); } do { v.addElement(Double.valueOf(tokenizer.sval)); // Read & store 1st row. } while (tokenizer.nextToken() == StreamTokenizer.TT_WORD); int n = v.size(); // Now we've got the number of columns! double row[] = new double[n]; for (int j = 0; j < n; j++) { row[j] = ((Double) v.elementAt(j)).doubleValue(); } v.removeAllElements(); v.addElement(row); // Start storing rows instead of columns. while (tokenizer.nextToken() == StreamTokenizer.TT_WORD) { // While non-empty lines v.addElement(row = new double[n]); int j = 0; do { if (j >= n) { throw new java.io.IOException("Row " + v.size() + " is too long."); } row[j++] = Double.valueOf(tokenizer.sval).doubleValue(); } while (tokenizer.nextToken() == StreamTokenizer.TT_WORD); if (j < n) { throw new java.io.IOException("Row " + v.size() + " is too short."); } } int m = v.size(); // Now we've got the number of rows. double[][] A = new double[m][]; v.copyInto(A); // copy the rows out of the vector return new Matrix(A); } /** * Check if size(A) == size(B) */ private void checkMatrixDimensions(Matrix B) { if (B.m != m || B.n != n) { throw new IllegalArgumentException("Matrix dimensions must agree."); } } /** * Writes out a matrix. The format can be read via the Matrix(Reader) * constructor. (FracPete: taken from old weka.core.Matrix class) * * @param w the output Writer * @throws Exception if an error occurs * @see #Matrix(Reader) */ public void write(Writer w) throws Exception { w.write("% Rows\tColumns\n"); w.write("" + getRowDimension() + "\t" + getColumnDimension() + "\n"); w.write("% Matrix elements\n"); for (int i = 0; i < getRowDimension(); i++) { for (int j = 0; j < getColumnDimension(); j++) { w.write("" + get(i, j) + "\t"); } w.write("\n"); } w.flush(); } /** * Converts a matrix to a string. (FracPete: taken from old weka.core.Matrix * class) * * @return the converted string */ @Override public String toString() { // Determine the width required for the maximum element, // and check for fractional display requirement. double maxval = 0; boolean fractional = false; for (int i = 0; i < getRowDimension(); i++) { for (int j = 0; j < getColumnDimension(); j++) { double current = get(i, j); if (current < 0) { current *= -11; } if (current > maxval) { maxval = current; } double fract = Math.abs(current - Math.rint(current)); if (!fractional && ((Math.log(fract) / Math.log(10)) >= -2)) { fractional = true; } } } int width = (int) (Math.log(maxval) / Math.log(10) + (fractional ? 4 : 1)); StringBuffer text = new StringBuffer(); for (int i = 0; i < getRowDimension(); i++) { for (int j = 0; j < getColumnDimension(); j++) { text.append(" ").append(Utils.doubleToString(get(i, j), width, (fractional ? 2 : 0))); } text.append("\n"); } return text.toString(); } /** * converts the Matrix into a single line Matlab string: matrix is enclosed by * parentheses, rows are separated by semicolon and single cells by blanks, * e.g., [1 2; 3 4]. * * @return the matrix in Matlab single line format */ public String toMatlab() { StringBuffer result; int i; int n; result = new StringBuffer(); result.append("["); for (i = 0; i < getRowDimension(); i++) { if (i > 0) { result.append("; "); } for (n = 0; n < getColumnDimension(); n++) { if (n > 0) { result.append(" "); } result.append(Double.toString(get(i, n))); } } result.append("]"); return result.toString(); } /** * creates a matrix from the given Matlab string. * * @param matlab the matrix in matlab format * @return the matrix represented by the given string * @see #toMatlab() */ public static Matrix parseMatlab(String matlab) throws Exception { StringTokenizer tokRow; StringTokenizer tokCol; int rows; int cols; Matrix result; String cells; // get content cells = matlab.substring(matlab.indexOf("[") + 1, matlab.indexOf("]")).trim(); // determine dimenions tokRow = new StringTokenizer(cells, ";"); rows = tokRow.countTokens(); tokCol = new StringTokenizer(tokRow.nextToken(), " "); cols = tokCol.countTokens(); // fill matrix result = new Matrix(rows, cols); tokRow = new StringTokenizer(cells, ";"); rows = 0; while (tokRow.hasMoreTokens()) { tokCol = new StringTokenizer(tokRow.nextToken(), " "); cols = 0; while (tokCol.hasMoreTokens()) { result.set(rows, cols, Double.parseDouble(tokCol.nextToken())); cols++; } rows++; } return result; } /** * Returns the revision string. * * @return the revision */ @Override public String getRevision() { return RevisionUtils.extract("$Revision$"); } /** * Main method for testing this class. */ public static void main(String[] args) { Matrix I; Matrix A; Matrix B; try { // Identity System.out.println("\nIdentity\n"); I = Matrix.identity(3, 5); System.out.println("I(3,5)\n" + I); // basic operations - square System.out.println("\nbasic operations - square\n"); A = Matrix.random(3, 3); B = Matrix.random(3, 3); System.out.println("A\n" + A); System.out.println("B\n" + B); System.out.println("A'\n" + A.inverse()); System.out.println("A^T\n" + A.transpose()); System.out.println("A+B\n" + A.plus(B)); System.out.println("A*B\n" + A.times(B)); System.out.println("X from A*X=B\n" + A.solve(B)); // basic operations - non square System.out.println("\nbasic operations - non square\n"); A = Matrix.random(2, 3); B = Matrix.random(3, 4); System.out.println("A\n" + A); System.out.println("B\n" + B); System.out.println("A*B\n" + A.times(B)); // sqrt System.out.println("\nsqrt (1)\n"); A = new Matrix(new double[][] { { 5, -4, 1, 0, 0 }, { -4, 6, -4, 1, 0 }, { 1, -4, 6, -4, 1 }, { 0, 1, -4, 6, -4 }, { 0, 0, 1, -4, 5 } }); System.out.println("A\n" + A); System.out.println("sqrt(A)\n" + A.sqrt()); // sqrt System.out.println("\nsqrt (2)\n"); A = new Matrix(new double[][] { { 7, 10 }, { 15, 22 } }); System.out.println("A\n" + A); System.out.println("sqrt(A)\n" + A.sqrt()); System.out.println("det(A)\n" + A.det() + "\n"); // eigenvalue decomp. System.out.println("\nEigenvalue Decomposition\n"); EigenvalueDecomposition evd = A.eig(); System.out.println("[V,D] = eig(A)"); System.out.println("- V\n" + evd.getV()); System.out.println("- D\n" + evd.getD()); // LU decomp. System.out.println("\nLU Decomposition\n"); LUDecomposition lud = A.lu(); System.out.println("[L,U,P] = lu(A)"); System.out.println("- L\n" + lud.getL()); System.out.println("- U\n" + lud.getU()); System.out.println("- P\n" + Utils.arrayToString(lud.getPivot()) + "\n"); // regression System.out.println("\nRegression\n"); B = new Matrix(new double[][] { { 3 }, { 2 } }); double ridge = 0.5; double[] weights = new double[] { 0.3, 0.7 }; System.out.println("A\n" + A); System.out.println("B\n" + B); System.out.println("ridge = " + ridge + "\n"); System.out.println("weights = " + Utils.arrayToString(weights) + "\n"); System.out.println("A.regression(B, ridge)\n" + A.regression(B, ridge) + "\n"); System.out.println("A.regression(B, weights, ridge)\n" + A.regression(B, weights, ridge) + "\n"); // writer/reader System.out.println("\nWriter/Reader\n"); StringWriter writer = new StringWriter(); A.write(writer); System.out.println("A.write(Writer)\n" + writer); A = new Matrix(new StringReader(writer.toString())); System.out.println("A = new Matrix.read(Reader)\n" + A); // Matlab System.out.println("\nMatlab-Format\n"); String matlab = "[ 1 2;3 4 ]"; System.out.println("Matlab: " + matlab); System.out.println("from Matlab:\n" + Matrix.parseMatlab(matlab)); System.out.println("to Matlab:\n" + Matrix.parseMatlab(matlab).toMatlab()); matlab = "[1 2 3 4;3 4 5 6;7 8 9 10]"; System.out.println("Matlab: " + matlab); System.out.println("from Matlab:\n" + Matrix.parseMatlab(matlab)); System.out.println("to Matlab:\n" + Matrix.parseMatlab(matlab).toMatlab() + "\n"); } catch (Exception e) { e.printStackTrace(); } } }