/* A library for mathematical operations on matrices of doubles. This file was automatically generated with a preprocessor, so that similar matrix operations are supported on ints, longs, floats, and doubles. Copyright (c) 1998-2014 The Regents of the University of California. All rights reserved. Permission is hereby granted, without written agreement and without license or royalty fees, to use, copy, modify, and distribute this software and its documentation for any purpose, provided that the above copyright notice and the following two paragraphs appear in all copies of this software. IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. PT_COPYRIGHT_VERSION_2 COPYRIGHTENDKEY */ package ptolemy.math; /////////////////////////////////////////////////////////////////// //// DoubleMatrixMath /** This class provides a library for mathematical operations on matrices of doubles. All calls expect matrix arguments to be non-null. In addition, all rows of the matrix are expected to have the same number of columns.
Some algorithms are from:
[1] Embree, Paul M. and Bruce Kimble. "C Language Algorithms for Digital Signal Processing," Prentice Hall, Englewood Cliffs, NJ, 1991. @author Jeff Tsay @version $Id: DoubleMatrixMath.java 70402 2014-10-23 00:52:20Z cxh $ @since Ptolemy II 1.0 @Pt.ProposedRating Yellow (ctsay) @Pt.AcceptedRating Yellow (ctsay) */ public class DoubleMatrixMath { // private constructor prevents construction of this class. private DoubleMatrixMath() { } /////////////////////////////////////////////////////////////////// //// public methods //// /** Return a new matrix that is constructed from the argument by * adding the second argument to every element. * @param matrix A matrix of doubles. * @param z The double number to add. * @return A new matrix of doubles. */ public static final double[][] add(double[][] matrix, double z) { double[][] returnValue = new double[_rows(matrix)][_columns(matrix)]; for (int i = 0; i < _rows(matrix); i++) { for (int j = 0; j < _columns(matrix); j++) { returnValue[i][j] = matrix[i][j] + z; } } return returnValue; } /** Return a new matrix that is constructed from the argument by * adding the second matrix to the first one. If the two * matrices are not the same size, throw an * IllegalArgumentException. * * @param matrix1 The first matrix of doubles. * @param matrix2 The second matrix of doubles. * @return A new matrix of doubles. */ public static final double[][] add(final double[][] matrix1, final double[][] matrix2) { _checkSameDimension("add", matrix1, matrix2); double[][] returnValue = new double[_rows(matrix1)][_columns(matrix1)]; for (int i = 0; i < _rows(matrix1); i++) { for (int j = 0; j < _columns(matrix1); j++) { returnValue[i][j] = matrix1[i][j] + matrix2[i][j]; } } return returnValue; } /** Return a new matrix that is a copy of the matrix argument. * @param matrix A matrix of doubles. * @return A new matrix of doubles. */ public static final double[][] allocCopy(final double[][] matrix) { return crop(matrix, 0, 0, _rows(matrix), _columns(matrix)); } /** Return a new array that is formed by applying an instance of a * DoubleBinaryOperation to each element in the input matrix, * using z as the left operand in all cases and the matrix elements * as the right operands (op.operate(z, matrix[i][j])). */ public static final double[][] applyBinaryOperation( DoubleBinaryOperation op, final double z, final double[][] matrix) { int rows = _rows(matrix); int columns = _columns(matrix); double[][] returnValue = new double[rows][columns]; for (int i = 0; i < rows; i++) { for (int j = 0; j < columns; j++) { returnValue[i][j] = op.operate(z, matrix[i][j]); } } return returnValue; } /** Return a new array that is formed by applying an instance of a * DoubleBinaryOperation to each element in the input matrix, * using the matrix elements as the left operands and z as the right * operand in all cases (op.operate(matrix[i][j], z)). */ public static final double[][] applyBinaryOperation( DoubleBinaryOperation op, final double[][] matrix, final double z) { int rows = _rows(matrix); int columns = _columns(matrix); double[][] returnValue = new double[rows][columns]; for (int i = 0; i < rows; i++) { for (int j = 0; j < columns; j++) { returnValue[i][j] = op.operate(matrix[i][j], z); } } return returnValue; } /** Return a new array that is formed by applying an instance of a * DoubleBinaryOperation to the two matrices, element by element, * using the elements of the first matrix as the left operands * and the elements of the second matrix as the right operands. * (op.operate(matrix1[i][j], matrix2[i][j])). If the matrices * are not the same size, throw an IllegalArgumentException. */ public static final double[][] applyBinaryOperation( DoubleBinaryOperation op, final double[][] matrix1, final double[][] matrix2) { int rows = _rows(matrix1); int columns = _columns(matrix1); _checkSameDimension("applyBinaryOperation", matrix1, matrix2); double[][] returnValue = new double[rows][columns]; for (int i = 0; i < rows; i++) { for (int j = 0; j < columns; j++) { returnValue[i][j] = op.operate(matrix1[i][j], matrix2[i][j]); } } return returnValue; } /** Return a new array that is formed by applying an instance of a * DoubleUnaryOperation to each element in the input matrix * (op.operate(matrix[i][j])). */ public static final double[][] applyUnaryOperation( final DoubleUnaryOperation op, final double[][] matrix) { int rows = _rows(matrix); int columns = _columns(matrix); double[][] returnValue = new double[rows][columns]; for (int i = 0; i < rows; i++) { for (int j = 0; j < columns; j++) { returnValue[i][j] = op.operate(matrix[i][j]); } } return returnValue; } /** Return a new matrix that is a sub-matrix of the input * matrix argument. The row and column from which to start * and the number of rows and columns to span are specified. * @param matrix A matrix of doubles. * @param rowStart An int specifying which row to start on. * @param colStart An int specifying which column to start on. * @param rowSpan An int specifying how many rows to copy. * @param colSpan An int specifying how many columns to copy. */ public static final double[][] crop(final double[][] matrix, final int rowStart, final int colStart, final int rowSpan, final int colSpan) { double[][] returnValue = new double[rowSpan][colSpan]; for (int i = 0; i < rowSpan; i++) { System.arraycopy(matrix[rowStart + i], colStart, returnValue[i], 0, colSpan); } return returnValue; } /** Return the determinant of a square matrix. * If the matrix is not square, throw an IllegalArgumentException. * This algorithm uses LU decomposition, and is taken from [1] */ public static final double determinant(final double[][] matrix) { _checkSquare("determinant", matrix); double[][] a; double det = 1.0; int n = _rows(matrix); a = allocCopy(matrix); for (int pivot = 0; pivot < n - 1; pivot++) { // find the biggest absolute pivot double big = Math.abs(a[pivot][pivot]); int swapRow = 0; // initialize for no swap for (int row = pivot + 1; row < n; row++) { double absElement = Math.abs(a[row][pivot]); if (absElement > big) { swapRow = row; big = absElement; } } // unless swapRow is still zero we must swap two rows if (swapRow != 0) { double[] aPtr = a[pivot]; a[pivot] = a[swapRow]; a[swapRow] = aPtr; // change sign of determinant because of swap det *= -a[pivot][pivot]; } else { // calculate the determinant by the product of the pivots det *= a[pivot][pivot]; } // if almost singular matrix, give up now // If almost singular matrix, give up now. if (Math.abs(det) <= Complex.EPSILON) { return det; } double pivotInverse = 1.0 / a[pivot][pivot]; for (int col = pivot + 1; col < n; col++) { a[pivot][col] *= pivotInverse; } for (int row = pivot + 1; row < n; row++) { double temp = a[row][pivot]; for (int col = pivot + 1; col < n; col++) { a[row][col] -= a[pivot][col] * temp; } } } // Last pivot, no reduction required. det *= a[n - 1][n - 1]; return det; } /** Return a new matrix that is constructed by placing the * elements of the input array on the diagonal of the square * matrix, starting from the top left corner down to the bottom * right corner. All other elements are zero. The size of of the * matrix is n x n, where n is the length of the input array. */ public static final double[][] diag(final double[] array) { int n = array.length; double[][] returnValue = new double[n][n]; // Assume the matrix is zero-filled. for (int i = 0; i < n; i++) { returnValue[i][i] = array[i]; } return returnValue; } /** Return a new matrix that is constructed from the argument by * dividing the second argument to every element. * @param matrix A matrix of doubles. * @param z The double number to divide. * @return A new matrix of doubles. */ public static final double[][] divide(double[][] matrix, double z) { double[][] returnValue = new double[_rows(matrix)][_columns(matrix)]; for (int i = 0; i < _rows(matrix); i++) { for (int j = 0; j < _columns(matrix); j++) { returnValue[i][j] = matrix[i][j] / z; } } return returnValue; } /** Return a new matrix that is constructed by element by element * division of the two matrix arguments. Each element of the * first matrix is divided by the corresponding element of the * second matrix. If the two matrices are not the same size, * throw an IllegalArgumentException. */ public static final double[][] divideElements(final double[][] matrix1, final double[][] matrix2) { int rows = _rows(matrix1); int columns = _columns(matrix1); _checkSameDimension("divideElements", matrix1, matrix2); double[][] returnValue = new double[rows][columns]; for (int i = 0; i < rows; i++) { for (int j = 0; j < columns; j++) { returnValue[i][j] = matrix1[i][j] / matrix2[i][j]; } } return returnValue; } /** Return a new array that is filled with the contents of the matrix. * The doubles are stored row by row, i.e. using the notation * (row, column), the entries of the array are in the following order * for a (m, n) matrix : * (0, 0), (0, 1), (0, 2), ... , (0, n-1), (1, 0), (1, 1), ..., (m-1)(n-1) * @param matrix A matrix of doubles. * @return A new array of doubles. */ public static final double[] fromMatrixToArray(final double[][] matrix) { return fromMatrixToArray(matrix, _rows(matrix), _columns(matrix)); } /** Return a new array that is filled with the contents of the matrix. * The maximum numbers of rows and columns to copy are specified so * that entries lying outside of this range can be ignored. The * maximum rows to copy cannot exceed the number of rows in the matrix, * and the maximum columns to copy cannot exceed the number of columns * in the matrix. * The doubles are stored row by row, i.e. using the notation * (row, column), the entries of the array are in the following order * for a matrix, limited to m rows and n columns : * (0, 0), (0, 1), (0, 2), ... , (0, n-1), (1, 0), (1, 1), ..., (m-1)(n-1) * @param matrix A matrix of doubles. * @return A new array of doubles. */ public static final double[] fromMatrixToArray(final double[][] matrix, int maxRow, int maxCol) { double[] returnValue = new double[maxRow * maxCol]; for (int i = 0; i < maxRow; i++) { System.arraycopy(matrix[i], 0, returnValue, i * maxCol, maxCol); } return returnValue; } /** Return a new matrix, which is defined by Aij = 1/(i+j+1), * the Hilbert matrix. The matrix is square with one * dimension specifier required. This matrix is useful because it always * has an inverse. */ public static final double[][] hilbert(final int dim) { double[][] returnValue = new double[dim][dim]; for (int i = 0; i < dim; i++) { for (int j = 0; j < dim; j++) { returnValue[i][j] = 1.0 / (i + j + 1); } } return returnValue; } /** Return an new identity matrix with the specified dimension. The * matrix is square, so only one dimension specifier is needed. * Note that this method does the same thing as identityDouble(), * but the latter is more useful in the expression language. */ public static final double[][] identity(final int dim) { double[][] returnValue = new double[dim][dim]; // we rely on the fact Java fills the allocated matrix with 0's for (int i = 0; i < dim; i++) { returnValue[i][i] = 1.0; } return returnValue; } /** Return an new identity matrix with the specified dimension. The * matrix is square, so only one dimension specifier is needed. */ public static final double[][] identityMatrixDouble(final int dim) { return identity(dim); } /** Return a new matrix that is constructed by inverting the input * matrix. If the input matrix is singular, throw an exception. * This method is from [1]. * @exception IllegalArgumentException If the matrix is singular. */ public static final double[][] inverse(final double[][] A) { _checkSquare("inverse", A); int n = _rows(A); double[][] Ai = allocCopy(A); // We depend on each of the elements being initialized to 0 int[] pivotFlag = new int[n]; int[] swapCol = new int[n]; int[] swapRow = new int[n]; int irow = 0; int icol = 0; for (int i = 0; i < n; i++) { // n iterations of pivoting // find the biggest pivot element double big = 0.0; for (int row = 0; row < n; row++) { if (pivotFlag[row] == 0) { for (int col = 0; col < n; col++) { if (pivotFlag[col] == 0) { double absElement = Math.abs(Ai[row][col]); if (absElement >= big) { big = absElement; irow = row; icol = col; } } } } } pivotFlag[icol]++; // swap rows to make this diagonal the biggest absolute pivot if (irow != icol) { for (int col = 0; col < n; col++) { double temp = Ai[irow][col]; Ai[irow][col] = Ai[icol][col]; Ai[icol][col] = temp; } } // store what we swapped swapRow[i] = irow; swapCol[i] = icol; // if the pivot is zero, the matrix is singular if (Ai[icol][icol] == 0.0) { throw new IllegalArgumentException( "Attempt to invert a singular matrix."); } // divide the row by the pivot double pivotInverse = 1.0 / Ai[icol][icol]; Ai[icol][icol] = 1.0; // pivot = 1 to avoid round off for (int col = 0; col < n; col++) { Ai[icol][col] *= pivotInverse; } // fix the other rows by subtracting for (int row = 0; row < n; row++) { if (row != icol) { double temp = Ai[row][icol]; Ai[row][icol] = 0.0; for (int col = 0; col < n; col++) { Ai[row][col] -= Ai[icol][col] * temp; } } } } // fix the effect of all the swaps for final answer for (int swap = n - 1; swap >= 0; swap--) { if (swapRow[swap] != swapCol[swap]) { for (int row = 0; row < n; row++) { double temp = Ai[row][swapRow[swap]]; Ai[row][swapRow[swap]] = Ai[row][swapCol[swap]]; Ai[row][swapCol[swap]] = temp; } } } return Ai; } /** Replace the destinationMatrix argument elements with the values of * the sourceMatrix argument. The destinationMatrix argument must be * large enough to hold all the values of sourceMatrix argument. * @param destinationMatrix A matrix of doubles, used as the destination. * @param sourceMatrix A matrix of doubles, used as the source. */ public static final void matrixCopy(final double[][] sourceMatrix, final double[][] destinationMatrix) { matrixCopy(sourceMatrix, 0, 0, destinationMatrix, 0, 0, _rows(sourceMatrix), _columns(sourceMatrix)); } /** Replace the destinationMatrix argument's values, in the specified row * and column range, with the sourceMatrix argument's values, starting * from specified row and column of the second matrix. * @param sourceMatrix A matrix of doubles, used as the destination. * @param sourceRowStart An int specifying the starting row of the source. * @param sourceColStart An int specifying the starting column of the * source. * @param destinationMatrix A matrix of doubles, used as the destination. * @param destinationRowStart An int specifying the starting row of the * destination. * @param destinationColumnStart An int specifying the starting column * of the destination. * @param rowSpan An int specifying how many rows to copy. * @param columnSpan An int specifying how many columns to copy. */ public static final void matrixCopy(final double[][] sourceMatrix, final int sourceRowStart, final int sourceColStart, final double[][] destinationMatrix, final int destinationRowStart, final int destinationColumnStart, final int rowSpan, final int columnSpan) { // We should verify the parameters here for (int i = 0; i < rowSpan; i++) { System.arraycopy(sourceMatrix[sourceRowStart + i], sourceColStart, destinationMatrix[destinationRowStart + i], destinationColumnStart, columnSpan); } } /** Return a new matrix that is constructed from the argument after * performing modulo operation by the second argument to every element. * @param matrix A matrix of doubles. * @param z The double number to divide. * @return A new matrix of doubles. */ public static final double[][] modulo(double[][] matrix, double z) { double[][] returnValue = new double[_rows(matrix)][_columns(matrix)]; for (int i = 0; i < _rows(matrix); i++) { for (int j = 0; j < _columns(matrix); j++) { returnValue[i][j] = matrix[i][j] % z; } } return returnValue; } /** Return a new matrix that is constructed by multiplying the matrix * by a scaleFactor. */ public static final double[][] multiply(final double[][] matrix, final double scaleFactor) { int rows = _rows(matrix); int columns = _columns(matrix); double[][] returnValue = new double[rows][columns]; for (int i = 0; i < rows; i++) { for (int j = 0; j < columns; j++) { returnValue[i][j] = matrix[i][j] * scaleFactor; } } return returnValue; } /** Return a new array that is constructed from the argument by * pre-multiplying the array (treated as a row vector) by a matrix. * The number of rows of the matrix must equal the number of elements * in the array. The returned array will have a length equal to the number * of columns of the matrix. */ public static final double[] multiply(final double[][] matrix, final double[] array) { int rows = _rows(matrix); int columns = _columns(matrix); if (rows != array.length) { throw new IllegalArgumentException( "preMultiply : array does not have the same number of " + "elements (" + array.length + ") as the number of rows " + "of the matrix (" + rows + ")"); } double[] returnValue = new double[columns]; for (int i = 0; i < columns; i++) { double sum = 0.0; for (int j = 0; j < rows; j++) { sum += matrix[j][i] * array[j]; } returnValue[i] = sum; } return returnValue; } /** Return a new array that is constructed from the argument by * post-multiplying the matrix by an array (treated as a row vector). * The number of columns of the matrix must equal the number of elements * in the array. The returned array will have a length equal to the number * of rows of the matrix. */ public static final double[] multiply(final double[] array, final double[][] matrix) { int rows = _rows(matrix); int columns = _columns(matrix); if (columns != array.length) { throw new IllegalArgumentException( "postMultiply() : array does not have the same number " + "of elements (" + array.length + ") as the number of " + "columns of the matrix (" + columns + ")"); } double[] returnValue = new double[rows]; for (int i = 0; i < rows; i++) { double sum = 0.0; for (int j = 0; j < columns; j++) { sum += matrix[i][j] * array[j]; } returnValue[i] = sum; } return returnValue; } /** Return a new matrix that is constructed from the argument by * multiplying the first matrix by the second one. * Note this operation is not commutative, * so care must be taken in the ordering of the arguments. * The number of columns of matrix1 * must equal the number of rows of matrix2. If matrix1 is of * size m x n, and matrix2 is of size n x p, the returned matrix * will have size m x p. *
Note that this method is different from the other multiply() * methods in that this method does not do pointwise multiplication. * * @see #multiplyElements(double[][], double[][]) * @param matrix1 The first matrix of doubles. * @param matrix2 The second matrix of doubles. * @return A new matrix of doubles. */ public static final double[][] multiply(double[][] matrix1, double[][] matrix2) { double[][] returnValue = new double[_rows(matrix1)][matrix2[0].length]; for (int i = 0; i < _rows(matrix1); i++) { for (int j = 0; j < matrix2[0].length; j++) { double sum = 0.0; for (int k = 0; k < matrix2.length; k++) { sum += matrix1[i][k] * matrix2[k][j]; } returnValue[i][j] = sum; } } return returnValue; } /** Return a new matrix that is constructed by element by element * multiplication of the two matrix arguments. If the two * matrices are not the same size, throw an * IllegalArgumentException. *
Note that this method does pointwise matrix multiplication.
* See {@link #multiply(double[][], double[][])} for standard
* matrix multiplication.
*/
public static final double[][] multiplyElements(final double[][] matrix1,
final double[][] matrix2) {
int rows = _rows(matrix1);
int columns = _columns(matrix1);
_checkSameDimension("multiplyElements", matrix1, matrix2);
double[][] returnValue = new double[rows][columns];
for (int i = 0; i < rows; i++) {
for (int j = 0; j < columns; j++) {
returnValue[i][j] = matrix1[i][j] * matrix2[i][j];
}
}
return returnValue;
}
/** Return a new matrix that is the additive inverse of the
* argument matrix.
*/
public static final double[][] negative(final double[][] matrix) {
int rows = _rows(matrix);
int columns = _columns(matrix);
double[][] returnValue = new double[rows][columns];
for (int i = 0; i < rows; i++) {
for (int j = 0; j < columns; j++) {
returnValue[i][j] = -matrix[i][j];
}
}
return returnValue;
}
/** Return a new matrix that is formed by orthogonalizing the
* columns of the input matrix (the column vectors are
* orthogonal). If not all columns are linearly independent, the
* output matrix will contain a column of zeros for all redundant
* input columns.
*/
public static final double[][] orthogonalizeColumns(final double[][] matrix) {
Object[] orthoInfo = _orthogonalizeRows(transpose(matrix));
return transpose((double[][]) orthoInfo[0]);
}
/** Return a new matrix that is formed by orthogonalizing the rows of the
* input matrix (the row vectors are orthogonal). If not all rows are
* linearly independent, the output matrix will contain a row of zeros
* for all redundant input rows.
*/
public static final double[][] orthogonalizeRows(final double[][] matrix) {
Object[] orthoInfo = _orthogonalizeRows(matrix);
return (double[][]) orthoInfo[0];
}
/** Return a new matrix that is formed by orthonormalizing the
* columns of the input matrix (the column vectors are orthogonal
* and have norm 1). If not all columns are linearly independent,
* the output matrix will contain a column of zeros for all
* redundant input columns.
*/
public static final double[][] orthonormalizeColumns(final double[][] matrix) {
return transpose(orthogonalizeRows(transpose(matrix)));
}
/** Return a new matrix that is formed by orthonormalizing the
* rows of the input matrix (the row vectors are orthogonal and
* have norm 1). If not all rows are linearly independent, the
* output matrix will contain a row of zeros for all redundant
* input rows.
*/
public static final double[][] orthonormalizeRows(final double[][] matrix) {
int rows = _rows(matrix);
Object[] orthoInfo = _orthogonalizeRows(matrix);
double[][] orthogonalMatrix = (double[][]) orthoInfo[0];
double[] oneOverNormSquaredArray = (double[]) orthoInfo[2];
for (int i = 0; i < rows; i++) {
orthogonalMatrix[i] = DoubleArrayMath.scale(orthogonalMatrix[i],
Math.sqrt(oneOverNormSquaredArray[i]));
}
return orthogonalMatrix;
}
/** Return a pair of matrices that are the decomposition of the
* input matrix (which must have linearly independent column
* vectors), which is m x n, into the matrix product of Q, which
* is m x n with orthonormal column vectors, and R, which is an
* invertible n x n upper triangular matrix.
*
* @param matrix The input matrix of doubles.
* @return The pair of newly allocated matrices of doubles,
* out[0] = Q, out[1] = R.
* @exception IllegalArgumentException if the columns vectors of the input
* matrix are not linearly independent.
*/
public static final double[][][] qr(final double[][] matrix) {
int columns = _columns(matrix);
/* Find an orthogonal basis using _orthogonalizeRows(). Note
* that _orthogonalizeRows() orthogonalizes row vectors, so
* we have use the transpose of input matrix to orthogonalize
* its columns vectors. The output will be the transpose of
* Q.
*/
Object[] orthoRowInfo = _orthogonalizeRows(transpose(matrix));
double[][] qT = (double[][]) orthoRowInfo[0];
// get the dot product matrix, dp[j][i] =
*
* Orthogonalization is done with the Gram-Schmidt process.
*/
protected static final Object[] _orthogonalizeRows(
final double[][] rowArrays) {
int rows = rowArrays.length;
int columns = rowArrays[0].length;
int nullity = 0;
double[][] orthogonalMatrix = new double[rows][];
double[] oneOverNormSquaredArray = new double[rows];
// A matrix containing the dot products of the input row
// vectors and output row vectors, dotProductMatrix[j][i] =
//
*
*
*
*
*
*