Added a "rectangular" Cholesky decomposition for positive semidefinite matrices.
JIRA: MATH-541 git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1096496 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
6ac4c4d226
commit
bb49301290
|
|
@ -0,0 +1,66 @@
|
|||
/*
|
||||
* Licensed to the Apache Software Foundation (ASF) under one or more
|
||||
* contributor license agreements. See the NOTICE file distributed with
|
||||
* this work for additional information regarding copyright ownership.
|
||||
* The ASF licenses this file to You 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 org.apache.commons.math.linear;
|
||||
|
||||
import org.apache.commons.math.random.CorrelatedRandomVectorGenerator;
|
||||
|
||||
|
||||
/**
|
||||
* An interface to classes that implement an algorithm to calculate a
|
||||
* rectangular variation of Cholesky decomposition of a real symmetric
|
||||
* positive semidefinite matrix.
|
||||
* <p>The rectangular Cholesky decomposition of a real symmetric positive
|
||||
* semidefinite matrix A consists of a rectangular matrix B with the same
|
||||
* number of rows such that: A is almost equal to BB<sup>T</sup>, depending
|
||||
* on a user-defined tolerance. In a sense, this is the square root of A.</p>
|
||||
* <p>The difference with respect to the regular {@link CholeskyDecomposition}
|
||||
* is that rows/columns may be permuted (hence the rectangular shape instead
|
||||
* of the traditional triangular shape) and there is a threshold to ignore
|
||||
* small diagonal elements. This is used for example to generate {@link
|
||||
* CorrelatedRandomVectorGenerator correlated random n-dimensions vectors}
|
||||
* in a p-dimension subspace (p < n). In other words, it allows generating
|
||||
* random vectors from a covariance matrix that is only positive semidefinite,
|
||||
* and not positive definite.</p>
|
||||
* <p>Rectangular Cholesky decomposition is <em>not</em> suited for solving
|
||||
* linear systems, so it does not provide any {@link DecompositionSolver
|
||||
* decomposition solver}.</p>
|
||||
* @see CholeskyDecomposition
|
||||
* @see CorrelatedRandomVectorGenerator
|
||||
* @version $Revision$ $Date$
|
||||
* @since 3.0
|
||||
*/
|
||||
public interface RectangularCholeskyDecomposition {
|
||||
|
||||
/** Get the root of the covariance matrix.
|
||||
* The root is the rectangular matrix <code>B</code> such that
|
||||
* the covariance matrix is equal to <code>B.B<sup>T</sup></code>
|
||||
* @return root of the square matrix
|
||||
* @see #getRank()
|
||||
*/
|
||||
RealMatrix getRootMatrix();
|
||||
|
||||
/** Get the rank of the symmetric positive semidefinite matrix.
|
||||
* The r is the number of independent rows in the symmetric positive semidefinite
|
||||
* matrix, it is also the number of columns of the rectangular
|
||||
* matrix of the decomposition.
|
||||
* @return r of the square matrix.
|
||||
* @see #getRootMatrix()
|
||||
*/
|
||||
int getRank();
|
||||
|
||||
}
|
||||
|
|
@ -0,0 +1,152 @@
|
|||
/*
|
||||
* Licensed to the Apache Software Foundation (ASF) under one or more
|
||||
* contributor license agreements. See the NOTICE file distributed with
|
||||
* this work for additional information regarding copyright ownership.
|
||||
* The ASF licenses this file to You 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 org.apache.commons.math.linear;
|
||||
|
||||
import org.apache.commons.math.util.FastMath;
|
||||
|
||||
/**
|
||||
* Calculates the rectangular Cholesky decomposition of a matrix.
|
||||
* <p>The rectangular Cholesky decomposition of a real symmetric positive
|
||||
* semidefinite matrix A consists of a rectangular matrix B with the same
|
||||
* number of rows such that: A is almost equal to BB<sup>T</sup>, depending
|
||||
* on a user-defined tolerance. In a sense, this is the square root of A.</p>
|
||||
*
|
||||
* @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a>
|
||||
* @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a>
|
||||
* @version $Revision$ $Date$
|
||||
* @since 2.0
|
||||
*/
|
||||
public class RectangularCholeskyDecompositionImpl implements RectangularCholeskyDecomposition {
|
||||
|
||||
/** Permutated Cholesky root of the symmetric positive semidefinite matrix. */
|
||||
private final RealMatrix root;
|
||||
|
||||
/** Rank of the symmetric positive semidefinite matrix. */
|
||||
private int rank;
|
||||
|
||||
/**
|
||||
* Decompose a symmetric positive semidefinite matrix.
|
||||
*
|
||||
* @param matrix Symmetric positive semidefinite matrix.
|
||||
* @param small Diagonal elements threshold under which column are
|
||||
* considered to be dependent on previous ones and are discarded.
|
||||
* @exception NonPositiveDefiniteMatrixException if the matrix is not
|
||||
* positive semidefinite.
|
||||
*/
|
||||
public RectangularCholeskyDecompositionImpl(RealMatrix matrix, double small)
|
||||
throws NonPositiveDefiniteMatrixException {
|
||||
|
||||
int order = matrix.getRowDimension();
|
||||
double[][] c = matrix.getData();
|
||||
double[][] b = new double[order][order];
|
||||
|
||||
int[] swap = new int[order];
|
||||
int[] index = new int[order];
|
||||
for (int i = 0; i < order; ++i) {
|
||||
index[i] = i;
|
||||
}
|
||||
|
||||
int r = 0;
|
||||
for (boolean loop = true; loop;) {
|
||||
|
||||
// find maximal diagonal element
|
||||
swap[r] = r;
|
||||
for (int i = r + 1; i < order; ++i) {
|
||||
int ii = index[i];
|
||||
int isi = index[swap[i]];
|
||||
if (c[ii][ii] > c[isi][isi]) {
|
||||
swap[r] = i;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// swap elements
|
||||
if (swap[r] != r) {
|
||||
int tmp = index[r];
|
||||
index[r] = index[swap[r]];
|
||||
index[swap[r]] = tmp;
|
||||
}
|
||||
|
||||
// check diagonal element
|
||||
int ir = index[r];
|
||||
if (c[ir][ir] < small) {
|
||||
|
||||
if (r == 0) {
|
||||
throw new NonPositiveDefiniteMatrixException(ir, small);
|
||||
}
|
||||
|
||||
// check remaining diagonal elements
|
||||
for (int i = r; i < order; ++i) {
|
||||
if (c[index[i]][index[i]] < -small) {
|
||||
// there is at least one sufficiently negative diagonal element,
|
||||
// the symmetric positive semidefinite matrix is wrong
|
||||
throw new NonPositiveDefiniteMatrixException(i, small);
|
||||
}
|
||||
}
|
||||
|
||||
// all remaining diagonal elements are close to zero, we consider we have
|
||||
// found the rank of the symmetric positive semidefinite matrix
|
||||
++r;
|
||||
loop = false;
|
||||
|
||||
} else {
|
||||
|
||||
// transform the matrix
|
||||
double sqrt = FastMath.sqrt(c[ir][ir]);
|
||||
b[r][r] = sqrt;
|
||||
double inverse = 1 / sqrt;
|
||||
for (int i = r + 1; i < order; ++i) {
|
||||
int ii = index[i];
|
||||
double e = inverse * c[ii][ir];
|
||||
b[i][r] = e;
|
||||
c[ii][ii] -= e * e;
|
||||
for (int j = r + 1; j < i; ++j) {
|
||||
int ij = index[j];
|
||||
double f = c[ii][ij] - e * b[j][r];
|
||||
c[ii][ij] = f;
|
||||
c[ij][ii] = f;
|
||||
}
|
||||
}
|
||||
|
||||
// prepare next iteration
|
||||
loop = ++r < order;
|
||||
}
|
||||
}
|
||||
|
||||
// build the root matrix
|
||||
rank = r;
|
||||
root = MatrixUtils.createRealMatrix(order, r);
|
||||
for (int i = 0; i < order; ++i) {
|
||||
for (int j = 0; j < r; ++j) {
|
||||
root.setEntry(index[i], j, b[i][j]);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix getRootMatrix() {
|
||||
return root;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public int getRank() {
|
||||
return rank;
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -18,10 +18,9 @@
|
|||
package org.apache.commons.math.random;
|
||||
|
||||
import org.apache.commons.math.exception.DimensionMismatchException;
|
||||
import org.apache.commons.math.linear.NonPositiveDefiniteMatrixException;
|
||||
import org.apache.commons.math.linear.MatrixUtils;
|
||||
import org.apache.commons.math.linear.RealMatrix;
|
||||
import org.apache.commons.math.util.FastMath;
|
||||
import org.apache.commons.math.linear.RectangularCholeskyDecomposition;
|
||||
import org.apache.commons.math.linear.RectangularCholeskyDecompositionImpl;
|
||||
|
||||
/**
|
||||
* A {@link RandomVectorGenerator} that generates vectors with with
|
||||
|
|
@ -68,10 +67,8 @@ public class CorrelatedRandomVectorGenerator
|
|||
private final NormalizedRandomGenerator generator;
|
||||
/** Storage for the normalized vector. */
|
||||
private final double[] normalized;
|
||||
/** Permutated Cholesky root of the covariance matrix. */
|
||||
private RealMatrix root;
|
||||
/** Rank of the covariance matrix. */
|
||||
private int rank;
|
||||
/** Root of the covariance matrix. */
|
||||
private final RealMatrix root;
|
||||
|
||||
/**
|
||||
* Builds a correlated random vector generator from its mean
|
||||
|
|
@ -97,10 +94,13 @@ public class CorrelatedRandomVectorGenerator
|
|||
}
|
||||
this.mean = mean.clone();
|
||||
|
||||
decompose(covariance, small);
|
||||
final RectangularCholeskyDecomposition decomposition =
|
||||
new RectangularCholeskyDecompositionImpl(covariance, small);
|
||||
root = decomposition.getRootMatrix();
|
||||
|
||||
this.generator = generator;
|
||||
normalized = new double[rank];
|
||||
normalized = new double[decomposition.getRank()];
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
|
|
@ -123,10 +123,13 @@ public class CorrelatedRandomVectorGenerator
|
|||
mean[i] = 0;
|
||||
}
|
||||
|
||||
decompose(covariance, small);
|
||||
final RectangularCholeskyDecomposition decomposition =
|
||||
new RectangularCholeskyDecompositionImpl(covariance, small);
|
||||
root = decomposition.getRootMatrix();
|
||||
|
||||
this.generator = generator;
|
||||
normalized = new double[rank];
|
||||
normalized = new double[decomposition.getRank()];
|
||||
|
||||
}
|
||||
|
||||
/** Get the underlying normalized components generator.
|
||||
|
|
@ -136,6 +139,16 @@ public class CorrelatedRandomVectorGenerator
|
|||
return generator;
|
||||
}
|
||||
|
||||
/** Get the rank of the covariance matrix.
|
||||
* The rank is the number of independent rows in the covariance
|
||||
* matrix, it is also the number of columns of the root matrix.
|
||||
* @return rank of the square matrix.
|
||||
* @see #getRootMatrix()
|
||||
*/
|
||||
public int getRank() {
|
||||
return normalized.length;
|
||||
}
|
||||
|
||||
/** Get the root of the covariance matrix.
|
||||
* The root is the rectangular matrix <code>B</code> such that
|
||||
* the covariance matrix is equal to <code>B.B<sup>T</sup></code>
|
||||
|
|
@ -146,122 +159,6 @@ public class CorrelatedRandomVectorGenerator
|
|||
return root;
|
||||
}
|
||||
|
||||
/** Get the rank of the covariance matrix.
|
||||
* The rank is the number of independent rows in the covariance
|
||||
* matrix, it is also the number of columns of the rectangular
|
||||
* matrix of the decomposition.
|
||||
* @return rank of the square matrix.
|
||||
* @see #getRootMatrix()
|
||||
*/
|
||||
public int getRank() {
|
||||
return rank;
|
||||
}
|
||||
|
||||
/** Decompose the original square matrix.
|
||||
* <p>The decomposition is based on a Choleski decomposition
|
||||
* where additional transforms are performed:
|
||||
* <ul>
|
||||
* <li>the rows of the decomposed matrix are permuted</li>
|
||||
* <li>columns with the too small diagonal element are discarded</li>
|
||||
* <li>the matrix is permuted</li>
|
||||
* </ul>
|
||||
* This means that rather than computing M = U<sup>T</sup>.U where U
|
||||
* is an upper triangular matrix, this method computed M=B.B<sup>T</sup>
|
||||
* where B is a rectangular matrix.
|
||||
* @param covariance covariance matrix
|
||||
* @param small diagonal elements threshold under which column are
|
||||
* considered to be dependent on previous ones and are discarded
|
||||
* @throws org.apache.commons.math.linear.NonPositiveDefiniteMatrixException
|
||||
* if the covariance matrix is not strictly positive definite.
|
||||
*/
|
||||
private void decompose(RealMatrix covariance, double small) {
|
||||
int order = covariance.getRowDimension();
|
||||
double[][] c = covariance.getData();
|
||||
double[][] b = new double[order][order];
|
||||
|
||||
int[] swap = new int[order];
|
||||
int[] index = new int[order];
|
||||
for (int i = 0; i < order; ++i) {
|
||||
index[i] = i;
|
||||
}
|
||||
|
||||
rank = 0;
|
||||
for (boolean loop = true; loop;) {
|
||||
|
||||
// find maximal diagonal element
|
||||
swap[rank] = rank;
|
||||
for (int i = rank + 1; i < order; ++i) {
|
||||
int ii = index[i];
|
||||
int isi = index[swap[i]];
|
||||
if (c[ii][ii] > c[isi][isi]) {
|
||||
swap[rank] = i;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// swap elements
|
||||
if (swap[rank] != rank) {
|
||||
int tmp = index[rank];
|
||||
index[rank] = index[swap[rank]];
|
||||
index[swap[rank]] = tmp;
|
||||
}
|
||||
|
||||
// check diagonal element
|
||||
int ir = index[rank];
|
||||
if (c[ir][ir] < small) {
|
||||
|
||||
if (rank == 0) {
|
||||
throw new NonPositiveDefiniteMatrixException(ir, small);
|
||||
}
|
||||
|
||||
// check remaining diagonal elements
|
||||
for (int i = rank; i < order; ++i) {
|
||||
if (c[index[i]][index[i]] < -small) {
|
||||
// there is at least one sufficiently negative diagonal element,
|
||||
// the covariance matrix is wrong
|
||||
throw new NonPositiveDefiniteMatrixException(i, small);
|
||||
}
|
||||
}
|
||||
|
||||
// all remaining diagonal elements are close to zero,
|
||||
// we consider we have found the rank of the covariance matrix
|
||||
++rank;
|
||||
loop = false;
|
||||
|
||||
} else {
|
||||
|
||||
// transform the matrix
|
||||
double sqrt = FastMath.sqrt(c[ir][ir]);
|
||||
b[rank][rank] = sqrt;
|
||||
double inverse = 1 / sqrt;
|
||||
for (int i = rank + 1; i < order; ++i) {
|
||||
int ii = index[i];
|
||||
double e = inverse * c[ii][ir];
|
||||
b[i][rank] = e;
|
||||
c[ii][ii] -= e * e;
|
||||
for (int j = rank + 1; j < i; ++j) {
|
||||
int ij = index[j];
|
||||
double f = c[ii][ij] - e * b[j][rank];
|
||||
c[ii][ij] = f;
|
||||
c[ij][ii] = f;
|
||||
}
|
||||
}
|
||||
|
||||
// prepare next iteration
|
||||
loop = ++rank < order;
|
||||
}
|
||||
}
|
||||
|
||||
// build the root matrix
|
||||
root = MatrixUtils.createRealMatrix(order, rank);
|
||||
for (int i = 0; i < order; ++i) {
|
||||
for (int j = 0; j < rank; ++j) {
|
||||
root.setEntry(index[i], j, b[i][j]);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/** Generate a correlated random vector.
|
||||
* @return a random vector as an array of double. The returned array
|
||||
* is created at each call, the caller can do what it wants with it.
|
||||
|
|
@ -269,7 +166,7 @@ public class CorrelatedRandomVectorGenerator
|
|||
public double[] nextVector() {
|
||||
|
||||
// generate uncorrelated vector
|
||||
for (int i = 0; i < rank; ++i) {
|
||||
for (int i = 0; i < normalized.length; ++i) {
|
||||
normalized[i] = generator.nextNormalizedDouble();
|
||||
}
|
||||
|
||||
|
|
@ -277,11 +174,13 @@ public class CorrelatedRandomVectorGenerator
|
|||
double[] correlated = new double[mean.length];
|
||||
for (int i = 0; i < correlated.length; ++i) {
|
||||
correlated[i] = mean[i];
|
||||
for (int j = 0; j < rank; ++j) {
|
||||
for (int j = 0; j < root.getColumnDimension(); ++j) {
|
||||
correlated[i] += root.getEntry(i, j) * normalized[j];
|
||||
}
|
||||
}
|
||||
|
||||
return correlated;
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -52,9 +52,12 @@ The <action> type attribute can be add,update,fix,remove.
|
|||
If the output is not quite correct, check for invisible trailing spaces!
|
||||
-->
|
||||
<release version="3.0" date="TBD" description="TBD">
|
||||
<action dev="luc" type="add" issue="MATH-541" >
|
||||
Added a "rectangular" Cholesky decomposition for positive semidefinite matrices.
|
||||
</action>
|
||||
<action dev="luc" type="add" issue="MATH-563" >
|
||||
Added setters allowing to change the step size control parameters of adaptive
|
||||
step size ODE integrators
|
||||
step size ODE integrators.
|
||||
</action>
|
||||
<action dev="luc" type="add" issue="MATH-557" >
|
||||
Added a compareTo method to MathUtils that uses a number of ulps as a
|
||||
|
|
|
|||
Loading…
Reference in New Issue