EconomyQRDecomposer.java
/*
* Copyright (C) 2012 Alberto Irurueta Carro (alberto@irurueta.com)
*
* Licensed 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 com.irurueta.algebra;
/**
* This decomposer computes economy QR decomposition, which is faster than
* typical QR decomposition.
*/
@SuppressWarnings("DuplicatedCode")
public class EconomyQRDecomposer extends Decomposer {
/**
* Constant defining default round error when determining full rank of
* matrices. This value is zero by default
*/
public static final double DEFAULT_ROUND_ERROR = 0.0;
/**
* Constant defining minimum allowed round error value when determining full
* rank of matrices.
*/
public static final double MIN_ROUND_ERROR = 0.0;
/**
* Internal matrix containing results of decomposition.
*/
private Matrix qr;
/**
* Internal array containing diagonal of R.
*/
private double[] rDiag;
/**
* Constructor of this class.
*/
public EconomyQRDecomposer() {
super();
qr = null;
rDiag = null;
}
/**
* Constructor of this class.
*
* @param inputMatrix Reference to input matrix to be decomposed.
*/
public EconomyQRDecomposer(final Matrix inputMatrix) {
super(inputMatrix);
qr = null;
rDiag = null;
}
/**
* Returns decomposer type corresponding to Economy QR decomposition.
*
* @return Decomposer type.
*/
@Override
public DecomposerType getDecomposerType() {
return DecomposerType.QR_ECONOMY_DECOMPOSITION;
}
/**
* Sets reference to input matrix to be decomposed.
*
* @param inputMatrix Reference to input matrix to be decomposed.
* @throws LockedException Exception thrown if attempting to call this
* method while this instance remains locked.
*/
@Override
public void setInputMatrix(final Matrix inputMatrix) throws LockedException {
super.setInputMatrix(inputMatrix);
qr = null;
rDiag = null;
}
/**
* Returns boolean indicating whether decomposition has been computed and
* results can be retrieved.
* Attempting to retrieve decomposition results when not available, will
* probably raise a NotAvailableException.
*
* @return Boolean indicating whether decomposition has been computed and
* results can be retrieved.
*/
@Override
public boolean isDecompositionAvailable() {
return qr != null;
}
/**
* This method computes QR matrix decomposition, which consists on factoring
* provided input matrix into an orthogonal matrix (Q) and an upper
* triangular matrix (R).
* In other words, if input matrix is A, then: A = Q * R
* Note: During execution of this method, this instance will be locked.
* Note: After execution of this method, QR decomposition will be available
* and operations such as retrieving Q and R matrices or solving systems of
* linear equations will be able to be done. Attempting to call any of such
* operations before calling this method will raise a NotAvailableException
* because they require computation of QR decomposition first.
*
* @throws NotReadyException Exception thrown if attempting to call this
* method when this instance is not ready (i.e. no input matrix has been
* provided).
* @throws LockedException Exception thrown if this decomposer is already
* locked before calling this method. Notice that this method will actually
* lock this instance while it is being executed.
*/
@Override
public void decompose() throws NotReadyException, LockedException {
if (!isReady()) {
throw new NotReadyException();
}
if (isLocked()) {
throw new LockedException();
}
locked = true;
final var rows = inputMatrix.getRows();
final var columns = inputMatrix.getColumns();
qr = new Matrix(inputMatrix);
rDiag = new double[columns];
double nrm;
double s;
// Main loop
for (var k = 0; k < columns; k++) {
// Compute 2-norm of k-th column without under/overflow
nrm = 0.0;
for (var i = k; i < rows; i++) {
nrm = Math.sqrt(nrm * nrm + Math.pow(qr.getElementAt(i, k), 2.0));
}
if (nrm != 0.0) {
// Form k-th Householder vector
if (qr.getElementAt(k, k) < 0.0) {
nrm = -nrm;
}
for (var i = k; i < rows; i++) {
qr.setElementAt(i, k, qr.getElementAt(i, k) / nrm);
}
qr.setElementAt(k, k, qr.getElementAt(k, k) + 1.0);
// Apply transformation to remaining columns
for (var j = k + 1; j < columns; j++) {
s = 0.0;
for (var i = k; i < rows; i++) {
s += qr.getElementAt(i, k) * qr.getElementAt(i, j);
}
s = -s / qr.getElementAt(k, k);
for (var i = k; i < rows; i++) {
qr.setElementAt(i, j, qr.getElementAt(i, j) + s * qr.getElementAt(i, k));
}
}
}
rDiag[k] = -nrm;
}
locked = false;
}
/**
* Returns boolean indicating whether provided input matrix has full rank or
* not.
* Squared matrices having full rank also have determinant different from
* zero.
* Note: Because of rounding errors, testing whether a matrix has full rank
* or not, might obtain unreliable results especially for non-square
* matrices. In such cases matrices usually tend to be considered as full
* rank even when they are not.
*
* @return Boolean indicating whether provided input matrix has full rank or
* not.
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing QR decomposition. To avoid this exception call
* decompose() method first.
* @throws WrongSizeException Exception thrown if provided rounding error is
* lower than minimum allowed value (MIN_ROUND_ERROR).
* @see #decompose()
*/
public boolean isFullRank() throws NotAvailableException, WrongSizeException {
return isFullRank(DEFAULT_ROUND_ERROR);
}
/**
* Returns boolean indicating whether provided input matrix has full rank or
* not.
* Squared matrices having full rank also have determinant different from
* zero.
* Note: Because of rounding errors, testing whether a matrix has full rank
* or not, might obtain unreliable results especially for non-square
* matrices. In such cases matrices usually tend to be considered as full
* rank even when they are not.
*
* @param roundingError Determines the amount of margin given to determine
* whether a matrix has full rank or not due to rounding errors. If not
* provided, by default rounding error is set to zero, but this value can
* be relaxed if needed.
* @return Boolean indicating whether provided input matrix has full rank or
* not.
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing QR decomposition. To avoid this exception call
* decompose() method first.
* @throws WrongSizeException Exception thrown if provided input matrix has
* less rows than columns.
* @throws IllegalArgumentException Exception thrown if provided rounding
* error is lower than minimum allowed value (MIN_ROUND_ERROR).
* @see #decompose()
*/
public boolean isFullRank(final double roundingError) throws NotAvailableException, WrongSizeException {
if (!isDecompositionAvailable()) {
throw new NotAvailableException();
}
if (roundingError < MIN_ROUND_ERROR) {
throw new IllegalArgumentException();
}
final var rows = inputMatrix.getRows();
final var columns = inputMatrix.getColumns();
if (rows < columns) {
throw new WrongSizeException();
}
for (var j = 0; j < columns; j++) {
if (Math.abs(rDiag[j]) < roundingError) {
return false;
}
}
return true;
}
/**
* Computes the Householder vectors and store them in provided matrix.
* Provided matrix will be resized if needed
*
* @param h Matrix where Householder vectors will be stored
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing QR decomposition. To avoid this exception call
* decompose() method first.
* @see #decompose()
*/
public void getH(final Matrix h) throws NotAvailableException {
if (!isDecompositionAvailable()) {
throw new NotAvailableException();
}
final var rows = inputMatrix.getRows();
final var columns = inputMatrix.getColumns();
if (h.getRows() != rows || h.getColumns() != columns) {
try {
h.resize(rows, columns);
} catch (final WrongSizeException ignore) {
// never happens
}
}
for (var i = 0; i < rows; i++) {
for (var j = 0; j < columns; j++) {
if (i >= j) {
h.setElementAt(i, j, qr.getElementAt(i, j));
} else {
h.setElementAt(i, j, 0.0);
}
}
}
}
/**
* Returns the Householder vectors.
*
* @return Lower trapezoidal matrix whose columns define the reflections.
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing QR decomposition. To avoid this exception call
* decompose() method first.
* @see #decompose()
*/
public Matrix getH() throws NotAvailableException {
final var rows = inputMatrix.getRows();
final var columns = inputMatrix.getColumns();
Matrix out = null;
try {
out = new Matrix(rows, columns);
} catch (final WrongSizeException ignore) {
// never happens
}
getH(out);
return out;
}
/**
* Computes upper triangular factor matrix and stores it into provided
* matrix.
* QR decomposition decomposes input matrix into Q (orthogonal matrix) and
* R, which is an upper triangular matrix.
*
* @param r Upper triangular factor matrix
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing QR decomposition. To avoid this exception call
* decompose() method first.
* @see #decompose()
*/
public void getR(final Matrix r) throws NotAvailableException {
if (!isDecompositionAvailable()) {
throw new NotAvailableException();
}
final var columns = inputMatrix.getColumns();
final var rows = inputMatrix.getRows();
if (r.getRows() != columns || r.getColumns() != columns) {
try {
r.resize(columns, columns);
} catch (final WrongSizeException ignore) {
// never happens
}
}
for (var i = 0; i < columns; i++) {
if (i < rows) {
for (var j = 0; j < columns; j++) {
if (i < j) {
r.setElementAt(i, j, qr.getElementAt(i, j));
} else if (i == j) {
r.setElementAt(i, j, rDiag[i]);
} else {
r.setElementAt(i, j, 0.0);
}
}
} else {
for (var j = 0; j < columns; j++) {
r.setElementAt(i, j, 0.0);
}
}
}
}
/**
* Return upper triangular factor matrix.
* QR decomposition decomposes input matrix into Q (orthogonal matrix) and
* R, which is an upper triangular matrix.
*
* @return Upper triangular factor matrix
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing QR decomposition. To avoid this exception call
* decompose() method first.
* @see #decompose()
*/
public Matrix getR() throws NotAvailableException {
final var columns = inputMatrix.getColumns();
Matrix out = null;
try {
out = new Matrix(columns, columns);
} catch (final WrongSizeException ignore) {
// never happens
}
getR(out);
return out;
}
/**
* Computes the economy-sized orthogonal factor matrix and stores it into
* provided matrix.
* QR decomposition decomposes input matrix into Q, which is an orthogonal
* matrix and R (upper triangular matrix).
*
* @param q Orthogonal factor matrix.
* @throws NotAvailableException Exception thrown if attempting to call
* @throws WrongSizeException Exception thrown if provided input matrix has
* less rows than columns.
* @see #decompose()
*/
public void getQ(final Matrix q) throws NotAvailableException, WrongSizeException {
if (!isDecompositionAvailable()) {
throw new NotAvailableException();
}
final var rows = inputMatrix.getRows();
final var columns = inputMatrix.getColumns();
double s;
if (rows < columns) {
throw new WrongSizeException();
}
if (q.getRows() != rows || q.getColumns() != columns) {
try {
q.resize(rows, columns);
} catch (final WrongSizeException ignore) {
// never happens
}
}
for (var k = columns - 1; k >= 0; k--) {
for (var i = 0; i < rows; i++) {
q.setElementAt(i, k, 0.0);
}
q.setElementAt(k, k, 1.0);
for (var j = k; j < columns; j++) {
if (qr.getElementAt(k, k) != 0) {
s = 0.0;
for (var i = k; i < rows; i++) {
s += qr.getElementAt(i, k) * q.getElementAt(i, j);
}
s = -s / qr.getElementAt(k, k);
for (var i = k; i < rows; i++) {
q.setElementAt(i, j, q.getElementAt(i, j) + s * qr.getElementAt(i, k));
}
}
}
}
}
/**
* Return the economy-sized orthogonal factor matrix.
* QR decomposition decomposes input matrix into Q, which is an orthogonal
* matrix and R (upper triangular matrix).
*
* @return Orthogonal factor matrix.
* @throws NotAvailableException Exception thrown if attempting to call
* @throws WrongSizeException Exception thrown if provided input matrix has
* fewer rows than columns.
* @see #decompose()
*/
public Matrix getQ() throws NotAvailableException, WrongSizeException {
final var rows = inputMatrix.getRows();
final var columns = inputMatrix.getColumns();
final var out = new Matrix(rows, columns);
getQ(out);
return out;
}
/**
* Solves a linear system of equations of the following form:
* A * X = B, where A is the input matrix provided for QR decomposition,
* X is the solution to the system of equations, and B is the parameters
* vector/matrix.
* Note: This method can be reused for different b vector/matrices without
* having to recompute QR decomposition on the same input matrix.
* Note: Provided b matrix must have the same number of rows as provided
* input matrix A, otherwise a WrongSizeException will be raised.
* Note: Provided input matrix "A" must have at least as many rows as columns,
* otherwise a WrongSizeException will be raised as well. For input matrices
* having a higher number of rows than columns, the system of equations will
* be overdetermined and the least squares solution will be found.
* Note: If provided input matrix A is rank deficient, a
* RankDeficientMatrixException will be thrown.
* Note: In order to execute this method, a QR decomposition must be
* available, otherwise a NotAvailableException will be raised. In order to
* avoid this exception call decompose() method first.
*
* @param b Parameters matrix that determine a linear system of equations.
* Provided matrix must have the same number of rows as provided input
* matrix for QR decomposition. Besides, each column on parameters matrix
* will represent a new system of equations, whose solution will be returned
* on appropriate column as an output of this method.
* @param result Matrix containing solution of linear system of equations on
* each column for each column of provided matrix of parameters b. Provided
* matrix will be resized if needed
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing QR decomposition. To avoid this exception call
* decompose() method first.
* @throws WrongSizeException Exception thrown if attempting to call this
* method using an input matrix with fewer rows than columns; or if provided
* parameters matrix (b) does not have the same number of rows as input
* matrix being QR decomposed.
* @throws RankDeficientMatrixException Exception thrown if provided input
* matrix to be QR decomposed is rank deficient. In this case linear system
* of equations cannot be solved.
*/
public void solve(final Matrix b, final Matrix result) throws NotAvailableException, WrongSizeException,
RankDeficientMatrixException {
solve(b, DEFAULT_ROUND_ERROR, result);
}
/**
* Solves a linear system of equations of the following form:
* A * X = B, where A is the input matrix provided for QR decomposition,
* X is the solution to the system of equations, and B is the parameters
* vector/matrix.
* Note: This method can be reused for different b vector/matrices without
* having to recompute QR decomposition on the same input matrix.
* Note: Provided b matrix must have the same number of rows as provided
* input matrix A, otherwise a WrongSizeException will be raised.
* Note: Provided input matrix "A" must have at least as many rows as columns,
* otherwise a WrongSizeException will be raised as well. For input matrices
* having a higher number of rows than columns, the system of equations will
* be overdetermined and the least squares solution will be found.
* Note: If provided input matrix A is rank deficient, a
* RankDeficientMatrixException will be thrown.
* Note: In order to execute this method, a QR decomposition must be
* available, otherwise a NotAvailableException will be raised. In order to
* avoid this exception call decompose() method first.
*
* @param b Parameters matrix that determine a linear system of equations.
* Provided matrix must have the same number of rows as provided input
* matrix for QR decomposition. Besides, each column on parameters matrix
* will represent a new system of equations, whose solution will be returned
* on appropriate column as an output of this method.
* @param roundingError threshold to determine whether matrix b has full rank or not.
* By default, this is typically a tiny value close to zero.
* @param result Matrix containing solution of linear system of equations on
* each column for each column of provided matrix of parameters b. Provided
* matrix will be resized if needed
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing QR decomposition. To avoid this exception call
* decompose() method first.
* @throws WrongSizeException Exception thrown if attempting to call this
* method using an input matrix with fewer rows than columns; or if provided
* parameters matrix (b) does not have the same number of rows as input
* matrix being QR decomposed.
* @throws RankDeficientMatrixException Exception thrown if provided input
* matrix to be QR decomposed is rank deficient. In this case linear system
* of equations cannot be solved.
* @throws IllegalArgumentException if provided rounding error is negative.
*/
public void solve(final Matrix b, final double roundingError, final Matrix result)
throws NotAvailableException, WrongSizeException, RankDeficientMatrixException {
if (!isDecompositionAvailable()) {
throw new NotAvailableException();
}
final var rows = inputMatrix.getRows();
final var columns = inputMatrix.getColumns();
final var rowsB = b.getRows();
final var colsB = b.getColumns();
double s;
if (rowsB != rows) {
throw new WrongSizeException();
}
if (roundingError < MIN_ROUND_ERROR) {
throw new IllegalArgumentException();
}
if (rows < columns) {
throw new WrongSizeException();
}
if (!isFullRank(roundingError)) {
throw new RankDeficientMatrixException();
}
// Copy b into X
final var x = new Matrix(b);
// Compute Y = transpose(Q) * B
for (var k = 0; k < columns; k++) {
for (var j = 0; j < colsB; j++) {
s = 0.0;
for (var i = k; i < rows; i++) {
s += qr.getElementAt(i, k) * x.getElementAt(i, j);
}
s = -s / qr.getElementAt(k, k);
for (var i = k; i < rows; i++) {
x.setElementAt(i, j, x.getElementAt(i, j) + s * qr.getElementAt(i, k));
}
}
}
// Solve R * X = Y
for (var k = columns - 1; k >= 0; k--) {
for (var j = 0; j < colsB; j++) {
x.setElementAt(k, j, x.getElementAt(k, j) / rDiag[k]);
}
for (var i = 0; i < k; i++) {
for (var j = 0; j < colsB; j++) {
x.setElementAt(i, j, x.getElementAt(i, j) - x.getElementAt(k, j) * qr.getElementAt(i, k));
}
}
}
// Pick only first columns rows of X in case of overdetermined systems
// (where rows > columns), otherwise rows == columns and we pick them all
if (result.getRows() != columns || result.getColumns() != colsB) {
// resize result
result.resize(columns, colsB);
}
result.setSubmatrix(0, 0, columns - 1, colsB - 1, x,
0, 0, columns - 1,
colsB - 1);
}
/**
* Solves a linear system of equations of the following form:
* A * X = B, where A is the input matrix provided for QR decomposition,
* X is the solution to the system of equations, and B is the parameters
* vector/matrix.
* Note: This method can be reused for different b vector/matrices without
* having to recompute QR decomposition on the same input matrix.
* Note: Provided b matrix must have the same number of rows as provided
* input matrix A, otherwise a WrongSizeException will be raised.
* Note: Provided input matrix "A" must have at least as many rows as columns,
* otherwise a WrongSizeException will be raised as well. For input matrices
* having a higher number of rows than columns, the system of equations will
* be overdetermined and the least squares solution will be found.
* Note: If provided input matrix A is rank deficient, a
* RankDeficientMatrixException will be thrown.
* Note: In order to execute this method, a QR decomposition must be
* available, otherwise a NotAvailableException will be raised. In order to
* avoid this exception call decompose() method first.
*
* @param b Parameters matrix that determine a linear system of equations.
* Provided matrix must have the same number of rows as provided input
* matrix for QR decomposition. Besides, each column on parameters matrix
* will represent a new system of equations, whose solution will be returned
* on appropriate column as an output of this method.
* @return Matrix containing solution of linear system of equations on each
* column for each column of provided matrix of parameters b.
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing QR decomposition. To avoid this exception call
* decompose() method first.
* @throws WrongSizeException Exception thrown if attempting to call this
* method using an input matrix with fewer rows than columns; or if provided
* parameters matrix (b) does not have the same number of rows as input
* matrix being QR decomposed.
* @throws RankDeficientMatrixException Exception thrown if provided input
* matrix to be QR decomposed is rank deficient. In this case linear system
* of equations cannot be solved.
*/
public Matrix solve(final Matrix b) throws NotAvailableException, WrongSizeException, RankDeficientMatrixException {
return solve(b, DEFAULT_ROUND_ERROR);
}
/**
* Solves a linear system of equations of the following form:
* A * X = B, where A is the input matrix provided for QR decomposition,
* X is the solution to the system of equations, and B is the parameters
* vector/matrix.
* Note: This method can be reused for different b vector/matrices without
* having to recompute QR decomposition on the same input matrix.
* Note: Provided b matrix must have the same number of rows as provided
* input matrix A, otherwise a WrongSizeException will be raised.
* Note: Provided input matrix "A" must have at least as many rows as columns,
* otherwise a WrongSizeException will be raised as well. For input matrices
* having a higher number of rows than columns, the system of equations will
* be overdetermined and the least squares solution will be found.
* Note: If provided input matrix A is rank deficient, a
* RankDeficientMatrixException will be thrown.
* Note: In order to execute this method, a QR decomposition must be
* available, otherwise a NotAvailableException will be raised. In order to
* avoid this exception call decompose() method first.
*
* @param b Parameters matrix that determine a linear system of equations.
* Provided matrix must have the same number of rows as provided input
* matrix for QR decomposition. Besides, each column on parameters matrix
* will represent a new system of equations, whose solution will be returned
* on appropriate column as an output of this method.
* @param roundingError Determines the amount of margin given to determine
* whether a matrix has full rank or not due to rounding errors. If not
* provided, by default rounding error is set to zero, but this value can be
* relaxed if needed.
* @return Matrix containing solution of linear system of equations on each
* column for each column of provided matrix of parameters b.
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing QR decomposition. To avoid this exception call
* decompose() method first.
* @throws WrongSizeException Exception thrown if attempting to call this
* method using an input matrix with less rows than columns; or if provided
* parameters matrix (b) does not have the same number of rows as input
* matrix being QR decomposed.
* @throws RankDeficientMatrixException Exception thrown if provided input
* matrix to be QR decomposed is rank deficient. In this case linear system
* of equations cannot be solved.
* @throws IllegalArgumentException Exception thrown if provided rounding
* error is lower than minimum allowed value (MIN_ROUND_ERROR)
* @see #decompose()
*/
public Matrix solve(final Matrix b, final double roundingError) throws NotAvailableException, WrongSizeException,
RankDeficientMatrixException {
final var columns = inputMatrix.getColumns();
final var colsB = b.getColumns();
final var out = new Matrix(columns, colsB);
solve(b, roundingError, out);
return out;
}
}