Java tutorial
/* * Licensed to the Apache Software Foundation (ASF) under one * or more contributor license agreements. See the NOTICE file * distributed with this work for additional information * regarding copyright ownership. The ASF licenses this file * to you under the Apache License, Version 2.0 (the * "License"); you may not use this file except in compliance * with the License. You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, * software distributed under the License is distributed on an * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY * KIND, either express or implied. See the License for the * specific language governing permissions and limitations * under the License. */ package hivemall.utils.math; import hivemall.utils.collections.lists.DoubleArrayList; import hivemall.utils.lang.Preconditions; import java.util.Arrays; import javax.annotation.Nonnull; import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.ArrayRealVector; import org.apache.commons.math3.linear.BlockRealMatrix; import org.apache.commons.math3.linear.DecompositionSolver; import org.apache.commons.math3.linear.DefaultRealMatrixPreservingVisitor; import org.apache.commons.math3.linear.LUDecomposition; import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.RealMatrixPreservingVisitor; import org.apache.commons.math3.linear.RealVector; import org.apache.commons.math3.linear.SingularMatrixException; import org.apache.commons.math3.linear.SingularValueDecomposition; public final class MatrixUtils { private MatrixUtils() { } /** * Solve Yule-walker equation by Levinson-Durbin Recursion. * * <pre> * R_j = _{i=1}^{k} A_i R_{j-i} where j = 1..k, R_{-i} = R'_i * </pre> * * @see <a href= * "http://www.emptyloop.com/technotes/a%20tutorial%20on%20linear%20prediction%20and%20levinson-durbin.pdf">Cedrick * Collomb: A tutorial on linear prediction and Levinson-Durbin</a> * @param R autocovariance where |R| >= order * @param A coefficient to be solved where |A| >= order + 1 * @return E variance of prediction error */ @Nonnull public static double[] aryule(@Nonnull final double[] R, @Nonnull final double[] A, final int order) { Preconditions.checkArgument(R.length > order, "|R| MUST be greater than or equals to " + order + ": " + R.length); Preconditions.checkArgument(A.length >= order + 1, "|A| MUST be greater than or equals to " + (order + 1) + ": " + A.length); final double[] E = new double[order + 1]; A[0] = 1.0d; E[0] = R[0]; for (int k = 0; k < order; k++) { double lambda = 0.d; for (int j = 0; j <= k; j++) { lambda -= A[j] * R[k + 1 - j]; } final double Ek = E[k]; if (Ek == 0.d) { lambda = 0.d; } else { lambda /= Ek; } for (int n = 0, last = (k + 1) / 2; n <= last; n++) { final int i = k + 1 - n; double tmp = A[i] + lambda * A[n]; A[n] += lambda * A[i]; A[i] = tmp; } E[k + 1] = Ek * (1.0d - lambda * lambda); } for (int i = 0; i < order + 1; i++) { A[i] = -A[i]; } return E; } @Deprecated @Nonnull public static double[] aryule2(@Nonnull final double[] R, @Nonnull final double[] A, final int order) { Preconditions.checkArgument(R.length > order, "|C| MUST be greater than or equals to " + order + ": " + R.length); Preconditions.checkArgument(A.length >= order + 1, "|A| MUST be greater than or equals to " + (order + 1) + ": " + A.length); final double[] E = new double[order + 1]; A[0] = E[0] = 1.0d; A[1] = -R[1] / R[0]; E[1] = R[0] + R[1] * A[1]; for (int k = 1; k < order; k++) { double lambda = 0.d; for (int j = 0; j <= k; j++) { lambda -= A[j] * R[k + 1 - j]; } lambda /= E[k]; final double[] U = new double[k + 2]; final double[] V = new double[k + 2]; U[0] = 1.0; // V[0] = 0.0; for (int i = 1; i <= k; i++) { U[i] = A[i]; V[k + 1 - i] = A[i]; } V[k + 1] = 1.0; // U[k + 1] = 0.0; for (int i = 0, threshold = k + 2; i < threshold; i++) { A[i] = U[i] + lambda * V[i]; } E[k + 1] = E[k] * (1.0d - lambda * lambda); } for (int i = 0; i < order + 1; i++) { A[i] = -A[i]; } return E; } /** * Fit an AR(order) model using the Burg's method. * * @see https://searchcode.com/codesearch/view/9503568/ * @param X data vector to estimate where |X| >= order * @param A coefficient to be solved where |A| >= order + 1 * @return E variance of white noise */ @Nonnull public static double[] arburg(@Nonnull final double[] X, @Nonnull final double[] A, final int order) { Preconditions.checkArgument(X.length > order, "|X| MUST be greater than or equals to " + order + ": " + X.length); Preconditions.checkArgument(A.length >= order + 1, "|A| MUST be greater than or equals to " + (order + 1) + ": " + A.length); final int nDataPoints = X.length; final double[] E = new double[order + 1]; E[0] = 0.0d; for (int i = 0; i < nDataPoints; i++) { E[0] += X[i] * X[i]; } // f and b are the forward and backward error sequences int currentErrorSequenceSize = nDataPoints - 1; double[] F = new double[currentErrorSequenceSize]; double[] B = new double[currentErrorSequenceSize]; for (int i = 0; i < currentErrorSequenceSize; i++) { F[i] = X[i + 1]; B[i] = X[i]; } A[0] = 1.0d; // remaining stages i=2 to p for (int i = 0; i < order; i++) { // get the i-th reflection coefficient double numerator = 0.0d; double denominator = 0.0d; for (int j = 0; j < currentErrorSequenceSize; j++) { numerator += F[j] * B[j]; denominator += F[j] * F[j] + B[j] * B[j]; } numerator *= 2.0d; double g = 0.0d; if (denominator != 0.0d) { g = numerator / denominator; } // generate next filter order final double[] prevA = new double[i]; for (int j = 0; j < i; j++) { // No need to copy A[0] = 1.0 prevA[j] = A[j + 1]; } A[1] = g; for (int j = 1; j < i + 1; j++) { A[j + 1] = prevA[j - 1] - g * prevA[i - j]; } // keep track of the error E[i + 1] = E[i] * (1 - g * g); // update the prediction error sequences final double[] prevF = new double[currentErrorSequenceSize]; for (int j = 0; j < currentErrorSequenceSize; j++) { prevF[j] = F[j]; } final int nextErrorSequenceSize = nDataPoints - i - 2; for (int j = 0; j < nextErrorSequenceSize; j++) { F[j] = prevF[j + 1] - g * B[j + 1]; B[j] = B[j] - g * prevF[j]; } currentErrorSequenceSize = nextErrorSequenceSize; } for (int i = 1, mid = order / 2 + 1; i < mid; i++) { // Reverse 1..(order - 1)-th elements by swapping final double tmp = A[i]; A[i] = A[order + 1 - i]; A[order + 1 - i] = tmp; } for (int i = 0; i < order + 1; i++) { A[i] = -A[i]; } return E; } /** * Construct a Toeplitz matrix. */ @Nonnull public static RealMatrix[][] toeplitz(@Nonnull final RealMatrix[] c, final int dim) { Preconditions.checkArgument(dim >= 1, "Invalid dimension: " + dim); Preconditions.checkArgument(c.length >= dim, "|c| must be greater than " + dim + ": " + c.length); /* * Toeplitz matrix (symmetric, invertible, k*dimensions by k*dimensions) * * /C_0 |C_1' |C_2' | . . . |C_{k-1}' \ * |--------+--------+--------+ +---------| * |C_1 |C_0 |C_1' | . | * |--------+--------+--------+ . | * |C_2 |C_1 |C_0 | . | * |--------+--------+--------+ | * | . . | * | . . | * | . . | * |--------+ +--------| * \C_{k-1} | . . . |C_0 / */ final RealMatrix c0 = c[0]; final RealMatrix[][] toeplitz = new RealMatrix[dim][dim]; for (int row = 0; row < dim; row++) { toeplitz[row][row] = c0; for (int col = 0; col < dim; col++) { if (row < col) { toeplitz[row][col] = c[col - row].transpose(); } else if (row > col) { toeplitz[row][col] = c[row - col]; } } } return toeplitz; } /** * Construct a Toeplitz matrix. */ @Nonnull public static double[][] toeplitz(@Nonnull final double[] c) { return toeplitz(c, c.length); } /** * Construct a Toeplitz matrix. */ @Nonnull public static double[][] toeplitz(@Nonnull final double[] c, final int dim) { Preconditions.checkArgument(dim >= 1, "Invalid dimension: " + dim); Preconditions.checkArgument(c.length >= dim, "|c| must be greater than " + dim + ": " + c.length); /* * Toeplitz matrix (symmetric, invertible, k*dimensions by k*dimensions) * * /C_0 |C_1' |C_2' | . . . |C_{k-1}' \ * |--------+--------+--------+ +---------| * |C_1 |C_0 |C_1' | . | * |--------+--------+--------+ . | * |C_2 |C_1 |C_0 | . | * |--------+--------+--------+ | * | . . | * | . . | * | . . | * |--------+ +--------| * \C_{k-1} | . . . |C_0 / */ final double c0 = c[0]; final double[][] toeplitz = new double[dim][dim]; for (int row = 0; row < dim; row++) { toeplitz[row][row] = c0; for (int col = 0; col < dim; col++) { if (row < col) { toeplitz[row][col] = c[col - row]; } else if (row > col) { toeplitz[row][col] = c[row - col]; } } } return toeplitz; } @Nonnull public static double[] flatten(@Nonnull final RealMatrix[][] grid) { Preconditions.checkArgument(grid.length >= 1, "The number of rows must be greater than 1"); Preconditions.checkArgument(grid[0].length >= 1, "The number of cols must be greater than 1"); final int rows = grid.length; final int cols = grid[0].length; RealMatrix grid00 = grid[0][0]; Preconditions.checkNotNull(grid00); int cellRows = grid00.getRowDimension(); int cellCols = grid00.getColumnDimension(); final DoubleArrayList list = new DoubleArrayList(rows * cols * cellRows * cellCols); final RealMatrixPreservingVisitor visitor = new DefaultRealMatrixPreservingVisitor() { @Override public void visit(int row, int column, double value) { list.add(value); } }; for (int row = 0; row < rows; row++) { for (int col = 0; col < cols; col++) { RealMatrix cell = grid[row][col]; cell.walkInRowOrder(visitor); } } return list.toArray(); } @Nonnull public static double[] flatten(@Nonnull final RealMatrix[] grid) { Preconditions.checkArgument(grid.length >= 1, "The number of rows must be greater than 1"); final int rows = grid.length; RealMatrix grid0 = grid[0]; Preconditions.checkNotNull(grid0); int cellRows = grid0.getRowDimension(); int cellCols = grid0.getColumnDimension(); final DoubleArrayList list = new DoubleArrayList(rows * cellRows * cellCols); final RealMatrixPreservingVisitor visitor = new DefaultRealMatrixPreservingVisitor() { @Override public void visit(int row, int column, double value) { list.add(value); } }; for (int row = 0; row < rows; row++) { RealMatrix cell = grid[row]; cell.walkInRowOrder(visitor); } return list.toArray(); } @Nonnull public static RealMatrix[] unflatten(@Nonnull final double[] data, final int rows, final int cols, final int len) { final RealMatrix[] grid = new RealMatrix[len]; int offset = 0; for (int k = 0; k < len; k++) { RealMatrix cell = new BlockRealMatrix(rows, cols); grid[k] = cell; for (int i = 0; i < rows; i++) { for (int j = 0; j < cols; j++) { if (offset >= data.length) { throw new IndexOutOfBoundsException( "Offset " + offset + " exceeded data.length " + data.length); } double value = data[offset]; cell.setEntry(i, j, value); offset++; } } } if (offset != data.length) { throw new IllegalArgumentException("Invalid data for unflatten"); } return grid; } @Nonnull public static RealMatrix combinedMatrices(@Nonnull final RealMatrix[][] grid, final int dimensions) { Preconditions.checkArgument(grid.length >= 1, "The number of rows must be greater than 1"); Preconditions.checkArgument(grid[0].length >= 1, "The number of cols must be greater than 1"); Preconditions.checkArgument(dimensions > 0, "Dimension should be more than 0: ", dimensions); final int rows = grid.length; final int cols = grid[0].length; final RealMatrix combined = new BlockRealMatrix(rows * dimensions, cols * dimensions); for (int row = 0; row < grid.length; row++) { for (int col = 0; col < grid[row].length; col++) { combined.setSubMatrix(grid[row][col].getData(), row * dimensions, col * dimensions); } } return combined; } @Nonnull public static RealMatrix combinedMatrices(@Nonnull final RealMatrix[] grid) { Preconditions.checkArgument(grid.length >= 1, "The number of rows must be greater than 0: " + grid.length); final int rows = grid.length; final int rowDims = grid[0].getRowDimension(); final int colDims = grid[0].getColumnDimension(); final RealMatrix combined = new BlockRealMatrix(rows * rowDims, colDims); for (int row = 0; row < grid.length; row++) { RealMatrix cell = grid[row]; Preconditions.checkArgument(cell.getRowDimension() == rowDims, "Mismatch in row dimensions at row ", row); Preconditions.checkArgument(cell.getColumnDimension() == colDims, "Mismatch in col dimensions at row ", row); combined.setSubMatrix(cell.getData(), row * rowDims, 0); } return combined; } @Nonnull public static RealMatrix inverse(@Nonnull final RealMatrix m) throws SingularMatrixException { return inverse(m, true); } @Nonnull public static RealMatrix inverse(@Nonnull final RealMatrix m, final boolean exact) throws SingularMatrixException { LUDecomposition LU = new LUDecomposition(m); DecompositionSolver solver = LU.getSolver(); final RealMatrix inv; if (exact || solver.isNonSingular()) { inv = solver.getInverse(); } else { SingularValueDecomposition SVD = new SingularValueDecomposition(m); inv = SVD.getSolver().getInverse(); } return inv; } public static double det(@Nonnull final RealMatrix m) { LUDecomposition LU = new LUDecomposition(m); return LU.getDeterminant(); } /** * Return a 2-D array with ones on the diagonal and zeros elsewhere. */ @Nonnull public static double[][] eye(int n) { final double[][] eye = new double[n][n]; for (int i = 0; i < n; i++) { eye[i][i] = 1; } return eye; } /** * L = A x R * * @return a matrix A that minimizes A x R - L */ @Nonnull public static RealMatrix solve(@Nonnull final RealMatrix L, @Nonnull final RealMatrix R) throws SingularMatrixException { return solve(L, R, true); } /** * L = A x R * * @return a matrix A that minimizes A x R - L */ @Nonnull public static RealMatrix solve(@Nonnull final RealMatrix L, @Nonnull final RealMatrix R, final boolean exact) throws SingularMatrixException { LUDecomposition LU = new LUDecomposition(R); DecompositionSolver solver = LU.getSolver(); final RealMatrix A; if (exact || solver.isNonSingular()) { A = LU.getSolver().solve(L); } else { SingularValueDecomposition SVD = new SingularValueDecomposition(R); A = SVD.getSolver().solve(L); } return A; } /** * Find the first singular vector/value of a matrix A based on the Power method. * * @see http * ://www.cs.yale.edu/homes/el327/datamining2013aFiles/07_singular_value_decomposition.pdf * @param A target matrix * @param x0 initial vector * @param nIter number of iterations for the Power method * @param u 1st left singular vector * @param v 1st right singular vector * @return 1st singular value */ public static double power1(@Nonnull final RealMatrix A, @Nonnull final double[] x0, final int nIter, @Nonnull final double[] u, @Nonnull final double[] v) { Preconditions.checkArgument(A.getColumnDimension() == x0.length, "Column size of A and length of x should be same"); Preconditions.checkArgument(A.getRowDimension() == u.length, "Row size of A and length of u should be same"); Preconditions.checkArgument(x0.length == v.length, "Length of x and u should be same"); Preconditions.checkArgument(nIter >= 1, "Invalid number of iterations: " + nIter); RealMatrix AtA = A.transpose().multiply(A); RealVector x = new ArrayRealVector(x0); for (int i = 0; i < nIter; i++) { x = AtA.operate(x); } double xNorm = x.getNorm(); for (int i = 0, n = v.length; i < n; i++) { v[i] = x.getEntry(i) / xNorm; } RealVector Av = new ArrayRealVector(A.operate(v)); double s = Av.getNorm(); for (int i = 0, n = u.length; i < n; i++) { u[i] = Av.getEntry(i) / s; } return s; } /** * Lanczos tridiagonalization for a symmetric matrix C to make s * s tridiagonal matrix T. * * @see http://www.cas.mcmaster.ca/~qiao/publications/spie05.pdf * @param C target symmetric matrix * @param a initial vector * @param T result is stored here */ public static void lanczosTridiagonalization(@Nonnull final RealMatrix C, @Nonnull final double[] a, @Nonnull final RealMatrix T) { Preconditions.checkArgument(Arrays.deepEquals(C.getData(), C.transpose().getData()), "Target matrix C must be a symmetric matrix"); Preconditions.checkArgument(C.getColumnDimension() == a.length, "Column size of A and length of a should be same"); Preconditions.checkArgument(T.getRowDimension() == T.getColumnDimension(), "T must be a square matrix"); int s = T.getRowDimension(); // initialize T with zeros T.setSubMatrix(new double[s][s], 0, 0); RealVector a0 = new ArrayRealVector(a.length); RealVector r = new ArrayRealVector(a); double beta0 = 1.d; for (int i = 0; i < s; i++) { RealVector a1 = r.mapDivide(beta0); RealVector Ca1 = C.operate(a1); double alpha1 = a1.dotProduct(Ca1); r = Ca1.add(a1.mapMultiply(-1.d * alpha1)).add(a0.mapMultiply(-1.d * beta0)); double beta1 = r.getNorm(); T.setEntry(i, i, alpha1); if (i - 1 >= 0) { T.setEntry(i, i - 1, beta0); } if (i + 1 < s) { T.setEntry(i, i + 1, beta1); } a0 = a1.copy(); beta0 = beta1; } } /** * QR decomposition for a tridiagonal matrix T. * * @see https://gist.github.com/lightcatcher/8118181 * @see http://www.ericmart.in/blog/optimizing_julia_tridiag_qr * @param T target tridiagonal matrix * @param R output matrix for R which is the same shape as T * @param Qt output matrix for Q.T which is the same shape an T */ public static void tridiagonalQR(@Nonnull final RealMatrix T, @Nonnull final RealMatrix R, @Nonnull final RealMatrix Qt) { int n = T.getRowDimension(); Preconditions.checkArgument(n == R.getRowDimension() && n == R.getColumnDimension(), "T and R must be the same shape"); Preconditions.checkArgument(n == Qt.getRowDimension() && n == Qt.getColumnDimension(), "T and Qt must be the same shape"); // initial R = T R.setSubMatrix(T.getData(), 0, 0); // initial Qt = identity Qt.setSubMatrix(eye(n), 0, 0); for (int i = 0; i < n - 1; i++) { // Householder projection for a vector x // https://en.wikipedia.org/wiki/Householder_transformation RealVector x = T.getSubMatrix(i, i + 1, i, i).getColumnVector(0); x = unitL2norm(x); RealMatrix subR = R.getSubMatrix(i, i + 1, 0, n - 1); R.setSubMatrix(subR.subtract(x.outerProduct(subR.preMultiply(x)).scalarMultiply(2)).getData(), i, 0); RealMatrix subQt = Qt.getSubMatrix(i, i + 1, 0, n - 1); Qt.setSubMatrix(subQt.subtract(x.outerProduct(subQt.preMultiply(x)).scalarMultiply(2)).getData(), i, 0); } } @Nonnull static RealVector unitL2norm(@Nonnull final RealVector x) { double x0 = x.getEntry(0); double sign = MathUtils.sign(x0); x.setEntry(0, x0 + sign * x.getNorm()); return x.unitVector(); } /** * Find eigenvalues and eigenvectors of given tridiagonal matrix T. * * @see http://web.csulb.edu/~tgao/math423/s94.pdf * @see http://stats.stackexchange.com/questions/20643/finding-matrix-eigenvectors-using-qr- * decomposition * @param T target tridiagonal matrix * @param nIter number of iterations for the QR method * @param eigvals eigenvalues are stored here * @param eigvecs eigenvectors are stored here */ public static void tridiagonalEigen(@Nonnull final RealMatrix T, @Nonnull final int nIter, @Nonnull final double[] eigvals, @Nonnull final RealMatrix eigvecs) { Preconditions.checkArgument(Arrays.deepEquals(T.getData(), T.transpose().getData()), "Target matrix T must be a symmetric (tridiagonal) matrix"); Preconditions.checkArgument(eigvecs.getRowDimension() == eigvecs.getColumnDimension(), "eigvecs must be a square matrix"); Preconditions.checkArgument(T.getRowDimension() == eigvecs.getRowDimension(), "T and eigvecs must be the same shape"); Preconditions.checkArgument(eigvals.length == eigvecs.getRowDimension(), "Number of eigenvalues and eigenvectors must be same"); int nEig = eigvals.length; // initialize eigvecs as an identity matrix eigvecs.setSubMatrix(eye(nEig), 0, 0); RealMatrix T_ = T.copy(); for (int i = 0; i < nIter; i++) { // QR decomposition for the tridiagonal matrix T RealMatrix R = new Array2DRowRealMatrix(nEig, nEig); RealMatrix Qt = new Array2DRowRealMatrix(nEig, nEig); tridiagonalQR(T_, R, Qt); RealMatrix Q = Qt.transpose(); T_ = R.multiply(Q); eigvecs.setSubMatrix(eigvecs.multiply(Q).getData(), 0, 0); } // diagonal elements correspond to the eigenvalues for (int i = 0; i < nEig; i++) { eigvals[i] = T_.getEntry(i, i); } } }