mikera.matrixx.decompose.impl.chol.Cholesky Maven / Gradle / Ivy
Show all versions of vectorz Show documentation
/*
* 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.decompose.impl.chol;
import mikera.matrixx.AMatrix;
import mikera.matrixx.Matrix;
import mikera.matrixx.decompose.ICholeskyResult;
/**
* This is an implementation of Cholesky that processes internal submatrices as blocks. This is
* done to reduce the number of cache issues.
*
* @author Peter Abeles
*/
public class Cholesky extends CholeskyCommon {
private int blockWidth; // how wide the blocks should be
private Matrix B; // row rectangular matrix
private CholeskyHelper chol;
/**
* Default block width
*/
private static final int BLOCK_WIDTH = 60;
private Cholesky() { blockWidth = Cholesky.BLOCK_WIDTH; }
private Cholesky(int blockWidth) { this.blockWidth = blockWidth; }
/**
*
* Computes the Cholesky Decomposition (A = LU) of a matrix, taking
* default block width = 60.
*
*
* If the matrix is not positive definite then this function will return
* null since it can't complete its computations. Not all errors will be
* found. This is an efficient way to check for positive definiteness.
*
* @param mat A symmetric positive definite matrix
* @return A Cholesky Decomposition Result
*/
public static ICholeskyResult decompose(AMatrix mat) {
return decompose(mat, BLOCK_WIDTH);
}
/**
*
* Computes the Cholesky LDU Decomposition (A = LDU) of a matrix.
*
*
* If the matrix is not positive definite then this function will return
* null since it can't complete its computations. Not all errors will be
* found. This is an efficient way to check for positive definiteness.
*
* @param mat A symmetric positive definite matrix
* @param blockWidth The width of a block.
* @return ICholeskyResult if decomposition is successful, null otherwise.
*/
public static ICholeskyResult decompose(AMatrix mat, int blockWidth) {
Cholesky temp = new Cholesky(blockWidth);
return temp._decompose(mat);
}
/**
*
* Performs Choleksy decomposition on the provided matrix.
*
*
*
* If the matrix is not positive definite then this function will return
* null since it can't complete its computations. Not all errors will be
* found. This is an efficient way to check for positive definiteness.
*
* @param mat A symmetric positive definite matrix
* @return CholeskyResult if decomposition is successful, null otherwise
*/
@Override
protected ICholeskyResult _decompose( AMatrix mat ) {
int rc=mat.rowCount();
int cc=mat.columnCount();
if( rc != cc ) {
throw new IllegalArgumentException("Must be a square matrix.");
}
n = mat.rowCount();
this.vv = new double[n];
t=mat.toDoubleArray();
T = Matrix.wrap(rc, cc, t);
if(mat.rowCount() < blockWidth) {
B = Matrix.create(0,0);
}
else {
B = Matrix.create(blockWidth,n);
}
chol = new CholeskyHelper(blockWidth);
return decomposeLower();
}
/**
*
* Performs Choleksy decomposition on the provided matrix.
*
*
*
* If the matrix is not positive definite then this function will return
* null since it can't complete its computations. Not all errors will be
* found.
*
* @return CholeskyResult if decomposition is successful, null otherwise
*/
@Override
protected CholeskyResult decomposeLower() {
if( n < blockWidth)
// B.reshape(0,0, false);
B = Matrix.create(0, 0);
else
// B.reshape(blockWidth,n-blockWidth, false);
B = Matrix.create(blockWidth, n-blockWidth);
int numBlocks = n / blockWidth;
int remainder = n % blockWidth;
if( remainder > 0 ) {
numBlocks++;
}
/*
* In Ejml, DenseMatrix64F has a mutable public field numcols. b_numCols is used here, in place of it.
*/
// B.numCols = n;
int b_numCols = n;
for( int i = 0; i < numBlocks; i++ ) {
// B.numCols -= blockWidth;
b_numCols -= blockWidth;
// if( B.numCols > 0 ) {
if( b_numCols > 0 ) {
// apply cholesky to the current block
if( !chol.decompose(T,(i*blockWidth)* T.columnCount() + i*blockWidth,blockWidth) ) return null;
int indexSrc = i*blockWidth* T.columnCount() + (i+1)*blockWidth;
int indexDst = (i+1)*blockWidth* T.columnCount() + i*blockWidth;
// B = L^(-1) * B
solveL_special(chol.getL().toMatrix().data, T,indexSrc,indexDst,B, b_numCols);
int indexL = (i+1)*blockWidth*n + (i+1)*blockWidth;
// c = c - a^T*a
symmRankTranA_sub(B, T,indexL, b_numCols);
} else {
int width = remainder > 0 ? remainder : blockWidth;
if( !chol.decompose(T,(i*blockWidth)* T.columnCount() + i*blockWidth,width) ) return null;
}
}
// zero the top right corner.
for( int i = 0; i < n; i++ ) {
for( int j = i+1; j < n; j++ ) {
t[i*n+j] = 0.0;
}
}
return new CholeskyResult(T);
}
/*
* Private version of solveL_special that takes an argument, b_numCols which represents B.numCols
* field in DenseMatrix64F in the ejml code
*/
private static void solveL_special( final double L[] ,
final Matrix b_src,
final int indexSrc , final int indexDst ,
final Matrix B, int b_numCols )
{
final double dataSrc[] = b_src.data;
final double b[]= B.data;
final int m = B.rowCount();
final int n = b_numCols;
final int widthL = m;
// for( int j = 0; j < n; j++ ) {
// for( int i = 0; i < widthL; i++ ) {
// double sum = dataSrc[indexSrc+i*b_src.numCols+j];
// for( int k=0; k