no.uib.cipr.matrix.Matrix Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of mt-java Show documentation
Show all versions of mt-java Show documentation
Matrix data structures, linear solvers, least squares methods, eigenvalue,
and singular value decompositions. For larger random dense matrices (above ~ 350 x 350)
matrix-matrix multiplication C = A.B is about 50% faster than MTJ.
The newest version!
/*
* Copyright (C) 2003-2006 Bjørn-Ove Heimsund
*
* This file is part of MTJ.
*
* This library is free software; you can redistribute it and/or modify it
* under the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation; either version 2.1 of the License, or (at your
* option) any later version.
*
* This library is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License
* for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this library; if not, write to the Free Software Foundation,
* Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
package no.uib.cipr.matrix;
/**
* Basic matrix interface. It holds doubles in a rectangular 2D
* array, and it is used alongside Vector in numerical
* computations. Implementing classes decides on the actual storage.
*
* Basic operations
*
* Use numRows and numColumns to get the basic size of
* a matrix. get(int,int) gets an element, and there are
* corresponding set(int,int,double) and
* add(int,int,double) methods as well. Note that matrix indices
* are zero-based (typical for Java and C). This means that the row-indices
* range from 0 to numRows-1, likewise for the columns. It is legal
* to have numRows or numColumns equal zero.
*
*
* Other basic operations are zero which zeros all the entries of
* the matrix, which can be cheaper than either zeroing the matrix manually, or
* creating a new matrix, and the operation copy which creates a
* deep copy of the matrix. This copy has separate storage, but starts with the
* same contents as the current matrix.
*
*
* Iterators
*
* The matrix interface extends Iterable, and the iterator returns
* a MatrixEntry which contains current index and entry value. Note
* that the iterator may skip non-zero entries. Using an iterator, many simple
* and efficient algorithms can be created. The iterator also permits changing
* values in the matrix, however only non-zero entries can be changed.
*
*
* Basic linear algebra
*
* A large selection of basic linear algebra operations are available. To ensure
* high efficiency, little or no internal memory allocation is done, and the
* user is required to supply the output arguments.
*
*
* The operations available include:
*
*
* - Additions
* - Matrices can be added to each other, even if their underlying matrix
* structures are different.
* - Multiplications
* - A matrix can be multiplied with vectors and other matrices. For increased
* efficiency, a multiplication can be combined with addition and scaling, and
* transpose matrix multiplications are also available.
* - Rank-updates
* - A matrix can be efficiently updated using low-rank updates. The updates
* can be contained in both matrices or vectors.
* - Transpositions
* - In-place transpositions of square matrices is supported, and the
* transpose of a matrix can be stored in another matrix of compatible size
* (possibly non-rectangular)
* - Solvers
* - Many dense and structured sparse matrices have fast, direct solvers, and
* can be used to solve linear systems without creating a factorization. These
* solvers are typically backed by subroutines in LAPACK
*
*/
public interface Matrix extends Iterable {
/**
* Number of rows in the matrix
*/
int numRows();
/**
* Number of columns in the matrix
*/
int numColumns();
/**
* Returns true if the matrix is square
*/
boolean isSquare();
/**
* A(row,column) = value
*/
void set(int row, int column, double value);
/**
* A(row,column) += value
*/
void add(int row, int column, double value);
/**
* Returns A(row,column)
*/
double get(int row, int column);
/**
* Creates a deep copy of the matrix
*
* @return A
*/
Matrix copy();
/**
* Zeros all the entries in the matrix, while preserving any underlying
* structure. Useful for general, unstructured matrices.
*
* @return A
*/
Matrix zero();
/**
* y = A*x
*
* @param x
* Vector of size A.numColumns()
* @param y
* Vector of size A.numRows()
* @return y
*/
Vector mult(Vector x, Vector y);
/**
* y = alpha*A*x
*
* @param x
* Vector of size A.numColumns()
* @param y
* Vector of size A.numRows()
* @return y
*/
Vector mult(double alpha, Vector x, Vector y);
/**
* y = A*x + y
*
* @param x
* Vector of size A.numColumns()
* @param y
* Vector of size A.numRows()
* @return y
*/
Vector multAdd(Vector x, Vector y);
/**
* y = alpha*A*x + y
*
* @param x
* Vector of size A.numColumns()
* @param y
* Vector of size A.numRows()
* @return y
*/
Vector multAdd(double alpha, Vector x, Vector y);
/**
* y = AT*x
*
* @param x
* Vector of size A.numRows()
* @param y
* Vector of size A.numColumns()
* @return y
*/
Vector transMult(Vector x, Vector y);
/**
* y = alpha*AT*x
*
* @param x
* Vector of size A.numRows()
* @param y
* Vector of size A.numColumns()
* @return y
*/
Vector transMult(double alpha, Vector x, Vector y);
/**
* y = AT*x + y
*
* @param x
* Vector of size A.numRows()
* @param y
* Vector of size A.numColumns()
* @return y
*/
Vector transMultAdd(Vector x, Vector y);
/**
* y = alpha*AT*x + y
*
* @param x
* Vector of size A.numRows()
* @param y
* Vector of size A.numColumns()
* @return y
*/
Vector transMultAdd(double alpha, Vector x, Vector y);
/**
* x = A\b. Not all matrices support this operation, those that
* do not throw UnsupportedOperationException. Note that it is
* often more efficient to use a matrix decomposition and its associated
* solver
*
* @param b
* Vector of size A.numRows()
* @param x
* Vector of size A.numColumns()
* @return x
* @throws MatrixSingularException
* If the matrix is singular
* @throws MatrixNotSPDException
* If the solver assumes that the matrix is symmetrical,
* positive definite, but that that property does not hold
*/
Vector solve(Vector b, Vector x) throws MatrixSingularException,
MatrixNotSPDException;
/**
* x = AT\b. Not all matrices support this
* operation, those that do not throw
* UnsupportedOperationException. Note that it is often more
* efficient to use a matrix decomposition and its associated solver
*
* @param b
* Vector of size A.numColumns()
* @param x
* Vector of size A.numRows()
* @return x
* @throws MatrixSingularException
* If the matrix is singular
* @throws MatrixNotSPDException
* If the solver assumes that the matrix is symmetrical,
* positive definite, but that that property does not hold
*/
Vector transSolve(Vector b, Vector x) throws MatrixSingularException,
MatrixNotSPDException;
/**
* A = x*xT + A. The matrix must be square, and the
* vector of the same length
*
* @return A
*/
Matrix rank1(Vector x);
/**
* A = alpha*x*xT + A. The matrix must be square,
* and the vector of the same length
*
* @return A
*/
Matrix rank1(double alpha, Vector x);
/**
* A = x*yT + A. The matrix must be square, and the
* vectors of the same length
*
* @return A
*/
Matrix rank1(Vector x, Vector y);
/**
* A = alpha*x*yT + A. The matrix must be square,
* and the vectors of the same length
*
* @return A
*/
Matrix rank1(double alpha, Vector x, Vector y);
/**
* A = x*yT + y*xT + A. The matrix must
* be square, and the vectors of the same length
*
* @return A
*/
Matrix rank2(Vector x, Vector y);
/**
* A = alpha*x*yT + alpha*y*xT + A. The
* matrix must be square, and the vectors of the same length
*
* @return A
*/
Matrix rank2(double alpha, Vector x, Vector y);
/**
* C = A*B
*
* @param B
* Matrix such that B.numRows() == A.numColumns()
* and B.numColumns() == C.numColumns()
* @param C
* Matrix such that C.numRows() == A.numRows() and
* B.numColumns() == C.numColumns()
* @return C
*/
Matrix mult(Matrix B, Matrix C);
/**
* C = alpha*A*B
*
* @param B
* Matrix such that B.numRows() == A.numColumns()
* and B.numColumns() == C.numColumns()
* @param C
* Matrix such that C.numRows() == A.numRows() and
* B.numColumns() == C.numColumns()
* @return C
*/
Matrix mult(double alpha, Matrix B, Matrix C);
/**
* C = A*B + C
*
* @param B
* Matrix such that B.numRows() == A.numColumns()
* and B.numColumns() == C.numColumns()
* @param C
* Matrix such that C.numRows() == A.numRows() and
* B.numColumns() == C.numColumns()
* @return C
*/
Matrix multAdd(Matrix B, Matrix C);
/**
* C = alpha*A*B + C
*
* @param B
* Matrix such that B.numRows() == A.numColumns()
* and B.numColumns() == C.numColumns()
* @param C
* Matrix such that C.numRows() == A.numRows() and
* B.numColumns() == C.numColumns()
* @return C
*/
Matrix multAdd(double alpha, Matrix B, Matrix C);
/**
* C = AT*B
*
* @param B
* Matrix such that B.numRows() == A.numRows() and
* B.numColumns() == C.numColumns()
* @param C
* Matrix such that C.numRows() == A.numColumns()
* and B.numColumns() == C.numColumns()
* @return C
*/
Matrix transAmult(Matrix B, Matrix C);
/**
* C = alpha*AT*B
*
* @param B
* Matrix such that B.numRows() == A.numRows() and
* B.numColumns() == C.numColumns()
* @param C
* Matrix such that C.numRows() == A.numColumns()
* and B.numColumns() == C.numColumns()
* @return C
*/
Matrix transAmult(double alpha, Matrix B, Matrix C);
/**
* C = AT*B + C
*
* @param B
* Matrix such that B.numRows() == A.numRows() and
* B.numColumns() == C.numColumns()
* @param C
* Matrix such that C.numRows() == A.numColumns()
* and B.numColumns() == C.numColumns()
* @return C
*/
Matrix transAmultAdd(Matrix B, Matrix C);
/**
* C = alpha*AT*B + C
*
* @param B
* Matrix such that B.numRows() == A.numRows() and
* B.numColumns() == C.numColumns()
* @param C
* Matrix such that C.numRows() == A.numColumns()
* and B.numColumns() == C.numColumns()
* @return C
*/
Matrix transAmultAdd(double alpha, Matrix B, Matrix C);
/**
* C = A*BT
*
* @param B
* Matrix such that B.numRows() == A.numRows() and
* B.numColumns() == C.numColumns()
* @param C
* Matrix such that C.numRows() == A.numColumns()
* and B.numColumns() == C.numColumns()
* @return C
*/
Matrix transBmult(Matrix B, Matrix C);
/**
* C = alpha*A*BT
*
* @param B
* Matrix such that B.numRows() == A.numRows() and
* B.numColumns() == C.numColumns()
* @param C
* Matrix such that C.numRows() == A.numColumns()
* and B.numColumns() == C.numColumns()
* @return C
*/
Matrix transBmult(double alpha, Matrix B, Matrix C);
/**
* C = A*BT + C
*
* @param B
* Matrix such that B.numRows() == A.numRows() and
* B.numColumns() == C.numColumns()
* @param C
* Matrix such that C.numRows() == A.numColumns()
* and B.numColumns() == C.numColumns()
* @return C
*/
Matrix transBmultAdd(Matrix B, Matrix C);
/**
* C = alpha*A*BT + C
*
* @param B
* Matrix such that B.numRows() == A.numRows() and
* B.numColumns() == C.numColumns()
* @param C
* Matrix such that C.numRows() == A.numColumns()
* and B.numColumns() == C.numColumns()
* @return C
*/
Matrix transBmultAdd(double alpha, Matrix B, Matrix C);
/**
* C = AT*BT
*
* @param B
* Matrix such that B.numColumns() == A.numRows()
* and B.numRows() == C.numColumns()
* @param C
* Matrix such that C.numRows() == A.numColumns()
* and B.numRows() == C.numColumns()
* @return C
*/
Matrix transABmult(Matrix B, Matrix C);
/**
* C = alpha*AT*BT
*
* @param B
* Matrix such that B.numColumns() == A.numRows()
* and B.numRows() == C.numColumns()
* @param C
* Matrix such that C.numRows() == A.numColumns()
* and B.numRows() == C.numColumns()
* @return C
*/
Matrix transABmult(double alpha, Matrix B, Matrix C);
/**
* C = AT*BT + C
*
* @param B
* Matrix such that B.numColumns() == A.numRows()
* and B.numRows() == C.numColumns()
* @param C
* Matrix such that C.numRows() == A.numColumns()
* and B.numRows() == C.numColumns()
* @return C
*/
Matrix transABmultAdd(Matrix B, Matrix C);
/**
* C = alpha*AT*BT + C
*
* @param B
* Matrix such that B.numColumns() == A.numRows()
* and B.numRows() == C.numColumns()
* @param C
* Matrix such that C.numRows() == A.numColumns()
* and B.numRows() == C.numColumns()
* @return C
*/
Matrix transABmultAdd(double alpha, Matrix B, Matrix C);
/**
* X = A\B. Not all matrices support this operation, those that
* do not throw UnsupportedOperationException. Note that it is
* often more efficient to use a matrix decomposition and its associated
* solver
*
* @param B
* Matrix with the same number of rows as A, and the
* same number of columns as X
* @param X
* Matrix with a number of rows equal A.numColumns()
* , and the same number of columns as B
* @return X
* @throws MatrixSingularException
* If the matrix is singular
* @throws MatrixNotSPDException
* If the solver assumes that the matrix is symmetrical,
* positive definite, but that that property does not hold
*/
Matrix solve(Matrix B, Matrix X) throws MatrixSingularException,
MatrixNotSPDException;
/**
* X = AT\B. Not all matrices support this
* operation, those that do not throw
* UnsupportedOperationException. Note that it is often more
* efficient to use a matrix decomposition and its associated transpose
* solver
*
* @param B
* Matrix with a number of rows equal A.numColumns()
* , and the same number of columns as X
* @param X
* Matrix with the same number of rows as A, and the
* same number of columns as B
* @return X
* @throws MatrixSingularException
* If the matrix is singular
* @throws MatrixNotSPDException
* If the solver assumes that the matrix is symmetrical,
* positive definite, but that that property does not hold
*/
Matrix transSolve(Matrix B, Matrix X) throws MatrixSingularException,
MatrixNotSPDException;
/**
* A = C*CT + A. The matrices must be square and of
* the same size
*
* @return A
*/
Matrix rank1(Matrix C);
/**
* A = alpha*C*CT + A. The matrices must be square
* and of the same size
*
* @return A
*/
Matrix rank1(double alpha, Matrix C);
/**
* A = CT*C + A The matrices must be square and of
* the same size
*
* @return A
*/
Matrix transRank1(Matrix C);
/**
* A = alpha*CT*C + A The matrices must be square
* and of the same size
*
* @return A
*/
Matrix transRank1(double alpha, Matrix C);
/**
* A = B*CT + C*BT + A. This matrix must
* be square
*
* @param B
* Matrix with the same number of rows as A and the
* same number of columns as C
* @param C
* Matrix with the same number of rows as A and the
* same number of columns as B
* @return A
*/
Matrix rank2(Matrix B, Matrix C);
/**
* A = alpha*B*CT + alpha*C*BT + A. This
* matrix must be square
*
* @param B
* Matrix with the same number of rows as A and the
* same number of columns as C
* @param C
* Matrix with the same number of rows as A and the
* same number of columns as B
* @return A
*/
Matrix rank2(double alpha, Matrix B, Matrix C);
/**
* A = BT*C + CT*B + A. This matrix must
* be square
*
* @param B
* Matrix with the same number of rows as C and the
* same number of columns as A
* @param C
* Matrix with the same number of rows as B and the
* same number of columns as A
* @return A
*/
Matrix transRank2(Matrix B, Matrix C);
/**
* A = alpha*BT*C + alpha*CT*B + A. This
* matrix must be square
*
* @param B
* Matrix with the same number of rows as C and the
* same number of columns as A
* @param C
* Matrix with the same number of rows as B and the
* same number of columns as A
* @return A
*/
Matrix transRank2(double alpha, Matrix B, Matrix C);
/**
* A = alpha*A
*
* @return A
*/
Matrix scale(double alpha);
/**
* A=B. The matrices must be of the same size
*
* @return A
*/
Matrix set(Matrix B);
/**
* A=alpha*B. The matrices must be of the same size
*
* @return A
*/
Matrix set(double alpha, Matrix B);
/**
* A = B + A. The matrices must be of the same size
*
* @return A
*/
Matrix add(Matrix B);
/**
* A = alpha*B + A. The matrices must be of the same size
*
* @return A
*/
Matrix add(double alpha, Matrix B);
/**
* Transposes the matrix in-place. In most cases, the matrix must be square
* for this to work.
*
* @return This matrix
*/
Matrix transpose();
/**
* Sets the tranpose of this matrix into B. Matrix dimensions
* must be compatible
*
* @param B
* Matrix with as many rows as this matrix has columns, and as
* many columns as this matrix has rows
* @return The matrix B=AT
*/
Matrix transpose(Matrix B);
/**
* Computes the given norm of the matrix
*
* @param type
* The type of norm to compute
*/
double norm(Norm type);
/**
* Supported matrix-norms. Note that Maxvalue is not a proper
* matrix norm
*/
enum Norm {
/**
* Maximum absolute row sum
*/
One,
/**
* The root of sum of the sum of squares
*/
Frobenius,
/**
* Maximum column sum
*/
Infinity,
/**
* Largest entry in absolute value. Not a proper matrix norm
*/
Maxvalue;
/**
* @return the String as required by the netlib libraries to represent
* this norm.
*/
public String netlib() {
// TODO: this is a bit of a hack
// shouldn't need to know about the internals of netlib
if (this == One)
return "1";
else if (this == Infinity)
return "I";
else
throw new IllegalArgumentException(
"Norm must be the 1 or the Infinity norm");
}
}
}