com.joptimizer.algebra.CholeskySparseFactorization Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of joptimizer Show documentation
Show all versions of joptimizer Show documentation
JOptimizer is a java library for solving general convex optimization problems, like for example LP, QP, QCQP, SOCP, SDP, GP and many more.
/*
* Copyright 2011-2017 joptimizer.com
*
* This work is licensed under the Creative Commons Attribution-NoDerivatives 4.0
* International License. To view a copy of this license, visit
*
* http://creativecommons.org/licenses/by-nd/4.0/
*
* or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.
*/
package com.joptimizer.algebra;
import org.apache.commons.lang3.ArrayUtils;
import org.apache.commons.logging.Log;
import org.apache.commons.logging.LogFactory;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import com.joptimizer.exception.JOptimizerException;
import com.joptimizer.util.ColtUtils;
import com.joptimizer.util.Utils;
import cern.colt.function.tdouble.IntIntDoubleFunction;
import cern.colt.matrix.tdouble.DoubleFactory1D;
import cern.colt.matrix.tdouble.DoubleFactory2D;
import cern.colt.matrix.tdouble.DoubleMatrix1D;
import cern.colt.matrix.tdouble.DoubleMatrix2D;
import cern.colt.matrix.tdouble.impl.SparseDoubleMatrix2D;
/**
* Cholesky factorization and inverse for symmetric and positive sparse matrix
* Q = L.L[T], with L left-lower triangular.
*
* NOTE: this class allows to factorize a matrix the is filled only in its subdiagonal lower left part.
*
* @see "Computing the Cholesky Factorization of Sparse Matrices" online free available pdf
* @author alberto trivellato ([email protected])
*/
public class CholeskySparseFactorization {
private int dim;
private SparseDoubleMatrix2D Q;
private MatrixRescaler rescaler = null;
private DoubleMatrix1D U;
protected DoubleFactory2D F2 = DoubleFactory2D.dense;
protected DoubleFactory1D F1 = DoubleFactory1D.dense;
private double[][] LcolumnsValues = null;
private double[] Svalues = null;
private DoubleMatrix2D L;
private DoubleMatrix2D LT;
private static Log log = LogFactory.getLog(CholeskySparseFactorization.class.getName());
public long factorizeTime = 0l;
public long foreachTime = 0l;
public CholeskySparseFactorization(SparseDoubleMatrix2D Q) {
this(Q, null);
}
/**
* @param Q sparse symmetric positive definite matrix. Only its subdiagonal lower left part is relevant.
*/
public CholeskySparseFactorization(SparseDoubleMatrix2D Q, MatrixRescaler rescaler) {
//ColtUtils.dumpSparseMatrix(Q);
//log.debug(org.apache.commons.lang3.ArrayUtils.toString(Q.toArray()));
this.dim = Q.rows();
this.Q = Q;
this.rescaler = rescaler;
}
/**
* Q = L.L[T], L lower-left triangular. Construction of the matrix L.
*/
public void factorize() throws JOptimizerException {
long t0 = System.currentTimeMillis();
this.LcolumnsValues = new double[dim][];
if(this.rescaler != null){
double[] cn_00_original = null;
double[] cn_2_original = null;
double[] cn_00_scaled = null;
double[] cn_2_scaled = null;
if(log.isDebugEnabled()){
cn_00_original = ColtUtils.getConditionNumberRange(new Array2DRowRealMatrix(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()), Integer.MAX_VALUE);
log.debug("cn_00_original Q before scaling: " + ArrayUtils.toString(cn_00_original));
cn_2_original = ColtUtils.getConditionNumberRange(new Array2DRowRealMatrix(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()), 2);
log.debug("cn_2_original Q before scaling : " + ArrayUtils.toString(cn_2_original));
}
//scaling the Q matrix, we have:
//Q1 = U.Q.U[T] = U.L.L[T].U[T] = (U.L).(U.L)[T]
//and because U is diagonal it preserves the triangular form of U.L, so
//Q1 = U.Q.U[T] = L1.L1[T] is the new Cholesky decomposition
DoubleMatrix1D Uv = rescaler.getMatrixScalingFactorsSymm(Q);
if(log.isDebugEnabled()){
boolean checkOK = rescaler.checkScaling(ColtUtils.fillSubdiagonalSymmetricMatrix(Q), Uv, Uv);
if(!checkOK){
log.warn("Scaling failed (checkScaling = false)");
}
}
this.U = Uv;
this.Q = (SparseDoubleMatrix2D)ColtUtils.diagonalMatrixMult(Uv, Q, Uv);
if(log.isDebugEnabled()){
cn_00_scaled = ColtUtils.getConditionNumberRange(new Array2DRowRealMatrix(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()), Integer.MAX_VALUE);
log.debug("cn_00_scaled Q after scaling : " + ArrayUtils.toString(cn_00_scaled));
cn_2_scaled = ColtUtils.getConditionNumberRange(new Array2DRowRealMatrix(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()), 2);
log.debug("cn_2_scaled Q after scaling : " + ArrayUtils.toString(cn_2_scaled));
if(cn_00_original[0] < cn_00_scaled[0] || cn_2_original[0] < cn_2_scaled[0]){
//log.info("Q: " + ArrayUtils.toString(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()));
log.warn("Problematic scaling");
//throw new RuntimeException("Scaling failed");
}
}
}
final int[] currentColumnIndexHolder = new int[] { -1 };
IntIntDoubleFunction myFunct = new IntIntDoubleFunction() {
public double apply(int i, int s, double pis) {
int step = currentColumnIndexHolder[0];
//log.debug("i:" + i + ", step:" + step + ": " + pis);
if (i >= step) {
// at this point, step is the row index and j is the column index,
// but we take into account only the lower left subdiagonal part of Q (that is symmetric)
//log.debug("i:" + i + ", step:" + step + ": " + pis);
Svalues[i - step] = pis;
}
return pis;
}
};
//view Q column by column
for (int step = 0; step < dim; step++) {
//log.debug("column:" + step);
// reset and fill initial values
Svalues = new double[dim - step];
DoubleMatrix2D P = Q.viewPart(0, step, dim, 1);
currentColumnIndexHolder[0] = step;
P.forEachNonZero(myFunct);
doStep(step);
}
this.factorizeTime += (System.currentTimeMillis() - t0);
}
private void doStep(int step) throws JOptimizerException {
// subtract L[j:n,1:j-n]L[j,1:j-1][T] from S (j=step+1)
for (int k = 0; k < step; k++) {
int j = step;//just for easy reading
double[] LcolumnsValuesK = LcolumnsValues[k];
double LJK = LcolumnsValuesK[j - k];// array of variable dimension (dim-colIndex)
if (Double.compare(LJK, 0.) == 0) {
continue;
}
// subtract L[j:n,k]*L[k,j][T]=L[j,k]*L[j:n,k]
for (int i = j; i < dim; i++) {
double LIK = LcolumnsValuesK[i - k];// array of variable dimension (dim-colIndex)
double LJKLIK = LJK * LIK;
Svalues[i - step] -= LJKLIK;
}
}
// compute L[j:n,j]
if (!(Svalues[0] > Utils.getDoubleMachineEpsilon())) {
throw new JOptimizerException("not positive definite matrix");
}
double evStep = Math.sqrt(Svalues[0]);
double[] LcolumnsValuesS = new double[dim - step];// max length of SValues
LcolumnsValuesS[0] = evStep;
for (int i = 1; i < Svalues.length; i++) {
double SvaluesR = Svalues[i];
if (Double.compare(SvaluesR, 0.) != 0) {
SvaluesR = SvaluesR / evStep;
LcolumnsValuesS[i] = SvaluesR;
}
}
LcolumnsValues[step] = LcolumnsValuesS;
if (log.isDebugEnabled()) {
log.debug("step " + step + " situation:");
log.debug("LcolumnsValues: " + ArrayUtils.toString(LcolumnsValues));
log.debug("Svalues: " + ArrayUtils.toString(Svalues));
}
}
/**
* Solves Q.x = b
*/
public DoubleMatrix1D solve(DoubleMatrix1D b) {
if (b.size() != dim) {
log.error("wrong dimension of vector b: expected " + dim + ", actual " + b.size());
throw new RuntimeException("wrong dimension of vector b: expected " + dim + ", actual " + b.size());
}
// with scaling, we must solve U.Q.U.z = U.b, after that we have x = U.z
if (this.rescaler != null) {
// b = ALG.mult(this.U, b);
b = ColtUtils.diagonalMatrixMult(this.U, b);
}
final double[] y = new double[dim];// copy
System.arraycopy(b.toArray(), 0, y, 0, dim);
// Solve L.y = b
for (int j = 0; j < dim; j++) {
final double[] LTJ = LcolumnsValues[j];
y[j] /= LTJ[0];// the diagonal of the matrix L
final double yJ = y[j];
for (int i = j + 1; i < dim; i++) {
y[i] -= yJ * LTJ[i - j];
}
}
// Solve L[T].x = y
final DoubleMatrix1D x = F1.make(dim);
for (int i = dim - 1; i > -1; i--) {
final double[] LTI = LcolumnsValues[i];
double sum = 0;
for (int j = dim - 1; j > i; j--) {
sum += LTI[j - i] * x.getQuick(j);
}
x.setQuick(i, (y[i] - sum) / LTI[0]);
}
if (this.rescaler != null) {
// return ALG.mult(this.U, x);
return ColtUtils.diagonalMatrixMult(this.U, x);
} else {
return x;
}
}
/**
* Solves Q.X = B
*/
public DoubleMatrix2D solve(DoubleMatrix2D B) {
if (B.rows() != dim) {
log.error("wrong dimension of vector b: expected " + dim +", actual " + B.rows());
throw new RuntimeException("wrong dimension of vector b: expected " + dim +", actual " + B.rows());
}
// with scaling, we must solve U.Q.U.z = U.b, after that we have x = U.z
if (this.rescaler != null) {
// B = ALG.mult(this.U, B);
B = ColtUtils.diagonalMatrixMult(this.U, B);
}
int nOfColumns = B.columns();
// copy
final double[][] Y = B.copy().toArray();
// Solve LY = B (same as L.Yc = Bc for every column of Y and B)
for (int j = 0; j < dim; j++) {
final double[] LTJ = LcolumnsValues[j];
for (int col = 0; col < nOfColumns; col++) {
Y[j][col] /= LTJ[0];// the diagonal of the matrix L
final double YJCol = Y[j][col];
if(Double.compare(YJCol, 0.)!=0){
for (int i = j + 1; i < dim; i++) {
Y[i][col] -= YJCol * LTJ[i - j];
}
}
}
}
// Solve L[T].X = Y (same as L[T].Xc = Yc for every column of X and Y)
final DoubleMatrix2D X = F2.make(dim, nOfColumns);
for (int i = dim - 1; i > -1; i--) {
final double[] LTI = LcolumnsValues[i];
double[] sum = new double[nOfColumns];
for (int col = 0; col < nOfColumns; col++) {
for (int j = dim - 1; j > i; j--) {
sum[col] += LTI[j - i] * X.getQuick(j, col);
}
X.setQuick(i, col, (Y[i][col] - sum[col]) / LTI[0]);
}
}
if (this.rescaler != null) {
// return ALG.mult(this.U, X);
return ColtUtils.diagonalMatrixMult(this.U, X);
} else {
return X;
}
}
public DoubleMatrix2D getL() {
if(this.L == null){
this.L = getLT().viewDice();
}
return this.L;
}
//@TODO: check this
public DoubleMatrix2D getLT() {
if (this.LT == null) {
double[][] myLT = new double[dim][];
for (int i = 0; i < dim; i++) {
double[] LTI = new double[dim];
double[] valuesI = LcolumnsValues[i];
for (int j = i; j < dim; j++) {
LTI[j] = valuesI[j - i];
}
myLT[i] = LTI;
}
if (this.rescaler != null) {
//Q = UInv.Q1.UInv[T] = UInv.L1.L1[T].UInv[T] = (UInv.L1).(UInv.L1)[T]
//so
//L = UInv.L1
DoubleMatrix1D UInv = F1.make(dim);
for(int i=0; i