/*
* Copyright (c) 2009-2013, Peter Abeles. All Rights Reserved.
*
* This file is part of Efficient Java Matrix Library (EJML).
*
* 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 mikera.matrixx.solve.impl.qr;
import mikera.matrixx.AMatrix;
import mikera.matrixx.Matrix;
import mikera.matrixx.decompose.impl.qr.HouseholderColQR;
import mikera.matrixx.decompose.impl.qr.QRHelperFunctions;
import mikera.matrixx.decompose.impl.qr.QRResult;
/**
* <p>
* QR decomposition can be used to solve for systems. However, this is not as computationally efficient
* as LU decomposition and costs about 3n<sup>2</sup> flops.
* </p>
* <p>
* It solve for x by first multiplying b by the transpose of Q then solving for the result.
* <br>
* QRx=b<br>
* Rx=Q^T b<br>
* </p>
*
* <p>
* A column major decomposition is used in this solver.
* <p>
*
* @author Peter Abeles
*/
public class QRHouseColSolver {
protected AMatrix A;
protected int numRows;
protected int numCols;
public AMatrix getA() {
return A;
}
protected void _setA(AMatrix A) {
this.A = A;
this.numRows = A.rowCount();
this.numCols = A.columnCount();
}
private HouseholderColQR decomposer;
private Matrix a;
private Matrix temp;
protected int maxRows = -1;
protected int maxCols = -1;
private double[][] QR; // a column major QR matrix
private Matrix R;
private double gammas[];
private QRResult result;
/**
* Creates a linear solver that uses QR decomposition.
*/
public QRHouseColSolver() {
decomposer = new HouseholderColQR(true);
}
public void setMaxSize( int maxRows , int maxCols )
{
this.maxRows = maxRows; this.maxCols = maxCols;
}
/**
* Performs QR decomposition on A
*
* @param A not modified.
*/
public boolean setA(AMatrix A) {
if( A.rowCount() > maxRows || A.columnCount() > maxCols )
setMaxSize(A.rowCount(),A.columnCount());
a = Matrix.create(A.rowCount(),1);
temp = Matrix.create(A.rowCount(),1);
_setA(A);
result = decomposer.decompose(A);
if( result == null )
return false;
gammas = decomposer.getGammas();
QR = decomposer.getQR();
R = result.getR().toMatrix();
return true;
}
public double quality() {
return qualityTriangular(true, R);
}
/**
* Solves for X using the QR decomposition.
*
* @param B A matrix that is n by m. Not modified.
*/
public AMatrix solve(AMatrix B) {
if( B.rowCount() != numRows)
throw new IllegalArgumentException("Unexpected dimensions for B");
Matrix X = Matrix.create(numCols, B.columnCount());
int BnumCols = B.columnCount();
// solve each column one by one
for( int colB = 0; colB < BnumCols; colB++ ) {
// make a copy of this column in the vector
for( int i = 0; i < numRows; i++ ) {
// a.data[i] = B.data[i*BnumCols + colB];
a.data[i] = B.unsafeGet(i, colB);
}
// Solve Qa=b
// a = Q'b
// a = Q_{n-1}...Q_2*Q_1*b
//
// Q_n*b = (I-gamma*u*u^T)*b = b - u*(gamma*U^T*b)
for( int n = 0; n < numCols; n++ ) {
double []u = QR[n];
double vv = u[n];
u[n] = 1;
QRHelperFunctions.rank1UpdateMultR(a, u, gammas[n], 0, n, numRows, temp.data);
u[n] = vv;
}
// solve for Rx = b using the standard upper triangular solver
solveU(R.asDoubleArray(),a.asDoubleArray(),numCols);
// save the results
double[]data = X.asDoubleArray();
for( int i = 0; i < numCols; i++ ) {
// X.data[i*X.columnCount()+colB] = a.data[i];
data[i*X.columnCount()+colB] = a.data[i];
}
}
return X;
}
/**
* <p>
* This is a forward substitution solver for non-singular upper triangular matrices.
* <br>
* b = U<sup>-1</sup>b<br>
* <br>
* where b is a vector, U is an n by n matrix.<br>
* </p>
*
* @param U An n by n non-singular upper triangular matrix. Not modified.
* @param b A vector of length n. Modified.
* @param n The size of the matrices.
*/
private void solveU( double U[] , double []b , int n )
{
for( int i =n-1; i>=0; i-- ) {
double sum = b[i];
int indexU = i*n+i+1;
for( int j = i+1; j <n; j++ ) {
sum -= U[indexU++]* b[j];
}
b[i] = sum/U[i*n+i];
}
}
/**
* Computes the quality of a triangular matrix. In
* this situation the quality is the absolute value of the product of
* each diagonal element divided by the magnitude of the largest diagonal element.
* If all diagonal elements are zero then zero is returned.
*
* @param upper if it is upper triangular or not.
* @param T A matrix. @return product of the diagonal elements.
* @return the quality of the system.
*/
public static double qualityTriangular(boolean upper, AMatrix T)
{
int N = Math.min(T.rowCount(),T.columnCount());
// TODO make faster by just checking the upper triangular portion
double max = T.absCopy().elementMax();
if( max == 0.0d )
return 0.0d;
double quality = 1.0;
for( int i = 0; i < N; i++ ) {
quality *= T.unsafeGet(i,i)/max;
}
return Math.abs(quality);
}
}