Added Cholesky factorization.

src/Numerics/LinearAlgebra/Double/Decomposition/Cholesky.cs

@@ -0,0 +1,78 @@
+// <copyright file="Cholesky.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.LinearAlgebra.Double.Decomposition
+{
+    using System;
+    using Properties;
+
+    /// <summary>
+    /// <para>A class which encapsulates the functionality of a Cholesky decomposition.</para>
+    /// <para>For a symmetric, positive definite matrix A, the Cholesky decomposition
+    /// is an lower triangular matrix L so that A = L*L'.</para>
+    /// </summary>
+    /// <remarks>
+    /// The computation of the Cholesky decomposition is done at construction time. If the matrix is not symmetric
+    /// or positive definite, the constructor will throw an exception.
+    /// </remarks>
+    public abstract class Cholesky
+    {
+        /// <summary>
+        /// Internal method which routes the call to perform the Cholesky factorization to the appropriate class.
+        /// </summary>
+        /// <param name="matrix">The matrix to factor.</param>
+        /// <returns>A cholesky factorization object.</returns>
+        internal static Cholesky Create(Matrix matrix)
+        {
+            var dense = matrix as DenseMatrix;
+            if (dense != null)
+            {
+                return new DenseCholesky(dense);
+            }
+
+            throw new NotImplementedException();
+        }
+
+        /// <summary>
+        /// A reference to the lower triangular form of the Cholesky matrix.
+        /// </summary>
+        public abstract Matrix Factor { get; }
+
+        /// <summary>
+        /// The determinant of the matrix for which the Cholesky matrix was computed.
+        /// </summary>
+        public abstract double Determinant { get; }
+
+        /// <summary>
+        /// The log determinant of the matrix for which the Cholesky matrix was computed.
+        /// </summary>
+        public abstract double DeterminantLn { get; }
+    }
+}

src/Numerics/LinearAlgebra/Double/Decomposition/DenseCholesky.cs

@@ -0,0 +1,117 @@
+// <copyright file="DenseCholesky.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.LinearAlgebra.Double.Decomposition
+{
+    using System;
+    using Properties;
+
+    /// <summary>
+    /// <para>A class which encapsulates the functionality of a Cholesky decomposition for dense matrices.</para>
+    /// <para>For a symmetric, positive definite matrix A, the Cholesky decomposition
+    /// is an lower triangular matrix L so that A = L*L'.</para>
+    /// </summary>
+    /// <remarks>
+    /// The computation of the Cholesky decomposition is done at construction time. If the matrix is not symmetric
+    /// or positive definite, the constructor will throw an exception.
+    /// </remarks>
+    public class DenseCholesky : Cholesky
+    {
+        /// <summary>
+        /// Stores the Cholesky factor for the decomposition.
+        /// </summary>
+        private DenseMatrix _factor;
+
+        /// <summary>
+        /// 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.
+        /// </summary>
+        /// <param name="matrix">The matrix to factor.</param>
+        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <b>null</b>.</exception>
+        /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
+        /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not positive definite.</exception>
+        public DenseCholesky(DenseMatrix matrix)
+        {
+            if (matrix == null)
+            {
+                throw new ArgumentNullException("matrix");
+            }
+
+            if (matrix.RowCount != matrix.ColumnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSquare);
+            }
+
+            // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting).
+            _factor = (DenseMatrix)matrix.Clone();
+            Control.LinearAlgebraProvider.CholeskyFactor(_factor.Data, _factor.RowCount);
+        }
+
+        /// <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;
+            }
+        }
+    }
+}

src/Numerics/LinearAlgebra/Double/Decomposition/ExtensionMethods.cs

@@ -0,0 +1,51 @@
+// <copyright file="ExtensionMethods.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.LinearAlgebra.Double.Decomposition
+{
+    using System;
+    using Properties;
+
+    /// <summary>
+    /// Extension methods which return decompositions for the various matrix classes.
+    /// </summary>
+    public static class ExtensionMethods
+    {
+        /// <summary>
+        /// Computes the cholesky decomposition for a matrix.
+        /// </summary>
+        /// <param name="matrix">The matrix to factor.</param>
+        /// <returns>The Cholesky decomposition object.</returns>
+        public static Cholesky Cholesky(this Matrix matrix)
+        {
+            return Decomposition.Cholesky.Create(matrix);
+        }
+    }
+}

src/Numerics/Numerics.csproj

@@ -139,6 +139,9 @@
     <Compile Include="Interpolation\IInterpolation.cs" />
     <Compile Include="Interpolation\Interpolate.cs" />
     <Compile Include="Interpolation\SplineBoundaryCondition.cs" />
+    <Compile Include="LinearAlgebra\Double\Decomposition\Cholesky.cs" />
+    <Compile Include="LinearAlgebra\Double\Decomposition\ExtensionMethods.cs" />
+    <Compile Include="LinearAlgebra\Double\Decomposition\DenseCholesky.cs" />
     <Compile Include="LinearAlgebra\Double\Matrix.Arithmetic.cs" />
     <Compile Include="LinearAlgebra\Double\DenseMatrix.cs" />
     <Compile Include="LinearAlgebra\Double\DenseVector.cs" />
@@ -246,6 +249,7 @@
       <Install>true</Install>
     </BootstrapperPackage>
   </ItemGroup>
+  <ItemGroup />
   <Import Project="$(MSBuildToolsPath)\Microsoft.CSharp.targets" />
   <!-- To modify your build process, add your task inside one of the targets below and uncomment it. 
        Other similar extension points exist, see Microsoft.Common.targets.

src/UnitTests/LinearAlgebraTests/Double/Decomposition/CholeskyTests.cs

@@ -0,0 +1,98 @@
+// <copyright file="CholeskyTests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2010 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Decomposition
+{
+    using System.Collections.Generic;
+    using MbUnit.Framework;
+    using LinearAlgebra.Double;
+    using LinearAlgebra.Double.Decomposition;
+
+    public class CholeskyTests
+    {
+        [Test]
+        [Row(1)]
+        [Row(10)]
+        [Row(100)]
+        [Row(1000)]
+        public void CanDecomposeIdentity(int order)
+        {
+            var I = DenseMatrix.Identity(order);
+            var C = I.Cholesky();
+
+            for (var i = 0; i < C.Factor.RowCount; i++)
+            {
+                for (var j = 0; j < C.Factor.ColumnCount; j++)
+                {
+                    if (i == j)
+                    {
+                        Assert.AreEqual(1.0, C.Factor[i, j]);
+                    }
+                    else
+                    {
+                        Assert.AreEqual(0.0, C.Factor[i, j]);
+                    }
+                }
+            }
+        }
+
+        [Test]
+        [ExpectedArgumentException]
+        public void CholeskyFailsWithDiagonalNonPositiveDefiniteMatrix()
+        {
+            var I = DenseMatrix.Identity(10);
+            I[3, 3] = -4.0;
+            var C = I.Cholesky();
+        }
+
+        [Test]
+        [Row(3,5)]
+        [Row(5,3)]
+        [ExpectedArgumentException]
+        public void CholeskyFailsWithNonSquareMatrix(int row, int col)
+        {
+            var I = new DenseMatrix(row, col);
+            var C = I.Cholesky();
+        }
+
+        [Test]
+        [Row(1)]
+        [Row(10)]
+        [Row(100)]
+        [Row(1000)]
+        public void IdentityDeterminantIsOne(int order)
+        {
+            var I = DenseMatrix.Identity(order);
+            var C = I.Cholesky();
+            Assert.AreEqual(1.0, C.Determinant);
+            Assert.AreEqual(0.0, C.DeterminantLn);
+        }
+    }
+}

src/UnitTests/UnitTests.csproj

@@ -115,6 +115,7 @@
     <Compile Include="InterpolationTests\InterpolationFunctionalContract.cs" />
     <Compile Include="InterpolationTests\InterpolationInfrastructureContract.cs" />
     <Compile Include="InterpolationTests\InterpolationFunctionalTest.cs" />
+    <Compile Include="LinearAlgebraTests\Double\Decomposition\CholeskyTests.cs" />
     <Compile Include="LinearAlgebraTests\Double\MatrixLoader.cs" />
     <Compile Include="LinearAlgebraTests\Double\MatrixTests.Arithmetic.cs" />
     <Compile Include="LinearAlgebraTests\Double\LinearAlgebraProviderTests.cs" />