Permalink
Fetching contributors…
Cannot retrieve contributors at this time
753 lines (638 sloc) 26.8 KB
// Accord Math Library
// The Accord.NET Framework
// http://accord-framework.net
//
// Copyright © César Souza, 2009-2017
// cesarsouza at gmail.com
//
// Original work copyright © Lutz Roeder, 2000
// Adapted from Mapack for .NET, September 2000
// Adapted from Mapack for COM and Jama routines
// http://www.aisto.com/roeder/dotnet
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
//
namespace Accord.Math.Decompositions
{
using System;
using Accord.Compat;
using System.Threading.Tasks;
/// <summary>
/// Cholesky Decomposition of a symmetric, positive definite matrix.
/// </summary>
///
/// <remarks>
/// <para>
/// For a symmetric, positive definite matrix <c>A</c>, the Cholesky decomposition is a
/// lower triangular matrix <c>L</c> so that <c>A = L * L'</c>.
/// If the matrix is not positive definite, the constructor returns a partial
/// decomposition and sets two internal variables that can be queried using the
/// <see cref="IsUndefined"/> and <see cref="IsPositiveDefinite"/> properties.</para>
/// <para>
/// Any square matrix A with non-zero pivots can be written as the product of a
/// lower triangular matrix L and an upper triangular matrix U; this is called
/// the LU decomposition. However, if A is symmetric and positive definite, we
/// can choose the factors such that U is the transpose of L, and this is called
/// the Cholesky decomposition. Both the LU and the Cholesky decomposition are
/// used to solve systems of linear equations.</para>
/// <para>
/// When it is applicable, the Cholesky decomposition is twice as efficient
/// as the LU decomposition.</para>
/// </remarks>
///
[Serializable]
public sealed class CholeskyDecompositionF : ICloneable, ISolverMatrixDecomposition<Single>
{
private Single[,] L;
private Single[] D;
private int n;
private bool positiveDefinite;
private bool destroyed;
private bool robust;
private bool undefined;
// cache for lazy evaluation
private Single[,] leftTriangularFactor;
private Single[,] diagonalMatrix;
private Single? determinant;
private double? lndeterminant;
private bool? nonsingular;
/// <summary>
/// Constructs a new Cholesky Decomposition.
/// </summary>
///
/// <param name="value">
/// The symmetric matrix, given in upper triangular form, to be decomposed.</param>
/// <param name="robust">
/// True to perform a square-root free LDLt decomposition, false otherwise.</param>
/// <param name="inPlace">
/// True to perform the decomposition in place, storing the factorization in the
/// lower triangular part of the given matrix.</param>
/// <param name="valueType">
/// How to interpret the matrix given to be decomposed. Using this parameter, a lower or
/// upper-triangular matrix can be interpreted as a symmetric matrix by assuming both lower
/// and upper parts contain the same elements. Use this parameter in conjunction with inPlace
/// to save memory by storing the original matrix and its decomposition at the same memory
/// location (lower part will contain the decomposition's L matrix, upper part will contains
/// the original matrix).</param>
///
public CholeskyDecompositionF(Single[,] value, bool robust = false,
bool inPlace = false, MatrixType valueType = MatrixType.UpperTriangular)
{
if (value.Rows() != value.Columns())
throw new DimensionMismatchException("value", "Matrix is not square.");
if (!inPlace)
value = value.Copy();
this.n = value.Rows();
this.L = value.ToUpperTriangular(valueType, result: value);
this.robust = robust;
if (robust)
{
LDLt(); // Compute square-root free decomposition
}
else
{
LLt(); // Compute standard Cholesky decomposition
}
}
/// <summary>
/// Gets whether the decomposed matrix was positive definite.
/// </summary>
///
public bool IsPositiveDefinite
{
get { return this.positiveDefinite; }
}
/// <summary>
/// Gets a value indicating whether the LDLt factorization
/// has been computed successfully or if it is undefined.
/// </summary>
///
/// <value>
/// <c>true</c> if the factorization is not defined; otherwise, <c>false</c>.
/// </value>
///
public bool IsUndefined
{
get { return this.undefined; }
}
/// <summary>
/// Gets the left (lower) triangular factor
/// <c>L</c> so that <c>A = L * D * L'</c>.
/// </summary>
///
public Single[,] LeftTriangularFactor
{
get
{
if (leftTriangularFactor == null)
{
if (destroyed)
throw new InvalidOperationException("The decomposition has been destroyed.");
if (undefined)
throw new InvalidOperationException("The decomposition is undefined (zero in diagonal).");
leftTriangularFactor = L.GetLowerTriangle();
}
return leftTriangularFactor;
}
}
/// <summary>
/// Gets the block diagonal matrix of diagonal elements in a LDLt decomposition.
/// </summary>
///
public Single[,] DiagonalMatrix
{
get
{
if (diagonalMatrix == null)
{
if (destroyed)
throw new InvalidOperationException("The decomposition has been destroyed.");
diagonalMatrix = Matrix.Diagonal(D);
}
return diagonalMatrix;
}
}
/// <summary>
/// Gets the one-dimensional array of diagonal elements in a LDLt decomposition.
/// </summary>
///
public Single[] Diagonal
{
get { return D; }
}
/// <summary>
/// Gets the determinant of the decomposed matrix.
/// </summary>
///
public Single Determinant
{
get
{
if (!determinant.HasValue)
{
if (destroyed)
throw new InvalidOperationException("The decomposition has been destroyed.");
if (undefined)
throw new InvalidOperationException("The decomposition is undefined (zero in diagonal).");
Single detL = 1, detD = 1;
for (int i = 0; i < n; i++)
detL *= L[i, i];
if (D != null)
{
for (int i = 0; i < n; i++)
detD *= D[i];
}
determinant = detL * detL * detD;
}
return determinant.Value;
}
}
/// <summary>
/// If the matrix is positive-definite, gets the
/// log-determinant of the decomposed matrix.
/// </summary>
///
public double LogDeterminant
{
get
{
if (!lndeterminant.HasValue)
{
if (destroyed)
throw new InvalidOperationException("The decomposition has been destroyed.");
if (undefined)
throw new InvalidOperationException("The decomposition is undefined (zero in diagonal).");
double detL = 0, detD = 0;
for (int i = 0; i < n; i++)
detL += Math.Log((double)L[i, i]);
if (D != null)
{
for (int i = 0; i < D.Length; i++)
detD += Math.Log((double)D[i]);
}
lndeterminant = detL + detL + detD;
}
return lndeterminant.Value;
}
}
/// <summary>
/// Gets a value indicating whether the decomposed
/// matrix is non-singular (i.e. invertible).
/// </summary>
///
public bool Nonsingular
{
get
{
if (!nonsingular.HasValue)
{
if (destroyed)
throw new InvalidOperationException("The decomposition has been destroyed.");
if (undefined)
throw new InvalidOperationException("The decomposition is undefined (zero in diagonal).");
bool nonSingular = true;
for (int i = 0; i < n && nonSingular; i++)
if (L[i, i] == 0 || D[i] == 0) nonSingular = false;
nonsingular = nonSingular;
}
return nonsingular.Value;
}
}
private unsafe void LLt()
{
D = Vector.Ones<Single>(n);
this.positiveDefinite = true;
for (int j = 0; j < n; j++)
{
Single s = 0;
for (int k = 0; k < j; k++)
{
Single t = L[k, j];
for (int i = 0; i < k; i++)
t -= L[j, i] * L[k, i];
t = t / L[k, k];
L[j, k] = t;
s += t * t;
}
s = L[j, j] - s;
// Use a tolerance for positive-definiteness
this.positiveDefinite &= (s > (Single)1e-14 * Math.Abs(L[j, j]));
L[j, j] = (Single)Math.Sqrt((double)s);
}
}
private unsafe void LDLt()
{
D = new Single[n];
this.positiveDefinite = true;
Single[] v = new Single[n];
for (int i = 0; i < n; i++)
{
for (int j = 0; j < i; j++)
v[j] = L[i, j] * D[j];
Single d = 0;
for (int k = 0; k < i; k++)
d += L[i, k] * v[k];
d = D[i] = v[i] = L[i, i] - d;
// Use a tolerance for positive-definiteness
this.positiveDefinite &= (v[i] > (Single)1e-14 * Math.Abs(L[i, i]));
// If one of the diagonal elements is zero, the
// decomposition (without pivoting) is undefined.
if (v[i] == 0) { undefined = true; return; }
Parallel.For(i + 1, n, k =>
{
Single s = 0;
for (int j = 0; j < i; j++)
s += L[k, j] * v[j];
L[k, i] = (L[i, k] - s) / d;
});
}
for (int i = 0; i < n; i++)
L[i, i] = 1;
}
/// <summary>
/// Solves a set of equation systems of type <c>A * X = B</c>.
/// </summary>
///
/// <param name="value">Right hand side matrix with as many rows as <c>A</c> and any number of columns.</param>
/// <returns>Matrix <c>X</c> so that <c>L * L' * X = B</c>.</returns>
/// <exception cref="T:System.ArgumentException">Matrix dimensions do not match.</exception>
/// <exception cref="T:System.NonSymmetricMatrixException">Matrix is not symmetric.</exception>
/// <exception cref="T:System.NonPositiveDefiniteMatrixException">Matrix is not positive-definite.</exception>
///
public Single[,] Solve(Single[,] value)
{
return Solve(value, false);
}
/// <summary>
/// Solves a set of equation systems of type <c>A * X = B</c>.
/// </summary>
///
/// <param name="value">Right hand side matrix with as many rows as <c>A</c> and any number of columns.</param>
/// <returns>Matrix <c>X</c> so that <c>L * L' * X = B</c>.</returns>
/// <exception cref="T:System.ArgumentException">Matrix dimensions do not match.</exception>
/// <exception cref="T:System.NonSymmetricMatrixException">Matrix is not symmetric.</exception>
/// <exception cref="T:System.NonPositiveDefiniteMatrixException">Matrix is not positive-definite.</exception>
/// <param name="inPlace">True to compute the solving in place, false otherwise.</param>
///
public Single[,] Solve(Single[,] value, bool inPlace)
{
if (value == null)
throw new ArgumentNullException("value");
if (value.Rows() != n)
throw new ArgumentException("Argument matrix should have the same number of rows as the decomposed matrix.", "value");
if (!robust && !positiveDefinite)
throw new NonPositiveDefiniteMatrixException("Decomposed matrix is not positive definite.");
if (undefined)
throw new InvalidOperationException("The decomposition is undefined (zero in diagonal).");
if (destroyed)
throw new InvalidOperationException("The decomposition has been destroyed.");
// int count = value.Columns();
var B = inPlace ? value : value.MemberwiseClone();
int m = B.Columns();
// Solve L*Y = B;
for (int k = 0; k < n; k++)
{
for (int j = 0; j < m; j++)
{
for (int i = 0; i < k; i++)
B[k, j] -= B[i, j] * L[k, i];
B[k, j] /= L[k, k];
}
}
if (robust)
{
for (int k = 0; k < D.Length; k++)
for (int j = 0; j < m; j++)
B[k, j] /= D[k];
}
// Solve L'*X = Y;
for (int k = n - 1; k >= 0; k--)
{
for (int j = 0; j < m; j++)
{
for (int i = k + 1; i < n; i++)
B[k, j] -= B[i, j] * L[i, k];
B[k, j] /= L[k, k];
}
}
return B;
}
/// <summary>
/// Solves a set of equation systems of type <c>A * x = b</c>.
/// </summary>
///
/// <param name="value">Right hand side column vector with as many rows as <c>A</c>.</param>
/// <returns>Vector <c>x</c> so that <c>L * L' * x = b</c>.</returns>
/// <exception cref="T:System.ArgumentException">Matrix dimensions do not match.</exception>
/// <exception cref="T:System.NonSymmetricMatrixException">Matrix is not symmetric.</exception>
/// <exception cref="T:System.NonPositiveDefiniteMatrixException">Matrix is not positive-definite.</exception>
///
public Single[] Solve(Single[] value)
{
return Solve(value, false);
}
/// <summary>
/// Solves a set of equation systems of type <c>A * x = b</c>.
/// </summary>
///
/// <param name="value">Right hand side column vector with as many rows as <c>A</c>.</param>
/// <returns>Vector <c>x</c> so that <c>L * L' * x = b</c>.</returns>
/// <exception cref="T:System.ArgumentException">Matrix dimensions do not match.</exception>
/// <exception cref="T:System.NonSymmetricMatrixException">Matrix is not symmetric.</exception>
/// <exception cref="T:System.NonPositiveDefiniteMatrixException">Matrix is not positive-definite.</exception>
/// <param name="inPlace">True to compute the solving in place, false otherwise.</param>
///
public Single[] Solve(Single[] value, bool inPlace)
{
if (value == null)
throw new ArgumentNullException("value");
if (value.Length != n)
throw new ArgumentException("Argument vector should have the same length as rows in the decomposed matrix.", "value");
if (!robust && !positiveDefinite)
throw new NonPositiveDefiniteMatrixException("Decomposed matrix is not positive definite.");
if (undefined)
throw new InvalidOperationException("The decomposition is undefined (zero in diagonal).");
if (destroyed)
throw new InvalidOperationException("The decomposition has been destroyed.");
var B = inPlace ? value : value.Copy();
// Solve L*Y = B;
for (int k = 0; k < n; k++)
{
for (int i = 0; i < k; i++)
B[k] -= B[i] * L[k, i];
B[k] /= L[k, k];
}
if (robust)
{
for (int k = 0; k < D.Length; k++)
B[k] /= D[k];
}
// Solve L'*X = Y;
for (int k = n - 1; k >= 0; k--)
{
for (int i = k + 1; i < n; i++)
B[k] -= B[i] * L[i, k];
B[k] /= L[k, k];
}
return B;
}
/// <summary>
/// Solves a set of equation systems of type <c>A * X = I</c>.
/// </summary>
///
public Single[,] Inverse()
{
return Solve(Matrix.Identity<Single>(n));
}
/// <summary>
/// Computes the diagonal of the inverse of the decomposed matrix.
/// </summary>
///
public Single[] InverseDiagonal(bool destroy = false)
{
return InverseDiagonal(new Single[n], destroy);
}
/// <summary>
/// Computes the diagonal of the inverse of the decomposed matrix.
/// </summary>
///
/// <param name="destroy">True to conserve memory by reusing the
/// same space used to hold the decomposition, thus destroying
/// it in the process. Pass false otherwise.</param>
/// <param name="result">The array to hold the result of the
/// computation. Should be of same length as the the diagonal
/// of the original matrix.</param>
///
public Single[] InverseDiagonal(Single[] result, bool destroy = false)
{
if (!robust && !positiveDefinite)
throw new NonPositiveDefiniteMatrixException("Matrix is not positive definite.");
if (undefined)
throw new InvalidOperationException("The decomposition is undefined (zero in diagonal).");
if (destroyed)
throw new InvalidOperationException("The decomposition has been destroyed.");
Single[,] S;
if (destroy)
{
S = L; destroyed = true;
}
else
{
S = Matrix.Zeros<Single>(n, n);
}
// References:
// http://books.google.com/books?id=myzIPBwyBbcC&pg=PA119
// Compute the inverse S of the lower triangular matrix L
// and store in place of the upper triangular part of S.
for (int j = n - 1; j >= 0; j--)
{
S[j, j] = 1 / L[j, j];
for (int i = j - 1; i >= 0; i--)
{
Single sum = 0;
for (int k = i + 1; k <= j; k++)
sum += L[k, i] * S[k, j];
S[i, j] = -sum / L[i, i];
}
}
// Compute the 2-norm squared of the rows
// of the upper (right) triangular matrix S.
if (robust)
{
for (int i = 0; i < n; i++)
{
Single sum = 0;
for (int j = i; j < n; j++)
sum += S[i, j] * S[i, j] / D[j];
result[i] = sum;
}
}
else
{
for (int i = 0; i < n; i++)
{
Single sum = 0;
for (int j = i; j < n; j++)
sum += S[i, j] * S[i, j];
result[i] = sum;
}
}
return result;
}
/// <summary>
/// Computes the trace of the inverse of the decomposed matrix.
/// </summary>
///
/// <param name="destroy">True to conserve memory by reusing the
/// same space used to hold the decomposition, thus destroying
/// it in the process. Pass false otherwise.</param>
///
public Single InverseTrace(bool destroy = false)
{
if (!robust && !positiveDefinite)
throw new NonPositiveDefiniteMatrixException("Matrix is not positive definite.");
if (undefined)
throw new InvalidOperationException("The decomposition is undefined (zero in diagonal).");
if (destroyed)
throw new InvalidOperationException("The decomposition has been destroyed.");
Single[,] S;
if (destroy)
{
S = L; destroyed = true;
}
else
{
S = Matrix.Zeros<Single>(n, n);
}
// References:
// http://books.google.com/books?id=myzIPBwyBbcC&pg=PA119
// Compute the inverse S of the lower triangular matrix L
// and store in place of the upper triangular part of S.
for (int j = n - 1; j >= 0; j--)
{
S[j, j] = 1 / L[j, j];
for (int i = j - 1; i >= 0; i--)
{
Single sum = 0;
for (int k = i + 1; k <= j; k++)
sum += L[k, i] * S[k, j];
S[i, j] = -sum / L[i, i];
}
}
// Compute the 2-norm squared of the rows
// of the upper (right) triangular matrix S.
Single trace = 0;
if (robust)
{
for (int i = 0; i < n; i++)
for (int j = i; j < n; j++)
trace += S[i, j] * S[i, j] / D[j];
}
else
{
for (int i = 0; i < n; i++)
for (int j = i; j < n; j++)
trace += S[i, j] * S[i, j];
}
return trace;
}
/// <summary>
/// Reverses the decomposition, reconstructing the original matrix <c>X</c>.
/// </summary>
///
public Single[,] Reverse()
{
if (destroyed)
throw new InvalidOperationException("The decomposition has been destroyed.");
if (undefined)
throw new InvalidOperationException("The decomposition is undefined (zero in diagonal).");
if (robust)
return LeftTriangularFactor.Dot(DiagonalMatrix).DotWithTransposed(LeftTriangularFactor);
return LeftTriangularFactor.DotWithTransposed(LeftTriangularFactor);
}
/// <summary>
/// Computes <c>(Xt * X)^1</c> (the inverse of the covariance matrix). This
/// matrix can be used to determine standard errors for the coefficients when
/// solving a linear set of equations through any of the <see cref="Solve(Single[,])"/>
/// methods.
/// </summary>
///
public Single[,] GetInformationMatrix()
{
var X = Reverse();
return X.TransposeAndDot(X).Inverse();
}
/// <summary>
/// Creates a new Cholesky decomposition directly from
/// an already computed left triangular matrix <c>L</c>.
/// </summary>
/// <param name="leftTriangular">The left triangular matrix from a Cholesky decomposition.</param>
///
public static CholeskyDecompositionF FromLeftTriangularMatrix(Single[,] leftTriangular)
{
var chol = new CholeskyDecompositionF();
chol.n = leftTriangular.Rows();
chol.L = leftTriangular;
chol.positiveDefinite = true;
chol.robust = false;
chol.D = Vector.Ones<Single>(chol.n);
return chol;
}
#region ICloneable Members
private CholeskyDecompositionF()
{
}
/// <summary>
/// Creates a new object that is a copy of the current instance.
/// </summary>
/// <returns>
/// A new object that is a copy of this instance.
/// </returns>
///
public object Clone()
{
var clone = new CholeskyDecompositionF();
clone.L = L.MemberwiseClone();
clone.D = (Single[])D.Clone();
clone.destroyed = destroyed;
clone.n = n;
clone.undefined = undefined;
clone.robust = robust;
clone.positiveDefinite = positiveDefinite;
return clone;
}
#endregion
}
}