Browse Source

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.
la-knuth
Jurgen Van Gael 17 years ago
parent
commit
f837c295cd
  1. 18
      src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
  2. 46
      src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs
  3. 60
      src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs
  4. 6
      src/Numerics/LinearAlgebra/Double/Factorization/ExtensionMethods.cs
  5. 17
      src/Numerics/LinearAlgebra/Double/Matrix.cs
  6. 59
      src/UnitTests/LinearAlgebraTests/Double/Factorization/CholeskyTests.cs
  7. 24
      src/UnitTests/LinearAlgebraTests/Double/MatrixTests.cs

18
src/Numerics/LinearAlgebra/Double/DenseMatrix.cs

@ -196,6 +196,24 @@ namespace MathNet.Numerics.LinearAlgebra.Double
Array.Clear(Data, 0, Data.Length); Array.Clear(Data, 0, Data.Length);
} }
/// <summary>
/// Returns the transpose of this matrix.
/// </summary>
/// <returns>The transpose of this matrix.</returns>
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 #region Elementary operations
/// <summary> /// <summary>

46
src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs

@ -28,22 +28,27 @@
// OTHER DEALINGS IN THE SOFTWARE. // OTHER DEALINGS IN THE SOFTWARE.
// </copyright> // </copyright>
namespace MathNet.Numerics.LinearAlgebra.Double.Decomposition namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
{ {
using System; using System;
using Properties; using Properties;
/// <summary> /// <summary>
/// <para>A class which encapsulates the functionality of a Cholesky decomposition.</para> /// <para>A class which encapsulates the functionality of a Cholesky factorization.</para>
/// <para>For a symmetric, positive definite matrix A, the Cholesky decomposition /// <para>For a symmetric, positive definite matrix A, the Cholesky factorization
/// is an lower triangular matrix L so that A = L*L'.</para> /// is an lower triangular matrix L so that A = L*L'.</para>
/// </summary> /// </summary>
/// <remarks> /// <remarks>
/// 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. /// or positive definite, the constructor will throw an exception.
/// </remarks> /// </remarks>
public abstract class Cholesky public abstract class Cholesky
{ {
/// <summary>
/// Stores the Cholesky factor.
/// </summary>
protected Matrix mFactor;
/// <summary> /// <summary>
/// Internal method which routes the call to perform the Cholesky factorization to the appropriate class. /// Internal method which routes the call to perform the Cholesky factorization to the appropriate class.
/// </summary> /// </summary>
@ -61,18 +66,43 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Decomposition
} }
/// <summary> /// <summary>
/// A reference to the lower triangular form of the Cholesky matrix. /// Returns the lower triangular form of the Cholesky matrix.
/// </summary> /// </summary>
public abstract Matrix Factor { get; } public virtual Matrix Factor
{
get { return mFactor; }
}
/// <summary> /// <summary>
/// The determinant of the matrix for which the Cholesky matrix was computed. /// The determinant of the matrix for which the Cholesky matrix was computed.
/// </summary> /// </summary>
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;
}
}
/// <summary> /// <summary>
/// The log determinant of the matrix for which the Cholesky matrix was computed. /// The log determinant of the matrix for which the Cholesky matrix was computed.
/// </summary> /// </summary>
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;
}
}
} }
} }

60
src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs

@ -28,30 +28,25 @@
// OTHER DEALINGS IN THE SOFTWARE. // OTHER DEALINGS IN THE SOFTWARE.
// </copyright> // </copyright>
namespace MathNet.Numerics.LinearAlgebra.Double.Decomposition namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
{ {
using System; using System;
using Properties; using Properties;
/// <summary> /// <summary>
/// <para>A class which encapsulates the functionality of a Cholesky decomposition for dense matrices.</para> /// <para>A class which encapsulates the functionality of a Cholesky factorization for dense matrices.</para>
/// <para>For a symmetric, positive definite matrix A, the Cholesky decomposition /// <para>For a symmetric, positive definite matrix A, the Cholesky factorization
/// is an lower triangular matrix L so that A = L*L'.</para> /// is an lower triangular matrix L so that A = L*L'.</para>
/// </summary> /// </summary>
/// <remarks> /// <remarks>
/// 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. /// or positive definite, the constructor will throw an exception.
/// </remarks> /// </remarks>
public class DenseCholesky : Cholesky public class DenseCholesky : Cholesky
{ {
/// <summary>
/// Stores the Cholesky factor for the decomposition.
/// </summary>
private DenseMatrix _factor;
/// <summary> /// <summary>
/// Initializes a new instance of the <see cref="DenseCholesky"/> class. This object will compute the /// Initializes a new instance of the <see cref="DenseCholesky"/> 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.
/// </summary> /// </summary>
/// <param name="matrix">The matrix to factor.</param> /// <param name="matrix">The matrix to factor.</param>
/// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <b>null</b>.</exception> /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <b>null</b>.</exception>
@ -70,48 +65,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Decomposition
} }
// Create a new matrix for the Cholesky factor, then perform factorization (while overwriting). // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting).
_factor = (DenseMatrix)matrix.Clone(); var factor = (DenseMatrix)matrix.Clone();
Control.LinearAlgebraProvider.CholeskyFactor(_factor.Data, _factor.RowCount); Control.LinearAlgebraProvider.CholeskyFactor(factor.Data, factor.RowCount);
} mFactor = factor;
/// <summary>
/// Returns the lower triangular form of the Cholesky matrix.
/// </summary>
public override Matrix Factor
{
get { return _factor; }
}
/// <summary>
/// The determinant of the matrix for which the Cholesky matrix was computed.
/// </summary>
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;
}
}
/// <summary>
/// The log determinant of the matrix for which the Cholesky matrix was computed.
/// </summary>
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;
}
} }
} }
} }

6
src/Numerics/LinearAlgebra/Double/Factorization/ExtensionMethods.cs

@ -28,13 +28,13 @@
// OTHER DEALINGS IN THE SOFTWARE. // OTHER DEALINGS IN THE SOFTWARE.
// </copyright> // </copyright>
namespace MathNet.Numerics.LinearAlgebra.Double.Decomposition namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
{ {
using System; using System;
using Properties; using Properties;
/// <summary> /// <summary>
/// Extension methods which return decompositions for the various matrix classes. /// Extension methods which return factorizations for the various matrix classes.
/// </summary> /// </summary>
public static class ExtensionMethods public static class ExtensionMethods
{ {
@ -45,7 +45,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Decomposition
/// <returns>The Cholesky decomposition object.</returns> /// <returns>The Cholesky decomposition object.</returns>
public static Cholesky Cholesky(this Matrix matrix) public static Cholesky Cholesky(this Matrix matrix)
{ {
return Decomposition.Cholesky.Create(matrix); return Factorization.Cholesky.Create(matrix);
} }
} }
} }

17
src/Numerics/LinearAlgebra/Double/Matrix.cs

@ -639,5 +639,22 @@ namespace MathNet.Numerics.LinearAlgebra.Double
} }
} }
} }
/// <summary>
/// Returns the transpose of this matrix.
/// </summary>
/// <returns>The transpose of this matrix.</returns>
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;
}
} }
} }

59
src/UnitTests/LinearAlgebraTests/Double/Factorization/CholeskyTests.cs

@ -28,12 +28,12 @@
// OTHER DEALINGS IN THE SOFTWARE. // OTHER DEALINGS IN THE SOFTWARE.
// </copyright> // </copyright>
namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Decomposition namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
{ {
using System.Collections.Generic; using System.Collections.Generic;
using MbUnit.Framework; using MbUnit.Framework;
using LinearAlgebra.Double; using LinearAlgebra.Double;
using LinearAlgebra.Double.Decomposition; using LinearAlgebra.Double.Factorization;
public class CholeskyTests public class CholeskyTests
{ {
@ -41,8 +41,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Decomposition
[Row(1)] [Row(1)]
[Row(10)] [Row(10)]
[Row(100)] [Row(100)]
[Row(1000)] public void CanFactorizeIdentity(int order)
public void CanDecomposeIdentity(int order)
{ {
var I = DenseMatrix.Identity(order); var I = DenseMatrix.Identity(order);
var C = I.Cholesky(); var C = I.Cholesky();
@ -86,7 +85,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Decomposition
[Row(1)] [Row(1)]
[Row(10)] [Row(10)]
[Row(100)] [Row(100)]
[Row(1000)]
public void IdentityDeterminantIsOne(int order) public void IdentityDeterminantIsOne(int order)
{ {
var I = DenseMatrix.Identity(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(1.0, C.Determinant);
Assert.AreEqual(0.0, C.DeterminantLn); 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);
}
}
}
} }
} }

24
src/UnitTests/LinearAlgebraTests/Double/MatrixTests.cs

@ -457,5 +457,29 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
Vector result = CreateVector(matrix.RowCount - 1); Vector result = CreateVector(matrix.RowCount - 1);
matrix.GetColumn(0, 0, matrix.RowCount, result); 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]);
}
}
}
} }
} }
Loading…
Cancel
Save