package edu.stanford.rsl.conrad.numerics;
/**
* Implements a QR decomposition for arbitrary matrices.
*
* This functor object allows to decompose any matrix \f$ \mathbf{A} \f$ into matrices
* \f$ \mathbf{Q} \f$ and \f$ \mathbf{R} \f$ such that
* \f[
* \mathbf{A} = \mathbf{Q} \mathbf{R}.
* \f]
* Use it for solving systems of linear equations and for minimization problems
* (with more equations than variables). The QR decomposition is not suited for
* computing the set of solutions for an under-determined system of equations
* (although it may be used to compute one of those solutions).
*
* Assuming \f$ \mathbf A \f$ is a \f$ m \times n \f$ matrix, the decomposition is
* performed such that the matrix \f$ \mathbf Q \f$ has dimensions \f$ m \times m \f$ and
* is orthogonal and the matrix \f$ \mathbf{R} \f$ has dimensions \f$ m \times n \f$ and
* has upper triangular form (\f$ r_{i,j} = 0 \f$ for \f$ i>j \f$).
* Internally, only a compressed version of the factors \f$ \mathbf{Q} \f$ and
* \f$ \mathbf{R} \f$ is stored which together needs approximately as much space as the
* original matrix \f$ \mathbf{A} \f$.
*
* Usage:
* \code
* Nude::Decomposition::QR<> qr;
* qr(A);
* qr.solve(b);
* \endcode
*
* \sa RQ
* \sa Deuflhard/Hohmann: Numerische Mathematik I
* \sa The development of this code started out with the implementation of TNT/JAMA at
* http://math.nist.gov/tnt/. However, it was completely revised, optimized, and commented.
*
* \author Andreas Keil
*
* \todo Handle singularity / do pivoting when decomposing.
* \todo Improve method isFullRank() to check against an epsilon.
*/
public class DecompositionQR {
/**
* Constructor performing the actual decomposition of a matrix {@latex.inline $\\mathbf{A}$}
* and storing the result in an internal format.
*
* Decomposition is performed as soon as a Matrix is given to this
* constructor and the result is stored internally. Afterwards, the other
* other members ({@link #solve(SimpleVector b)}, {@link #solve(SimpleMatrix B)},
* {@link #getQ()}, and {@link #getR()}) can be used multiple times without having
* to recompute the decomposition.
*
* @param A The Matrix to be decomposed.
*/
public DecompositionQR(final SimpleMatrix A) {
// get dimensions
final int m = A.getRows();
final int n = A.getCols();
final int min_mn = Math.min(m, n); // number of Householder reflections used to construct an upper triangular R from A
// initialize internal variables
this.QR = new SimpleMatrix(A);
this.Rdiag = new SimpleVector(min_mn);
// main loop over columns of A which are to be reduced to upper triangular form
for (int k = 0; k < min_mn; ++k) {
// compute 2-norm of k-th column without under/overflow using C++'s _hypot for computing the sqrt(a^2 + b^2)
double alpha = 0.0;
for (int i = k; i < m; ++i) alpha += this.QR.getElement(i, k) * this.QR.getElement(i, k);
alpha = Math.sqrt(alpha);
if (alpha != 0.0) { //TODO: error handling or permutation for singular case when alpha == 0.0?
// choose the sign for alpha to prevent loss of significance
if (this.QR.getElement(k, k) > 0.0) alpha = -alpha;
// form k-th Householder vector v as (y/alpha - e_1), so that Qy = alpha e_1
for (int i = k; i < m; ++i) this.QR.divideElementBy(i, k, alpha);
this.QR.subtractFromElement(k, k, 1.0);
// apply transformation to remaining columns x, where Qx = x + (v'*x / v_1) * v
for (int j = k + 1; j < n; ++j) {
// compute v'*x
double s = 0.0;
for (int i = k; i < m; ++i) s += this.QR.getElement(i, k) * this.QR.getElement(i, j);
// compute v'*x / v_1
s /= this.QR.getElement(k, k);
// compute the final result Qx = x + (v'*x / v_1) * v
for (int i = k; i < m; ++i) this.QR.addToElement(i, j, s*this.QR.getElement(i, k));
}
}
// store diagonal entry of R
this.Rdiag.setElementValue(k, alpha);
}
}
/**
* Specifies whether the input Matrix {@latex.inline $A$} has full rank.
*
* @return Whether the input Matrix has full rank ({@latex.inline $\\min\\{m, n\\}$}).
*/
public boolean isFullRank() {
final int min_mn = this.Rdiag.getLen();
for (int l = 0; l < min_mn; ++l)
if (this.Rdiag.getElement(l) == 0.0) return false;
return true;
}
/**
* Computes solution Vector {@latex.inline $\\mathbf x$} for the right-hand-side {@latex.inline $\\mathbf b$}.
*
* Depending on the size of {@latex.inline $\\mathbf A \\in \\mathbb{R}^{m \\times n}$}, the problem task
* can either be the solution of a system of equations with a unique solution or a
* minimization task:
* task:
* {@latex.ilb %preamble{\\usepackage{amsmath}} \\begin{align}
* \\mathbf A \\mathbf x = \\mathbf b &: m = n \\quad\\text{(unique solution)} \\\\
* \\min_{\\mathbf x}\\| \\mathbf A \\mathbf x - \\mathbf b \\|_2 &: m > n \\quad\\text{(minimization)} \\\\
* \\mathbf A \\mathbf x = \\mathbf b &: m < n \\quad\\text{(one of the infinite set of solutions)}
* \\end{align} }
* For {@latex.inline $m < n$} the QR decomposition computes one of the infinite number of
* solutions to the under-determined system of equations. Better use the RQ
* decomposition in this case.
*
* <p><em>Remark:</em> If you want to further improve the accuracy of your solution in case (1),
* perform a correction step in the following manner:<br>
* {@code
* x = qr.solve(b);
* x += qr.solve(b - A*x);
* }
*
* @param b The right-hand-side Vector.
* @return The solution Vector {@latex.inline $\\mathbf x$}.
*
* @see #solve(SimpleMatrix)
*/
public SimpleVector solve(final SimpleVector b) {
final SimpleMatrix B = new SimpleMatrix(b.getLen(), 1);
B.setColValue(0, b);
final SimpleMatrix X = this.solve(B);
return X.getCol(0);
}
/**
* Computes solution Matrix {@latex.inline $\\mathbf X$} for the right-hand-side {@latex.inline $\\mathbf B$}.
*
* This method does the same as {@link #solve(SimpleVector)} but for multiple
* right-hand-side vectors, given as a Matrix.
*
* @param B The right-hand-side Matrix for which a solution is computed column-wise.
* @return The solution Matrix {@latex.inline $\\mathbf X$}.
*
* @see #solve(SimpleVector)
*/
public SimpleMatrix solve(final SimpleMatrix B) {
final int m = this.QR.getRows();
final int n = this.QR.getCols();
final int min_mn = Math.min(m, n); // number of Householder reflections used to construct an upper triangular R from A
final int nb = B.getCols(); // column number of right-hand-side (number of RHSs to compute a solution for)
// check for conformant arrays
if (B.getRows() != m) throw new IllegalArgumentException("Number of rows of A and B do not conform!");
// check whether we can actually use the QR decomposition for solving the given problem
//if (n > m) ; //TODO Maybe issue a warning here, since for an under-determined system, QR can compute a solution candidate, but not the one with the smallest norm and not the whole solution space (as can RQ)
// check for rank deficiency
if (!isFullRank()) throw new IllegalArgumentException("Matrix does not have full rank!");
// make internal copy of B for further inplace computations
SimpleMatrix X = B.clone();
// compute X = Q^T * B = Qmin_mn * ... * Q1 * B (inplace in X)
for (int k = 0; k < min_mn; ++k) {
// iterate through columns of B
for (int j = 0; j < nb; ++j) {
// compute v'*x
double s = 0.0;
for (int i = k; i < m; ++i) s += this.QR.getElement(i, k) * X.getElement(i, j);
// compute v'*x / v_1
s /= this.QR.getElement(k, k);
// compute the final result Qx = x + (v'*x / v_1) * v
for (int i = k; i < m; ++i) X.addToElement(i, j, s*this.QR.getElement(i,k));
}
}
// solve R*X' = X (inplace in X)
// iterate through columns of X to compute columns of X'
for (int j = 0; j < nb; ++j) {
// iterate through non-zero rows of R to compute non-zero rows of X'
for (int k = min_mn-1; k >= 0; --k) {
for (int i = k+1; i < min_mn; ++i) X.subtractFromElement(k, j, this.QR.getElement(k, i) * X.getElement(i, j));
X.divideElementBy(k, j, this.Rdiag.getElement(k));
}
}
// return solution (m == n), approximation (m > n), or a representative among the set of solutions (m < n)
if (m > n) {
X = X.getSubMatrix(0, 0, n, nb);
} else if (m < n) {
SimpleMatrix Xnew = new SimpleMatrix(n, nb);
Xnew.setSubMatrixValue(0, 0, X);
X = Xnew;
for (int i = m+1; i < n; ++i)
for (int j = 0; j < nb; ++j)
X.setElementValue(i, j, 0.0);
}
return X;
}
/**
* Compute R from the internal storages QR and Rdiag.
*
* R has the same dimensions as the original matrix.
* An economy-sized version would only return the non-zero submatrix.
*/
public SimpleMatrix getR() {
final int m = this.QR.getRows();
final int n = this.QR.getCols();
final int min_mn = Math.min(m, n); // number of Householder reflections used to construct an upper triangular R from A
// construct R as zero matrix
final SimpleMatrix R = new SimpleMatrix(m, n);
// fill non-zero positions of R with Rdiag and upper-right part of QR
for (int i = 0; i < min_mn; ++i) {
R.setElementValue(i, i, this.Rdiag.getElement(i));
for (int j = i+1; j < n; ++j) R.setElementValue(i, j, this.QR.getElement(i,j));
}
// return result
return R;
}
/**
* Compute Q from the internal storage QR.
*
* Q is a square matrix with both dimensions equaling the row number of the original matrix.
* An economy-sized version would only return the columns of Q corresponding to the economy-sized version of R.
*/
public SimpleMatrix getQ() {
final int m = this.QR.getRows();
final int n = this.QR.getCols();
final int min_mn = Math.min(m, n); // number of Householder reflections used to construct an upper triangular R from A
// start with an identity matrix
final SimpleMatrix Q = new SimpleMatrix(m, m);
Q.identity();
// iterate through Householder reflections to compute Q = Q1 * ... * Qmin_mn * I
for (int k = min_mn - 1; k >= 0; --k) {
if (this.QR.getElement(k, k) != 0.0) { // a zero entry in this position means we didn't have to reflect anything since we encountered a singularity => Q_k = I
// iterate through columns x of Q
for (int j = k; j < m; ++j) {
// compute v'*x
double s = 0.0;
for (int i = k; i < m; ++i) s += this.QR.getElement(i, k) * Q.getElement(i, j);
// compute v'*x / v_1
s /= this.QR.getElement(k, k);
// compute the final result Qx = x + (v'*x / v_1) * v
for (int i = k; i < m; ++i) Q.addToElement(i, j, s*this.QR.getElement(i,k));
}
}
}
// return result
return Q;
}
private final SimpleMatrix QR;
private final SimpleVector Rdiag;
}
/*
* Copyright (C) 2010-2014 Andreas Keil
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/