From f837c295cdc997c74e553eb873dae026fec76b1b Mon Sep 17 00:00:00 2001 From: Jurgen Van Gael Date: Sun, 16 May 2010 05:26:14 +0800 Subject: [PATCH] Added Matrix.Transpose with unit tests. Renamed all Decompositions to Factorizations. Moved some Cholesky methods from DenseCholesky.cs to Cholesky.cs. Added unit test which Cholesky factorizes a random pd matrix. --- .../LinearAlgebra/Double/DenseMatrix.cs | 18 ++++++ .../Double/Factorization/Cholesky.cs | 46 +++++++++++--- .../Double/Factorization/DenseCholesky.cs | 60 +++---------------- .../Double/Factorization/ExtensionMethods.cs | 6 +- src/Numerics/LinearAlgebra/Double/Matrix.cs | 17 ++++++ .../Double/Factorization/CholeskyTests.cs | 59 ++++++++++++++++-- .../LinearAlgebraTests/Double/MatrixTests.cs | 24 ++++++++ 7 files changed, 162 insertions(+), 68 deletions(-) diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs index c66c6506..a90ead04 100644 --- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs @@ -196,6 +196,24 @@ namespace MathNet.Numerics.LinearAlgebra.Double Array.Clear(Data, 0, Data.Length); } + /// + /// Returns the transpose of this matrix. + /// + /// The transpose of this matrix. + public override Matrix Transpose() + { + DenseMatrix ret = new DenseMatrix(ColumnCount, RowCount); + for (int j = 0; j < ColumnCount; j++) + { + int index = j * RowCount; + for (int i = 0; i < RowCount; i++) + { + ret.Data[i * ColumnCount + j] = Data[index + i]; + } + } + return ret; + } + #region Elementary operations /// diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs b/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs index 3246596c..47e0c652 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs @@ -28,22 +28,27 @@ // OTHER DEALINGS IN THE SOFTWARE. // -namespace MathNet.Numerics.LinearAlgebra.Double.Decomposition +namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { using System; using Properties; /// - /// A class which encapsulates the functionality of a Cholesky decomposition. - /// For a symmetric, positive definite matrix A, the Cholesky decomposition + /// A class which encapsulates the functionality of a Cholesky factorization. + /// For a symmetric, positive definite matrix A, the Cholesky factorization /// is an lower triangular matrix L so that A = L*L'. /// /// - /// The computation of the Cholesky decomposition is done at construction time. If the matrix is not symmetric + /// The computation of the Cholesky factorization is done at construction time. If the matrix is not symmetric /// or positive definite, the constructor will throw an exception. /// public abstract class Cholesky { + /// + /// Stores the Cholesky factor. + /// + protected Matrix mFactor; + /// /// Internal method which routes the call to perform the Cholesky factorization to the appropriate class. /// @@ -61,18 +66,43 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Decomposition } /// - /// A reference to the lower triangular form of the Cholesky matrix. + /// Returns the lower triangular form of the Cholesky matrix. /// - public abstract Matrix Factor { get; } + public virtual Matrix Factor + { + get { return mFactor; } + } /// /// The determinant of the matrix for which the Cholesky matrix was computed. /// - public abstract double Determinant { get; } + public virtual double Determinant + { + get + { + double det = 1.0; + for (int j = 0; j < mFactor.RowCount; j++) + { + det *= (mFactor[j, j] * mFactor[j, j]); + } + return det; + } + } /// /// The log determinant of the matrix for which the Cholesky matrix was computed. /// - public abstract double DeterminantLn { get; } + public virtual double DeterminantLn + { + get + { + double det = 0.0; + for (int j = 0; j < mFactor.RowCount; j++) + { + det += (2.0 * Math.Log(mFactor[j, j])); + } + return det; + } + } } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs index 1f03a90c..ad1e1771 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs @@ -28,30 +28,25 @@ // OTHER DEALINGS IN THE SOFTWARE. // -namespace MathNet.Numerics.LinearAlgebra.Double.Decomposition +namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { using System; using Properties; /// - /// A class which encapsulates the functionality of a Cholesky decomposition for dense matrices. - /// For a symmetric, positive definite matrix A, the Cholesky decomposition + /// A class which encapsulates the functionality of a Cholesky factorization for dense matrices. + /// For a symmetric, positive definite matrix A, the Cholesky factorization /// is an lower triangular matrix L so that A = L*L'. /// /// - /// The computation of the Cholesky decomposition is done at construction time. If the matrix is not symmetric + /// The computation of the Cholesky factorization is done at construction time. If the matrix is not symmetric /// or positive definite, the constructor will throw an exception. /// public class DenseCholesky : Cholesky { - /// - /// Stores the Cholesky factor for the decomposition. - /// - private DenseMatrix _factor; - /// /// Initializes a new instance of the class. This object will compute the - /// Cholesky decomposition when the constructor is called and cache it's factorization. + /// Cholesky factorization when the constructor is called and cache it's factorization. /// /// The matrix to factor. /// If is null. @@ -70,48 +65,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Decomposition } // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting). - _factor = (DenseMatrix)matrix.Clone(); - Control.LinearAlgebraProvider.CholeskyFactor(_factor.Data, _factor.RowCount); - } - - /// - /// Returns the lower triangular form of the Cholesky matrix. - /// - public override Matrix Factor - { - get { return _factor; } - } - - /// - /// The determinant of the matrix for which the Cholesky matrix was computed. - /// - public override double Determinant - { - get - { - double det = 1.0; - for (int j = 0; j < _factor.RowCount; j++) - { - det *= (_factor[j, j] * _factor[j, j]); - } - return det; - } - } - - /// - /// The log determinant of the matrix for which the Cholesky matrix was computed. - /// - public override double DeterminantLn - { - get - { - double det = 0.0; - for (int j = 0; j < _factor.RowCount; j++) - { - det += (2.0 * Math.Log(_factor[j, j])); - } - return det; - } + var factor = (DenseMatrix)matrix.Clone(); + Control.LinearAlgebraProvider.CholeskyFactor(factor.Data, factor.RowCount); + mFactor = factor; } } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/ExtensionMethods.cs b/src/Numerics/LinearAlgebra/Double/Factorization/ExtensionMethods.cs index d2653edb..fe56050e 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/ExtensionMethods.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/ExtensionMethods.cs @@ -28,13 +28,13 @@ // OTHER DEALINGS IN THE SOFTWARE. // -namespace MathNet.Numerics.LinearAlgebra.Double.Decomposition +namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { using System; using Properties; /// - /// Extension methods which return decompositions for the various matrix classes. + /// Extension methods which return factorizations for the various matrix classes. /// public static class ExtensionMethods { @@ -45,7 +45,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Decomposition /// The Cholesky decomposition object. public static Cholesky Cholesky(this Matrix matrix) { - return Decomposition.Cholesky.Create(matrix); + return Factorization.Cholesky.Create(matrix); } } } diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.cs b/src/Numerics/LinearAlgebra/Double/Matrix.cs index 9d271e9c..ca2f0a57 100644 --- a/src/Numerics/LinearAlgebra/Double/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Double/Matrix.cs @@ -639,5 +639,22 @@ namespace MathNet.Numerics.LinearAlgebra.Double } } } + + /// + /// Returns the transpose of this matrix. + /// + /// The transpose of this matrix. + public virtual Matrix Transpose() + { + Matrix ret = CreateMatrix(ColumnCount, RowCount); + for (int j = 0; j < ColumnCount; j++) + { + for (int i = 0; i < RowCount; i++) + { + ret.At(j, i, At(i, j)); + } + } + return ret; + } } } \ No newline at end of file diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/CholeskyTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/CholeskyTests.cs index e1a43973..bfeb81ca 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Factorization/CholeskyTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/CholeskyTests.cs @@ -28,12 +28,12 @@ // OTHER DEALINGS IN THE SOFTWARE. // -namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Decomposition +namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization { using System.Collections.Generic; using MbUnit.Framework; using LinearAlgebra.Double; - using LinearAlgebra.Double.Decomposition; + using LinearAlgebra.Double.Factorization; public class CholeskyTests { @@ -41,8 +41,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Decomposition [Row(1)] [Row(10)] [Row(100)] - [Row(1000)] - public void CanDecomposeIdentity(int order) + public void CanFactorizeIdentity(int order) { var I = DenseMatrix.Identity(order); var C = I.Cholesky(); @@ -86,7 +85,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Decomposition [Row(1)] [Row(10)] [Row(100)] - [Row(1000)] public void IdentityDeterminantIsOne(int order) { var I = DenseMatrix.Identity(order); @@ -94,5 +92,56 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Decomposition Assert.AreEqual(1.0, C.Determinant); Assert.AreEqual(0.0, C.DeterminantLn); } + + [Test] + [Row(1)] + [Row(2)] + [Row(5)] + [Row(10)] + [Row(50)] + [Row(100)] + [MultipleAsserts] + public void CanFactorizeRandomMatrix(int order) + { + // Fill a matrix with standard random numbers. + var normal = new Distributions.Normal(); + normal.RandomSource = new Random.MersenneTwister(1); + var A = new DenseMatrix(order); + for (int i = 0; i < order; i++) + { + for (int j = 0; j < order; j++) + { + A[i, j] = normal.Sample(); + } + } + + // Generate a matrix which is positive definite. + var X = A.Transpose() * A; + var chol = X.Cholesky(); + var C = chol.Factor; + + // Make sure the Cholesky factor has the right dimensions. + Assert.AreEqual(order, C.RowCount); + Assert.AreEqual(order, C.ColumnCount); + + // Make sure the Cholesky factor is lower triangular. + for (int i = 0; i < C.RowCount; i++) + { + for (int j = i+1; j < C.ColumnCount; j++) + { + Assert.AreEqual(0.0, C[i, j]); + } + } + + // Make sure the cholesky factor times it's transpose is the original matrix. + var XfromC = C * C.Transpose(); + for (int i = 0; i < XfromC.RowCount; i++) + { + for (int j = 0; j < XfromC.ColumnCount; j++) + { + Assert.AreApproximatelyEqual(X[i,j], XfromC[i, j], 1.0e-13); + } + } + } } } diff --git a/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.cs index d6b53203..14d683d1 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.cs @@ -457,5 +457,29 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double Vector result = CreateVector(matrix.RowCount - 1); matrix.GetColumn(0, 0, matrix.RowCount, result); } + + [Test] + [Row("Singular3x3")] + [Row("Square3x3")] + [Row("Square4x4")] + [Row("Tall3x2")] + [Row("Wide2x3")] + [MultipleAsserts] + public void CanTransposeMatrix(string name) + { + var matrix = CreateMatrix(testData2D[name]); + var transpose = matrix.Transpose(); + + Assert.AreNotSame(matrix, transpose); + Assert.AreEqual(matrix.RowCount, transpose.ColumnCount); + Assert.AreEqual(matrix.ColumnCount, transpose.RowCount); + for (var i = 0; i < matrix.RowCount; i++) + { + for (var j = 0; j < matrix.ColumnCount; j++) + { + Assert.AreEqual(matrix[i, j], transpose[j, i]); + } + } + } } } \ No newline at end of file