Merged RectangularCholeskyDecomposition and RectangularCholeskyDecompositionImpl (see MATH-662).
git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1175105 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Sebastien Brisard
parent
18323639c9
commit
b4764661a3
|
|
@ -17,11 +17,10 @@
|
|||
|
||||
package org.apache.commons.math.linear;
|
||||
|
||||
import org.apache.commons.math.util.FastMath;
|
||||
|
||||
/**
|
||||
* An interface to classes that implement an algorithm to calculate a
|
||||
* rectangular variation of Cholesky decomposition of a real symmetric
|
||||
* positive semidefinite matrix.
|
||||
* 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
|
||||
|
|
@ -38,12 +37,118 @@ package org.apache.commons.math.linear;
|
|||
* linear systems, so it does not provide any {@link DecompositionSolver
|
||||
* decomposition solver}.</p>
|
||||
*
|
||||
* @see CholeskyDecomposition
|
||||
* @see org.apache.commons.math.random.CorrelatedRandomVectorGenerator
|
||||
* @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a>
|
||||
* @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a>
|
||||
* @version $Id$
|
||||
* @since 3.0
|
||||
* @since 2.0 (changed to concrete class in 3.0)
|
||||
*/
|
||||
public interface RectangularCholeskyDecomposition {
|
||||
public class 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 RectangularCholeskyDecomposition(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(c[ir][ir], 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(c[index[i]][index[i]], 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]);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/** Get the root of the covariance matrix.
|
||||
* The root is the rectangular matrix <code>B</code> such that
|
||||
|
|
@ -51,7 +156,9 @@ public interface RectangularCholeskyDecomposition {
|
|||
* @return root of the square matrix
|
||||
* @see #getRank()
|
||||
*/
|
||||
RealMatrix getRootMatrix();
|
||||
public RealMatrix getRootMatrix() {
|
||||
return root;
|
||||
}
|
||||
|
||||
/** Get the rank of the symmetric positive semidefinite matrix.
|
||||
* The r is the number of independent rows in the symmetric positive semidefinite
|
||||
|
|
@ -60,6 +167,8 @@ public interface RectangularCholeskyDecomposition {
|
|||
* @return r of the square matrix.
|
||||
* @see #getRootMatrix()
|
||||
*/
|
||||
int getRank();
|
||||
public int getRank() {
|
||||
return rank;
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,152 +0,0 @@
|
|||
/*
|
||||
* 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 $Id$
|
||||
* @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(c[ir][ir], 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(c[index[i]][index[i]], 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;
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -20,7 +20,6 @@ package org.apache.commons.math.random;
|
|||
import org.apache.commons.math.exception.DimensionMismatchException;
|
||||
import org.apache.commons.math.linear.RealMatrix;
|
||||
import org.apache.commons.math.linear.RectangularCholeskyDecomposition;
|
||||
import org.apache.commons.math.linear.RectangularCholeskyDecompositionImpl;
|
||||
|
||||
/**
|
||||
* A {@link RandomVectorGenerator} that generates vectors with with
|
||||
|
|
@ -95,7 +94,7 @@ public class CorrelatedRandomVectorGenerator
|
|||
this.mean = mean.clone();
|
||||
|
||||
final RectangularCholeskyDecomposition decomposition =
|
||||
new RectangularCholeskyDecompositionImpl(covariance, small);
|
||||
new RectangularCholeskyDecomposition(covariance, small);
|
||||
root = decomposition.getRootMatrix();
|
||||
|
||||
this.generator = generator;
|
||||
|
|
@ -124,7 +123,7 @@ public class CorrelatedRandomVectorGenerator
|
|||
}
|
||||
|
||||
final RectangularCholeskyDecomposition decomposition =
|
||||
new RectangularCholeskyDecompositionImpl(covariance, small);
|
||||
new RectangularCholeskyDecomposition(covariance, small);
|
||||
root = decomposition.getRootMatrix();
|
||||
|
||||
this.generator = generator;
|
||||
|
|
|
|||
Loading…
Reference in New Issue