org.ejml.alg.dense.linsol.qr.BaseLinearSolverQrp Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of ejml Show documentation
Show all versions of ejml Show documentation
A fast and easy to use dense matrix linear algebra library written in Java.
/*
* Copyright (c) 2009-2012, Peter Abeles. All Rights Reserved.
*
* This file is part of Efficient Java Matrix Library (EJML).
*
* EJML 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 3
* of the License, or (at your option) any later version.
*
* EJML 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 EJML. If not, see
* Base class for QR pivot based pseudo inverse classes. It will return either the
* basic of minimal 2-norm solution. See [1] for details. The minimal 2-norm solution refers to the solution
* 'x' whose 2-norm is the smallest making it unique, not some other error function.
*
*
*
*
* R = [ R12 R12 ] r P^T*x = [ y ] r Q^T*b = [ c ] r
* [ 0 0 ] m-r [ z ] n -r [ d ] m-r
* r n-r
*
* where r is the rank of the matrix and (m,n) is the dimension of the linear system.
*
*
*
*
*
* The solution 'x' is found by solving the system below. The basic solution is found by setting z=0
*
* [ R_11^-1*(c - R12*z) ]
* x = [ z ]
*
*
*
*
* NOTE: The matrix rank is determined using the provided QR decomposition. [1] mentions that this will not always
* work and could cause some problems.
*
*
*
* [1] See page 258-259 in Gene H. Golub and Charles F. Van Loan "Matrix Computations" 3rd Ed, 1996
*
*
* @author Peter Abeles
*/
public abstract class BaseLinearSolverQrp extends LinearSolverAbstract {
QRPDecomposition decomposition;
// if true then only the basic solution will be found
protected boolean norm2Solution;
protected DenseMatrix64F Y = new DenseMatrix64F(1,1);
protected DenseMatrix64F R = new DenseMatrix64F(1,1);
// stores sub-matrices inside the R matrix
protected DenseMatrix64F R11 = new DenseMatrix64F(1,1);
// store an identity matrix for computing the inverse
protected DenseMatrix64F I = new DenseMatrix64F(1,1);
// rank of the system matrix
protected int rank;
protected LinearSolver internalSolver = LinearSolverFactory.leastSquares(1, 1);
// used to compute optimal 2-norm solution
private DenseMatrix64F W = new DenseMatrix64F(1,1);
/**
* Configures internal parameters.
*
* @param decomposition Used to solve the linear system.
* @param norm2Solution If true then the optimal 2-norm solution will be computed for degenerate systems.
*/
protected BaseLinearSolverQrp(QRPDecomposition decomposition,
boolean norm2Solution)
{
this.decomposition = decomposition;
this.norm2Solution = norm2Solution;
if( internalSolver.modifiesA() )
internalSolver = new LinearSolverSafe(internalSolver);
}
@Override
public boolean setA(DenseMatrix64F A) {
_setA(A);
decomposition.setSingularThreshold(CommonOps.elementMaxAbs(A)* UtilEjml.EPS );
if( !decomposition.decompose(A) )
return false;
rank = decomposition.getRank();
R.reshape(numRows,numCols);
decomposition.getR(R,false);
// extract the r11 triangle sub matrix
R11.reshape(rank, rank);
CommonOps.extract(R, 0, rank, 0, rank, R11, 0, 0);
if( norm2Solution && rank < numCols ) {
// extract the R12 sub-matrix
W.reshape(rank,numCols - rank);
CommonOps.extract(R,0,rank,rank,numCols,W,0,0);
// W=inv(R11)*R12
TriangularSolver.solveU(R11.data, 0, R11.numCols, R11.numCols, W.data, 0, W.numCols, W.numCols);
// set the identity matrix in the upper portion
W.reshape(numCols, W.numCols,true);
for( int i = 0; i < numCols-rank; i++ ) {
for( int j = 0; j < numCols-rank; j++ ) {
if( i == j )
W.set(i+rank,j,-1);
else
W.set(i+rank,j,0);
}
}
}
return true;
}
@Override
public double quality() {
return SpecializedOps.qualityTriangular(true, R);
}
/**
*
* Upgrades the basic solution to the optimal 2-norm solution.
*
*
*
* First solves for 'z'
*
* || x_b - P*[ R_11^-1 * R_12 ] * z ||2
* min z || [ - I_{n-r} ] ||
*
*
*
* @param X basic solution, also output solution
*/
protected void upgradeSolution( DenseMatrix64F X ) {
DenseMatrix64F z = Y; // recycle Y
// compute the z which will minimize the 2-norm of X
// because of the identity matrix tacked onto the end 'A' should never be singular
if( !internalSolver.setA(W) )
throw new RuntimeException("This should never happen. Is input NaN?");
z.reshape(numCols-rank,1);
internalSolver.solve(X, z);
// compute X by tweaking the original
CommonOps.multAdd(-1, W, z, X);
}
@Override
public void invert(DenseMatrix64F A_inv) {
if( A_inv.numCols != numRows || A_inv.numRows != numCols )
throw new IllegalArgumentException("Unexpected dimensions for A_inv");
I.reshape(numRows, numRows);
CommonOps.setIdentity(I);
solve(I, A_inv);
}
}