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