/*
* 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.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 http://www.emptyloop.com/technotes/a%20tutorial%20on%20linear%20prediction%20and%20levinson-durbin.pdf
* @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 greather 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 greather 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 greather than 1");
Preconditions.checkArgument(grid[0].length >= 1,
"The number of cols must be greather 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 greather 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 greather than 1");
Preconditions.checkArgument(grid[0].length >= 1,
"The number of cols must be greather 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 greather 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);
}
}
}