

Java example source code file (RectangularCholeskyDecomposition.java)
The RectangularCholeskyDecomposition.java Java example source code/* * 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/LICENSE2.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.math3.linear; import org.apache.commons.math3.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, depending * on a userdefined 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 * org.apache.commons.math3.random.CorrelatedRandomVectorGenerator correlated * random ndimensions vectors} in a pdimension 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 not suited for solving * linear systems, so it does not provide any {@link DecompositionSolver * decomposition solver}.</p> * * @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld * @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia * @since 2.0 (changed to concrete class in 3.0) */ 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. * <p> * <b>Note: this constructor follows the linpack method to detect dependent * columns by proceeding with the Cholesky algorithm until a nonpositive diagonal * element is encountered. * * @see <a href="http://eprints.ma.man.ac.uk/1193/01/covered/MIMS_ep2008_56.pdf"> * Analysis of the Cholesky Decomposition of a Semidefinite Matrix</a> * * @param matrix Symmetric positive semidefinite matrix. * @exception NonPositiveDefiniteMatrixException if the matrix is not * positive semidefinite. * @since 3.1 */ public RectangularCholeskyDecomposition(RealMatrix matrix) throws NonPositiveDefiniteMatrixException { this(matrix, 0); } /** * Decompose a symmetric positive semidefinite matrix. * * @param matrix Symmetric positive semidefinite matrix. * @param small Diagonal elements threshold under which columns 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 { final int order = matrix.getRowDimension(); final double[][] c = matrix.getData(); final double[][] b = new double[order][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 int swapR = r; for (int i = r + 1; i < order; ++i) { int ii = index[i]; int isr = index[swapR]; if (c[ii][ii] > c[isr][isr]) { swapR = i; } } // swap elements if (swapR != r) { final int tmpIndex = index[r]; index[r] = index[swapR]; index[swapR] = tmpIndex; final double[] tmpRow = b[r]; b[r] = b[swapR]; b[swapR] = tmpRow; } // 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 loop = false; } else { // transform the matrix final double sqrt = FastMath.sqrt(c[ir][ir]); b[r][r] = sqrt; final double inverse = 1 / sqrt; final double inverse2 = 1 / c[ir][ir]; for (int i = r + 1; i < order; ++i) { final int ii = index[i]; final double e = inverse * c[ii][ir]; b[i][r] = e; c[ii][ii] = c[ii][ir] * c[ii][ir] * inverse2; for (int j = r + 1; j < i; ++j) { final int ij = index[j]; final 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 such that * the covariance matrix is equal to <code>B.B^{T} * @return root of the square matrix * @see #getRank() */ 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 * matrix, it is also the number of columns of the rectangular * matrix of the decomposition. * @return r of the square matrix. * @see #getRootMatrix() */ public int getRank() { return rank; } } Other Java examples (source code examples)Here is a short list of links related to this Java RectangularCholeskyDecomposition.java source code file: 
Copyright 19982021 Alvin Alexander, alvinalexander.com
All Rights Reserved.
A percentage of advertising revenue from
pages under the /java/jwarehouse
URI on this website is
paid back to open source projects.