Menu
▾
▴
[Dnanalytics-commits] r181 - in trunk: LinearAlgebra LinearAlgebra/Decomposition UnitTests UnitTests/LinearAlgebra/Decomposition
|
From: <in...@dn...> - 2006-01-20 15:52:20
|
Author: marcus
Date: 2006-01-20 10:52:02 -0500 (Fri, 20 Jan 2006)
New Revision: 181
Added:
trunk/UnitTests/LinearAlgebra/Decomposition/AbstractQRTest.cs
trunk/UnitTests/LinearAlgebra/Decomposition/AbstractSvdTest.cs
trunk/UnitTests/LinearAlgebra/Decomposition/DenseQRTest.cs
trunk/UnitTests/LinearAlgebra/Decomposition/DenseSvdTest.cs
Modified:
trunk/LinearAlgebra/Decomposition/DenseSvd.cs
trunk/LinearAlgebra/Decomposition/ISvd.cs
trunk/LinearAlgebra/Decomposition/Svd.cs
trunk/LinearAlgebra/ISimplifiedLapack.cs
trunk/LinearAlgebra/ManagedLapack.cs
trunk/LinearAlgebra/NativeLapack.cs
trunk/UnitTests/UnitTests.csproj
Log:
started filling out Svd.
Modified: trunk/LinearAlgebra/Decomposition/DenseSvd.cs
===================================================================
--- trunk/LinearAlgebra/Decomposition/DenseSvd.cs 2006-01-20 11:10:32 UTC (rev 180)
+++ trunk/LinearAlgebra/Decomposition/DenseSvd.cs 2006-01-20 15:52:02 UTC (rev 181)
@@ -17,10 +17,17 @@
/// <summary>
/// Class for computing the Singular Value Decomposition of a <see cref="Matrix"/>.
/// </summary>
- public sealed class DenseSvd : Algorithm
+ public sealed class DenseSvd : Algorithm, ISvd
{
#region Fields
+ private readonly bool computeVectors;
+ private readonly int rows;
+ private readonly int columns;
private DenseMatrix matrix;
+ private DenseMatrix u;
+ private DenseMatrix v;
+ private DenseVector s;
+ private int rank;
#endregion Fields
#region Constructors
@@ -28,8 +35,9 @@
/// Constructs an QR object for the given matrix.
/// </summary>
/// <param name="matrix">The matrix to factor.</param>
+ /// <param name="computeVectors">Whether to compute the left and right signular vectors or not.</param>
/// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
- public DenseSvd(Matrix matrix)
+ public DenseSvd(Matrix matrix, bool computeVectors)
{
if (matrix == null)
{
@@ -37,6 +45,9 @@
}
this.matrix = new DenseMatrix(matrix);
+ this.computeVectors = computeVectors;
+ this.rows = matrix.Rows;
+ this.columns = matrix.Columns;
}
#endregion Constructors
@@ -50,18 +61,19 @@
get
{
Compute();
- throw new NotImplementedException();
+ return s[0];
}
}
///<summary>Returns the condition number <c>max(S) / min(S)</c>.</summary>
///<value>The condition number.</value>
- public double Condition
+ public double ConditionNumber
{
get
{
Compute();
- throw new NotImplementedException();
+ int tmp = System.Math.Min(rows, columns) - 1;
+ return s[0] / s[tmp];
}
}
@@ -72,7 +84,7 @@
get
{
Compute();
- throw new NotImplementedException();
+ return rank;
}
}
#endregion Properties
@@ -84,8 +96,15 @@
///if <c>computeVectors</c> in the constructor is set to false.</returns>
public Matrix U()
{
- Compute();
- throw new NotImplementedException();
+ if (computeVectors)
+ {
+ Compute();
+ return u.Clone();
+ }
+ else
+ {
+ return null;
+ }
}
/// <summary>
@@ -97,16 +116,19 @@
/// <remarks>The method does nothing if <c>computeVectors</c> in the constructor is set to false.</remarks>
public void U(Matrix result)
{
- if (result == null)
+ if (computeVectors)
{
- throw new ArgumentNullException("result", Strings.NullParameterException);
+ if (result == null)
+ {
+ throw new ArgumentNullException("result", Strings.NullParameterException);
+ }
+ if (result.Rows != rows || result.Columns != rows)
+ {
+ throw new NotConformableException("result", Strings.ParameterNotConformable);
+ }
+ Compute();
+ u.CopyTo(result);
}
- if (result.Rows != matrix.Rows || result.Columns != matrix.Columns)
- {
- throw new NotConformableException("result", Strings.ParameterNotConformable);
- }
- Compute();
- throw new NotImplementedException();
}
///<summary>Returns the right singular vectors as a <see cref="Matrix"/>.</summary>
@@ -114,8 +136,16 @@
///if <c>computeVectors</c> in the constructor is set to false.</returns>
public Matrix V()
{
- Compute();
- throw new NotImplementedException();
+ if (computeVectors)
+ {
+ Compute();
+ return v.Clone();
+ }
+ else
+ {
+ return null;
+ }
+
}
/// <summary>
@@ -127,16 +157,19 @@
/// <remarks>The method does nothing if <c>computeVectors</c> in the constructor is set to false.</remarks>
public void V(Matrix result)
{
- if (result == null)
+ if (computeVectors)
{
- throw new ArgumentNullException("result", Strings.NullParameterException);
+ if (result == null)
+ {
+ throw new ArgumentNullException("result", Strings.NullParameterException);
+ }
+ if (result.Rows != columns || result.Columns != columns)
+ {
+ throw new NotConformableException("result", Strings.ParameterNotConformable);
+ }
+ Compute();
+ v.CopyTo(result);
}
- if (result.Rows != matrix.Rows || result.Columns != matrix.Columns)
- {
- throw new NotConformableException("result", Strings.ParameterNotConformable);
- }
- Compute();
- throw new NotImplementedException();
}
///<summary>Returns the singular values as a diagonal <see cref="Matrix"/>.</summary>
@@ -144,7 +177,9 @@
public Matrix W()
{
Compute();
- throw new NotImplementedException();
+ Matrix result = MatrixBuilder.CreateMatrix(rows, columns);
+ W(result);
+ return result;
}
/// <summary>
@@ -159,12 +194,21 @@
{
throw new ArgumentNullException("result", Strings.NullParameterException);
}
- if (result.Rows != matrix.Rows || result.Columns != matrix.Columns)
+ if (result.Rows != rows || result.Columns != columns)
{
throw new NotConformableException("result", Strings.ParameterNotConformable);
}
Compute();
- throw new NotImplementedException();
+ for (int i = 0; i < rows; i++)
+ {
+ for (int j = 0; j < columns; j++)
+ {
+ if (i == j)
+ {
+ result.ValueAt(i, i, s[i]);
+ }
+ }
+ }
}
///<summary>Returns the singular values as a <see cref="Vector"/>.</summary>
@@ -172,7 +216,7 @@
public Vector S()
{
Compute();
- throw new NotImplementedException();
+ return s.Clone();
}
/// <summary>
@@ -187,13 +231,18 @@
{
throw new ArgumentNullException("result", Strings.NullParameterException);
}
- //if (result.Rows != rows || result.Columns != rows)
- //{
- // throw new NotConformableException("result", Strings.ParameterNotConformable);
- // }
+ if (result.Count != s.Count)
+ {
+ throw new NotConformableException("result", Strings.ParameterNotConformable);
+ }
+
Compute();
- throw new NotImplementedException();
+ for (int i = 0; i < s.Count; i++)
+ {
+ result[i] = s[i];
+ }
+
}
#endregion Public Methods
@@ -204,6 +253,35 @@
protected override void InternalCompute()
{
+ int nm = System.Math.Min(rows, columns);
+ s = (DenseVector)VectorBuilder.CreateVector(nm);
+
+ if (computeVectors)
+ {
+ u = (DenseMatrix)MatrixBuilder.CreateMatrix(rows);
+ v = (DenseMatrix)MatrixBuilder.CreateMatrix(columns);
+ matrix.lapack.SVDDecomposition(rows, columns, matrix.data, s.data, u.data, v.data);
+ v.Transpose();
+ }
+ else
+ {
+ matrix.lapack.SVDDecomposition(rows, columns, matrix.data, s.data, null, null);
+ }
+
+
+ double eps = System.Math.Pow(2.0, -52.0);
+ double tol = System.Math.Max(rows, columns) * s[0] * eps;
+ rank = 0;
+
+ for (int h = 0; h < nm; h++)
+ {
+ if (s[h] > tol)
+ {
+ rank++;
+ }
+ }
+
+ matrix = null;
}
#endregion Protected Members
}
Modified: trunk/LinearAlgebra/Decomposition/ISvd.cs
===================================================================
--- trunk/LinearAlgebra/Decomposition/ISvd.cs 2006-01-20 11:10:32 UTC (rev 180)
+++ trunk/LinearAlgebra/Decomposition/ISvd.cs 2006-01-20 15:52:02 UTC (rev 181)
@@ -76,7 +76,7 @@
///<summary>Returns the condition number <c>max(S) / min(S)</c>.</summary>
///<value>The condition number.</value>
- double Condition { get; }
+ double ConditionNumber { get; }
///<summary>Returns the effective numerical matrix rank.</summary>
///<value>The number of non-negligible singular values.</value>
Modified: trunk/LinearAlgebra/Decomposition/Svd.cs
===================================================================
--- trunk/LinearAlgebra/Decomposition/Svd.cs 2006-01-20 11:10:32 UTC (rev 180)
+++ trunk/LinearAlgebra/Decomposition/Svd.cs 2006-01-20 15:52:02 UTC (rev 181)
@@ -31,15 +31,16 @@
/// Constructs an QR object for the given matrix.
/// </summary>
/// <param name="matrix">The matrix to factor.</param>
+ /// <param name="computeVectors">Whether to compute the left and right signular vectors or not.</param>
/// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
- public Svd(Matrix matrix)
+ public Svd(Matrix matrix, bool computeVectors)
{
if (matrix == null)
{
throw new ArgumentNullException("matrix", Strings.NullParameterException);
}
- //decomp = new DenseSvd(matrix);
+ decomp = new DenseSvd(matrix, computeVectors);
rows = matrix.Rows;
columns = matrix.Columns;
}
@@ -60,11 +61,11 @@
///<summary>Returns the condition number <c>max(S) / min(S)</c>.</summary>
///<value>The condition number.</value>
- public double Condition
+ public double ConditionNumber
{
get
{
- return decomp.Condition;
+ return decomp.ConditionNumber;
}
}
@@ -102,7 +103,7 @@
{
throw new ArgumentNullException("result", Strings.NullParameterException);
}
- if (result.Rows != rows || result.Columns != columns)
+ if (result.Rows != rows || result.Columns != rows)
{
throw new NotConformableException("result", Strings.ParameterNotConformable);
}
@@ -130,7 +131,7 @@
{
throw new ArgumentNullException("result", Strings.NullParameterException);
}
- if (result.Rows != rows || result.Columns != columns)
+ if (result.Rows != columns || result.Columns != columns)
{
throw new NotConformableException("result", Strings.ParameterNotConformable);
}
@@ -182,13 +183,14 @@
{
throw new ArgumentNullException("result", Strings.NullParameterException);
}
- //if (result.Rows != rows || result.Columns != rows)
- //{
- // throw new NotConformableException("result", Strings.ParameterNotConformable);
- // }
+ int nm = System.Math.Min(rows, columns);
+ if (result.Count != nm)
+ {
+ throw new NotConformableException("result", Strings.ParameterNotConformable);
+ }
+
decomp.S(result);
- throw new NotImplementedException();
}
#endregion Public Methods
Modified: trunk/LinearAlgebra/ISimplifiedLapack.cs
===================================================================
--- trunk/LinearAlgebra/ISimplifiedLapack.cs 2006-01-20 11:10:32 UTC (rev 180)
+++ trunk/LinearAlgebra/ISimplifiedLapack.cs 2006-01-20 15:52:02 UTC (rev 181)
@@ -40,5 +40,7 @@
int CholeskyFactor(int order, double[] data);
void QRFactor(int m, int n, double[] q, double[] r);
+
+ void SVDDecomposition(int rows, int cols, double[] a, double[] s, double[] u, double [] v);
}
}
\ No newline at end of file
Modified: trunk/LinearAlgebra/ManagedLapack.cs
===================================================================
--- trunk/LinearAlgebra/ManagedLapack.cs 2006-01-20 11:10:32 UTC (rev 180)
+++ trunk/LinearAlgebra/ManagedLapack.cs 2006-01-20 15:52:02 UTC (rev 181)
@@ -312,6 +312,18 @@
}
}
+ public void SVDDecomposition(int m, int n, double[] a, double[] s, double[] u, double[] v)
+ {
+ char jobu = 'A';
+ char jobvt = 'A';
+ if (u == null)
+ {
+ jobu = 'N';
+ jobvt = 'N';
+ }
+ // lapack.Dgesvd(jobu, jobvt, m, n, a, m, s, u, vt, n);
+ }
+
#endregion Public Methods
}
}
\ No newline at end of file
Modified: trunk/LinearAlgebra/NativeLapack.cs
===================================================================
--- trunk/LinearAlgebra/NativeLapack.cs 2006-01-20 11:10:32 UTC (rev 180)
+++ trunk/LinearAlgebra/NativeLapack.cs 2006-01-20 15:52:02 UTC (rev 181)
@@ -123,6 +123,19 @@
lapack.Dorgqr(m, m, n, q, m, tau);
}
}
+
+ public void SVDDecomposition(int m, int n, double[] a, double[] s, double[] u, double[] v)
+ {
+ char jobu = 'A';
+ char jobvt = 'A';
+ if (u == null)
+ {
+ jobu = 'N';
+ jobvt = 'N';
+ }
+ lapack.Dgesvd(jobu, jobvt, m, n, a, m, s, u, m, v, n);
+ }
+
#endregion Public Methods
}
}
\ No newline at end of file
Added: trunk/UnitTests/LinearAlgebra/Decomposition/AbstractQRTest.cs
===================================================================
--- trunk/UnitTests/LinearAlgebra/Decomposition/AbstractQRTest.cs 2006-01-20 11:10:32 UTC (rev 180)
+++ trunk/UnitTests/LinearAlgebra/Decomposition/AbstractQRTest.cs 2006-01-20 15:52:02 UTC (rev 181)
@@ -0,0 +1,319 @@
+using System;
+using System.Collections;
+using System.Collections.Generic;
+using NUnit.Framework;
+using dnAnalytics.Math;
+using dnAnalytics.LinearAlgebra;
+using dnAnalytics.LinearAlgebra.Decomposition;
+using dnAnalytics.LinearAlgebra.IO;
+using dnAnalytics.Exceptions;
+
+namespace dnAnalytics.UnitTests.LinearAlgebra.Decomposition
+{
+ public abstract class AbstractQRTest
+ {
+ public abstract void GetDecompositionClass(String file, out Matrix matrix, out QR qr);
+ private QR squareQr;
+ private QR singularQr;
+ private QR wideQr;
+ private QR tallQr;
+
+ private Matrix squareMatrix;
+ private Matrix singularMatrix;
+ private Matrix wideMatrix;
+ private Matrix tallMatrix;
+
+ public abstract Matrix GetMatrix(int order);
+ public abstract Matrix GetMatrix(int rows, int columns);
+
+ [TestFixtureSetUp]
+ public void SetUp()
+ {
+ MatrixMarketReader mmr = new MatrixMarketReader();
+ GetDecompositionClass("..\\..\\TestMatrices\\random_real_general_array_100.matrix_market", out squareMatrix, out squareQr);
+ GetDecompositionClass("..\\..\\TestMatrices\\gear_integer_general_coordinate_100.matrix_market", out singularMatrix, out singularQr);
+ GetDecompositionClass("..\\..\\TestMatrices\\random_real_general_array_10_20.matrix_market", out wideMatrix, out wideQr);
+ GetDecompositionClass("..\\..\\TestMatrices\\random_real_general_array_20_10.matrix_market", out tallMatrix, out tallQr);
+ }
+
+ [Test]
+ [ExpectedException(typeof(ArgumentNullException))]
+ public void QRConstructorNull()
+ {
+ Matrix matrix = null;
+ QR lu = new QR(matrix);
+ }
+
+ [Test]
+ public void NotFullRank()
+ {
+ Assert.IsFalse(singularQr.IsFullRank);
+ }
+
+ [Test]
+ public void FullRank()
+ {
+ Assert.IsTrue(squareQr.IsFullRank);
+ Assert.IsTrue(wideQr.IsFullRank);
+ Assert.IsTrue(tallQr.IsFullRank);
+ }
+
+ [Test]
+ public void Factors()
+ {
+ Matrix q = squareQr.Q();
+ Matrix r = squareQr.R();
+
+ Matrix result = GetMatrix(squareMatrix.Rows, squareMatrix.Columns);
+ q.Multiply(r, result);
+ double error = 0;
+ for (int i = 0; i < result.Rows; i++)
+ {
+ for (int j = 0; j < result.Columns; j++)
+ {
+ error = 0;
+ if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
+ {
+ error = System.Math.Abs(result[i, j] - squareMatrix[i, j]);
+ }
+ else
+ {
+ error = System.Math.Abs((result[i, j] - squareMatrix[i, j]) / squareMatrix[i, j]);
+ }
+ Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
+ }
+ }
+
+
+ q = tallQr.Q();
+ r = tallQr.R();
+
+ result = GetMatrix(tallMatrix.Rows, tallMatrix.Columns);
+ q.Multiply(r, result);
+
+ error = 0;
+ for (int i = 0; i < result.Rows; i++)
+ {
+ for (int j = 0; j < result.Columns; j++)
+ {
+ error = 0;
+ if (System.Math.Abs(tallMatrix[i, j]) < 1e10)
+ {
+ error = System.Math.Abs(result[i, j] - tallMatrix[i, j]);
+ }
+ else
+ {
+ error = System.Math.Abs((result[i, j] - tallMatrix[i, j]) / tallMatrix[i, j]);
+ }
+ Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
+ }
+ }
+
+ q = wideQr.Q();
+ r = wideQr.R();
+
+ result = GetMatrix(wideMatrix.Rows, wideMatrix.Columns);
+ q.Multiply(r, result);
+ error = 0;
+ for (int i = 0; i < result.Rows; i++)
+ {
+ for (int j = 0; j < result.Columns; j++)
+ {
+ error = 0;
+ if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
+ {
+ error = System.Math.Abs(result[i, j] - wideMatrix[i, j]);
+ }
+ else
+ {
+ error = System.Math.Abs((result[i, j] - wideMatrix[i, j]) / wideMatrix[i, j]);
+ }
+ Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
+ }
+ }
+
+ q = singularQr.Q();
+ r = singularQr.R();
+ result = GetMatrix(singularMatrix.Rows, singularMatrix.Columns);
+ q.Multiply(r, result);
+ error = 0;
+ for (int i = 0; i < result.Rows; i++)
+ {
+ for (int j = 0; j < result.Columns; j++)
+ {
+ error = 0;
+ if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
+ {
+ error = System.Math.Abs(result[i, j] - singularMatrix[i, j]);
+ }
+ else
+ {
+ error = System.Math.Abs((result[i, j] - singularMatrix[i, j]) / singularMatrix[i, j]);
+ }
+ Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
+ }
+ }
+ }
+
+ [Test]
+ public void FactorsResults()
+ {
+ Matrix q = GetMatrix(squareMatrix.Rows, squareMatrix.Rows);
+ squareQr.Q(q);
+ Matrix r = GetMatrix(squareMatrix.Rows, squareMatrix.Columns);
+ squareQr.R(r);
+
+ Matrix result = GetMatrix(squareMatrix.Rows, squareMatrix.Columns);
+ q.Multiply(r, result);
+ double error = 0;
+ for (int i = 0; i < result.Rows; i++)
+ {
+ for (int j = 0; j < result.Columns; j++)
+ {
+ error = 0;
+ if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
+ {
+ error = System.Math.Abs(result[i, j] - squareMatrix[i, j]);
+ }
+ else
+ {
+ error = System.Math.Abs((result[i, j] - squareMatrix[i, j]) / squareMatrix[i, j]);
+ }
+ Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
+ }
+ }
+
+
+ q = GetMatrix(tallMatrix.Rows, tallMatrix.Rows);
+ tallQr.Q(q);
+ r = GetMatrix(tallMatrix.Rows, tallMatrix.Columns);
+ tallQr.R(r);
+
+ result = GetMatrix(tallMatrix.Rows, tallMatrix.Columns);
+ q.Multiply(r, result);
+ error = 0;
+ for (int i = 0; i < result.Rows; i++)
+ {
+ for (int j = 0; j < result.Columns; j++)
+ {
+ error = 0;
+ if (System.Math.Abs(tallMatrix[i, j]) < 1e10)
+ {
+ error = System.Math.Abs(result[i, j] - tallMatrix[i, j]);
+ }
+ else
+ {
+ error = System.Math.Abs((result[i, j] - tallMatrix[i, j]) / tallMatrix[i, j]);
+ }
+ Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
+ }
+ }
+
+ q = GetMatrix(wideMatrix.Rows, wideMatrix.Rows);
+ wideQr.Q(q);
+ r = GetMatrix(wideMatrix.Rows, wideMatrix.Columns);
+ wideQr.R(r);
+
+ result = GetMatrix(wideMatrix.Rows, wideMatrix.Columns);
+ q.Multiply(r, result);
+ error = 0;
+ for (int i = 0; i < result.Rows; i++)
+ {
+ for (int j = 0; j < result.Columns; j++)
+ {
+ error = 0;
+ if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
+ {
+ error = System.Math.Abs(result[i, j] - wideMatrix[i, j]);
+ }
+ else
+ {
+ error = System.Math.Abs((result[i, j] - wideMatrix[i, j]) / wideMatrix[i, j]);
+ }
+ Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
+ }
+ }
+
+ q = singularQr.Q();
+ r = singularQr.R();
+ result = GetMatrix(singularMatrix.Rows, singularMatrix.Columns);
+ q.Multiply(r, result);
+ error = 0;
+ for (int i = 0; i < result.Rows; i++)
+ {
+ for (int j = 0; j < result.Columns; j++)
+ {
+ error = 0;
+ if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
+ {
+ error = System.Math.Abs(result[i, j] - singularMatrix[i, j]);
+ }
+ else
+ {
+ error = System.Math.Abs((result[i, j] - singularMatrix[i, j]) / singularMatrix[i, j]);
+ }
+ Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
+ }
+ }
+ }
+
+ [Test]
+ [ExpectedException(typeof(ArgumentNullException))]
+ public void QResultNull()
+ {
+ Matrix result = null;
+ squareQr.Q(result);
+ }
+
+ [Test]
+ [ExpectedException(typeof(ArgumentNullException))]
+ public void RResultNull()
+ {
+ Matrix result = null;
+ squareQr.R(result);
+ }
+
+ [Test]
+ [ExpectedException(typeof(NotConformableException))]
+ public void QResultNotConformable1()
+ {
+ Matrix result = GetMatrix(squareMatrix.Rows + 1, squareMatrix.Rows);
+ squareQr.Q(result);
+ }
+
+ [Test]
+ [ExpectedException(typeof(NotConformableException))]
+ public void QResultNotConformable2()
+ {
+ Matrix result = GetMatrix(squareMatrix.Rows, squareMatrix.Rows + 1);
+ squareQr.Q(result);
+ }
+
+ [Test]
+ [ExpectedException(typeof(NotConformableException))]
+ public void RResultNotConformable1()
+ {
+ Matrix result = GetMatrix(squareMatrix.Rows + 1, squareMatrix.Rows);
+ squareQr.R(result);
+ }
+
+ [Test]
+ [ExpectedException(typeof(NotConformableException))]
+ public void RResultNotConformable2()
+ {
+ Matrix result = GetMatrix(squareMatrix.Rows, squareMatrix.Rows + 1);
+ squareQr.R(result);
+ }
+
+ [Test]
+ public void Determinant()
+ {
+ Assert.AreEqual(0, singularQr.Determinant);
+ double expected = 5.006656015927541e248;
+ double actual = squareQr.Determinant;
+ double error = System.Math.Abs((actual - expected) / expected);
+ Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
+ }
+
+ }
+}
Added: trunk/UnitTests/LinearAlgebra/Decomposition/AbstractSvdTest.cs
===================================================================
--- trunk/UnitTests/LinearAlgebra/Decomposition/AbstractSvdTest.cs 2006-01-20 11:10:32 UTC (rev 180)
+++ trunk/UnitTests/LinearAlgebra/Decomposition/AbstractSvdTest.cs 2006-01-20 15:52:02 UTC (rev 181)
@@ -0,0 +1,431 @@
+using System;
+using System.Collections;
+using System.Collections.Generic;
+using NUnit.Framework;
+using dnAnalytics.Math;
+using dnAnalytics.LinearAlgebra;
+using dnAnalytics.LinearAlgebra.Decomposition;
+using dnAnalytics.LinearAlgebra.IO;
+using dnAnalytics.Exceptions;
+
+namespace dnAnalytics.UnitTests.LinearAlgebra.Decomposition
+{
+ public abstract class AbstractSvdTest
+ {
+ private static double ACCEPTABLE_ERROR = 1e-11;
+ public abstract void GetDecompositionClass(String file, out Matrix matrix, out Svd svd);
+ private Svd squareSvd;
+ private Svd singularSvd;
+ private Svd wideSvd;
+ private Svd tallSvd;
+
+ private Matrix squareMatrix;
+ private Matrix singularMatrix;
+ private Matrix wideMatrix;
+ private Matrix tallMatrix;
+
+ public abstract Vector GetVector(int size);
+ public abstract Matrix GetMatrix(int order);
+ public abstract Matrix GetMatrix(int rows, int columns);
+
+ [TestFixtureSetUp]
+ public void SetUp()
+ {
+ Settings.UseNativeLibrary = true;
+ MatrixMarketReader mmr = new MatrixMarketReader();
+ GetDecompositionClass("..\\..\\TestMatrices\\random_real_general_array_100.matrix_market", out squareMatrix, out squareSvd);
+ GetDecompositionClass("..\\..\\TestMatrices\\gear_integer_general_coordinate_100.matrix_market", out singularMatrix, out singularSvd);
+ GetDecompositionClass("..\\..\\TestMatrices\\random_real_general_array_10_20.matrix_market", out wideMatrix, out wideSvd);
+ GetDecompositionClass("..\\..\\TestMatrices\\random_real_general_array_20_10.matrix_market", out tallMatrix, out tallSvd);
+ }
+
+ [Test]
+ [ExpectedException(typeof(ArgumentNullException))]
+ public void SvdConstructorNull()
+ {
+ Matrix matrix = null;
+ Svd lu = new Svd(matrix, false);
+ }
+
+ [Test]
+ public void Decomposition()
+ {
+ Matrix u = squareSvd.U();
+ Matrix w = squareSvd.W();
+ Matrix vt = squareSvd.V().GetTranspose();
+ Matrix result = u * w * vt;
+
+ //Console.WriteLine(result.ToString("0.#"));
+ //Console.WriteLine(squareMatrix.ToString("0.#"));
+ double error = 0;
+ for (int i = 0; i < result.Rows; i++)
+ {
+ for (int j = 0; j < result.Columns; j++)
+ {
+ error = 0;
+ if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
+ {
+ error = System.Math.Abs(result[i, j] - squareMatrix[i, j]);
+ }
+ else
+ {
+ error = System.Math.Abs((result[i, j] - squareMatrix[i, j]) / squareMatrix[i, j]);
+ }
+
+ Assert.IsTrue(error < ACCEPTABLE_ERROR);
+ }
+ }
+
+ u = tallSvd.U();
+ w = tallSvd.W();
+ vt = tallSvd.V().GetTranspose();
+ result = u * w * vt;
+
+ error = 0;
+ for (int i = 0; i < result.Rows; i++)
+ {
+ for (int j = 0; j < result.Columns; j++)
+ {
+ error = 0;
+ if (System.Math.Abs(tallMatrix[i, j]) < 1e10)
+ {
+ error = System.Math.Abs(result[i, j] - tallMatrix[i, j]);
+ }
+ else
+ {
+ error = System.Math.Abs((result[i, j] - tallMatrix[i, j]) / tallMatrix[i, j]);
+ }
+ Assert.IsTrue(error < ACCEPTABLE_ERROR);
+ }
+ }
+
+ u = wideSvd.U();
+ w = wideSvd.W();
+ vt = wideSvd.V().GetTranspose();
+ result = u * w * vt;
+
+ error = 0;
+ for (int i = 0; i < result.Rows; i++)
+ {
+ for (int j = 0; j < result.Columns; j++)
+ {
+ error = 0;
+ if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
+ {
+ error = System.Math.Abs(result[i, j] - wideMatrix[i, j]);
+ }
+ else
+ {
+ error = System.Math.Abs((result[i, j] - wideMatrix[i, j]) / wideMatrix[i, j]);
+ }
+ Assert.IsTrue(error < ACCEPTABLE_ERROR);
+ }
+ }
+
+ u = singularSvd.U();
+ w = singularSvd.W();
+ vt = singularSvd.V().GetTranspose();
+ result = u * w * vt;
+
+ error = 0;
+ for (int i = 0; i < result.Rows; i++)
+ {
+ for (int j = 0; j < result.Columns; j++)
+ {
+ error = 0;
+ if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
+ {
+ error = System.Math.Abs(result[i, j] - singularMatrix[i, j]);
+ }
+ else
+ {
+ error = System.Math.Abs((result[i, j] - singularMatrix[i, j]) / singularMatrix[i, j]);
+ }
+ Assert.IsTrue(error < ACCEPTABLE_ERROR);
+ }
+ }
+ }
+
+ [Test]
+ public void DecompositionResults()
+ {
+ Matrix u = GetMatrix(squareMatrix.Rows);
+ squareSvd.U(u);
+ Matrix w = GetMatrix(squareMatrix.Rows, squareMatrix.Columns);
+ squareSvd.W(w);
+
+ Matrix vt = GetMatrix(squareMatrix.Columns);
+ squareSvd.V(vt);
+ vt.Transpose();
+
+ Matrix result = u * w * vt;
+
+ double error = 0;
+ for (int i = 0; i < result.Rows; i++)
+ {
+ for (int j = 0; j < result.Columns; j++)
+ {
+ error = 0;
+ if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
+ {
+ error = System.Math.Abs(result[i, j] - squareMatrix[i, j]);
+ }
+ else
+ {
+ error = System.Math.Abs((result[i, j] - squareMatrix[i, j]) / squareMatrix[i, j]);
+ }
+
+ Assert.IsTrue(error < ACCEPTABLE_ERROR);
+ }
+ }
+
+ u = GetMatrix(tallMatrix.Rows);
+ tallSvd.U(u);
+ w = GetMatrix(tallMatrix.Rows, tallMatrix.Columns);
+ tallSvd.W(w);
+ vt = GetMatrix(tallMatrix.Columns);
+ tallSvd.V(vt);
+ vt.Transpose();
+
+ result = u * w * vt;
+
+ error = 0;
+ for (int i = 0; i < result.Rows; i++)
+ {
+ for (int j = 0; j < result.Columns; j++)
+ {
+ error = 0;
+ if (System.Math.Abs(tallMatrix[i, j]) < 1e10)
+ {
+ error = System.Math.Abs(result[i, j] - tallMatrix[i, j]);
+ }
+ else
+ {
+ error = System.Math.Abs((result[i, j] - tallMatrix[i, j]) / tallMatrix[i, j]);
+ }
+ Assert.IsTrue(error < ACCEPTABLE_ERROR);
+ }
+ }
+
+ u = GetMatrix(wideMatrix.Rows);
+ wideSvd.U(u);
+ w = GetMatrix(wideMatrix.Rows, wideMatrix.Columns);
+ wideSvd.W(w);
+ vt = GetMatrix(wideMatrix.Columns);
+ wideSvd.V(vt);
+ vt.Transpose();
+
+ result = u * w * vt;
+
+ error = 0;
+ for (int i = 0; i < result.Rows; i++)
+ {
+ for (int j = 0; j < result.Columns; j++)
+ {
+ error = 0;
+ if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
+ {
+ error = System.Math.Abs(result[i, j] - wideMatrix[i, j]);
+ }
+ else
+ {
+ error = System.Math.Abs((result[i, j] - wideMatrix[i, j]) / wideMatrix[i, j]);
+ }
+ Assert.IsTrue(error < ACCEPTABLE_ERROR);
+ }
+ }
+
+ u = GetMatrix(singularMatrix.Rows);
+ singularSvd.U(u);
+ w = GetMatrix(singularMatrix.Rows, singularMatrix.Columns);
+ singularSvd.W(w);
+ vt = GetMatrix(singularMatrix.Columns);
+ singularSvd.V(vt);
+ vt.Transpose();
+
+ result = u * w * vt;
+
+ error = 0;
+ for (int i = 0; i < result.Rows; i++)
+ {
+ for (int j = 0; j < result.Columns; j++)
+ {
+ error = 0;
+ if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
+ {
+ error = System.Math.Abs(result[i, j] - singularMatrix[i, j]);
+ }
+ else
+ {
+ error = System.Math.Abs((result[i, j] - singularMatrix[i, j]) / singularMatrix[i, j]);
+ }
+ Assert.IsTrue(error < ACCEPTABLE_ERROR);
+ }
+ }
+ }
+
+ [Test]
+ public void DontComputeVectors()
+ {
+ Svd svd = new Svd(squareMatrix, false);
+ Assert.IsNull(svd.U());
+ Assert.IsNull(svd.V());
+
+ Matrix result = GetMatrix(squareMatrix.Rows);
+ svd.U(result);
+ for (int i = 0; i < result.Rows; i++)
+ {
+ for (int j = 0; j < result.Columns; j++)
+ {
+ Assert.AreEqual(0, result.ValueAt(i, j));
+ }
+ }
+
+ svd.V(result);
+ for (int i = 0; i < result.Rows; i++)
+ {
+ for (int j = 0; j < result.Columns; j++)
+ {
+ Assert.AreEqual(0, result.ValueAt(i, j));
+ }
+ }
+ }
+
+ [Test]
+ [ExpectedException(typeof(ArgumentNullException))]
+ public void UResultNull()
+ {
+ Matrix result = null;
+ squareSvd.U(result);
+ }
+
+ [Test]
+ [ExpectedException(typeof(ArgumentNullException))]
+ public void VResultNull()
+ {
+ Matrix result = null;
+ squareSvd.V(result);
+ }
+
+ [Test]
+ [ExpectedException(typeof(ArgumentNullException))]
+ public void WResultNull()
+ {
+ Matrix result = null;
+ squareSvd.W(result);
+ }
+
+ [Test]
+ [ExpectedException(typeof(ArgumentNullException))]
+ public void SResultNull()
+ {
+ Vector result = null;
+ squareSvd.S(result);
+ }
+
+
+ [Test]
+ [ExpectedException(typeof(NotConformableException))]
+ public void UResultNotConformable1()
+ {
+ Matrix result = GetMatrix(squareMatrix.Rows + 1, squareMatrix.Rows);
+ squareSvd.U(result);
+ }
+
+ [Test]
+ [ExpectedException(typeof(NotConformableException))]
+ public void UResultNotConformable2()
+ {
+ Matrix result = GetMatrix(squareMatrix.Rows, squareMatrix.Rows + 1);
+ squareSvd.U(result);
+ }
+
+ [Test]
+ [ExpectedException(typeof(NotConformableException))]
+ public void VResultNotConformable1()
+ {
+ Matrix result = GetMatrix(squareMatrix.Columns + 1, squareMatrix.Columns);
+ squareSvd.V(result);
+ }
+
+ [Test]
+ [ExpectedException(typeof(NotConformableException))]
+ public void VResultNotConformable2()
+ {
+ Matrix result = GetMatrix(squareMatrix.Columns, squareMatrix.Columns + 1);
+ squareSvd.V(result);
+ }
+
+ [Test]
+ [ExpectedException(typeof(NotConformableException))]
+ public void WResultNotConformable1()
+ {
+ Matrix result = GetMatrix(squareMatrix.Rows + 1, squareMatrix.Columns);
+ squareSvd.W(result);
+ }
+
+ [Test]
+ [ExpectedException(typeof(NotConformableException))]
+ public void WResultNotConformable2()
+ {
+ Matrix result = GetMatrix(squareMatrix.Rows, squareMatrix.Columns + 1);
+ squareSvd.W(result);
+ }
+
+ [Test]
+ [ExpectedException(typeof(NotConformableException))]
+ public void SResultNotConformable1()
+ {
+ Vector result = GetVector(squareMatrix.Columns + 1);
+ squareSvd.S(result);
+ }
+
+ [Test]
+ public void Rank()
+ {
+ int rank = squareSvd.Rank;
+ Assert.AreEqual(squareMatrix.Rows, rank);
+
+ rank = tallSvd.Rank;
+ Assert.AreEqual(tallMatrix.Columns, rank);
+
+ rank = wideSvd.Rank;
+ Assert.AreEqual(wideMatrix.Rows, rank);
+
+ rank = singularSvd.Rank;
+ Assert.AreEqual(squareMatrix.Rows - 1, rank);
+ }
+
+ [Test]
+ public void Norm()
+ {
+ double norm = squareSvd.Norm2;
+ Assert.AreEqual(1011.4191, norm, 1e-4);
+
+ norm = tallSvd.Norm2;
+ // Assert.AreEqual(tallMatrix.Columns, norm);
+
+ norm = wideSvd.Norm2;
+ //Assert.AreEqual(wideMatrix.Rows, norm);
+
+ norm = singularSvd.Norm2;
+ Assert.AreEqual(2, norm);
+ }
+
+ [Test]
+ public void ConditionNumber()
+ {
+ double cn = squareSvd.ConditionNumber;
+ Assert.AreEqual(squareMatrix.Rows, cn);
+
+ cn = tallSvd.ConditionNumber;
+ Assert.AreEqual(tallMatrix.Columns, cn);
+
+ cn = wideSvd.ConditionNumber;
+ Assert.AreEqual(wideMatrix.Rows, cn);
+
+ cn = singularSvd.ConditionNumber;
+ Assert.AreEqual(squareMatrix.Rows - 1, cn);
+ }
+ }
+}
Added: trunk/UnitTests/LinearAlgebra/Decomposition/DenseQRTest.cs
===================================================================
--- trunk/UnitTests/LinearAlgebra/Decomposition/DenseQRTest.cs 2006-01-20 11:10:32 UTC (rev 180)
+++ trunk/UnitTests/LinearAlgebra/Decomposition/DenseQRTest.cs 2006-01-20 15:52:02 UTC (rev 181)
@@ -0,0 +1,33 @@
+using System;
+using System.Collections;
+using System.Collections.Generic;
+using NUnit.Framework;
+using dnAnalytics.Math;
+using dnAnalytics.LinearAlgebra;
+using dnAnalytics.LinearAlgebra.Decomposition;
+using dnAnalytics.LinearAlgebra.IO;
+using dnAnalytics.Exceptions;
+
+namespace dnAnalytics.UnitTests.LinearAlgebra.Decomposition
+{
+ [TestFixture]
+ public class DenseQRTest : AbstractQRTest
+ {
+ public override Matrix GetMatrix(int rows, int columns)
+ {
+ return MatrixBuilder.CreateMatrix(rows, columns, MatrixType.Dense);
+ }
+
+ public override Matrix GetMatrix(int order)
+ {
+ return MatrixBuilder.CreateMatrix(order, MatrixType.Dense);
+ }
+
+ public override void GetDecompositionClass(String file, out Matrix matrix, out QR qr)
+ {
+ MatrixMarketReader mmr = new MatrixMarketReader();
+ matrix = mmr.ReadMatrix(file, MatrixType.Dense);
+ qr = new QR(matrix);
+ }
+ }
+}
Added: trunk/UnitTests/LinearAlgebra/Decomposition/DenseSvdTest.cs
===================================================================
--- trunk/UnitTests/LinearAlgebra/Decomposition/DenseSvdTest.cs 2006-01-20 11:10:32 UTC (rev 180)
+++ trunk/UnitTests/LinearAlgebra/Decomposition/DenseSvdTest.cs 2006-01-20 15:52:02 UTC (rev 181)
@@ -0,0 +1,38 @@
+using System;
+using System.Collections;
+using System.Collections.Generic;
+using NUnit.Framework;
+using dnAnalytics.Math;
+using dnAnalytics.LinearAlgebra;
+using dnAnalytics.LinearAlgebra.Decomposition;
+using dnAnalytics.LinearAlgebra.IO;
+using dnAnalytics.Exceptions;
+
+namespace dnAnalytics.UnitTests.LinearAlgebra.Decomposition
+{
+ [TestFixture]
+ public class DenseSvdTest : AbstractSvdTest
+ {
+ public override Vector GetVector(int size)
+ {
+ return VectorBuilder.CreateVector(size, VectorType.Dense);
+ }
+
+ public override Matrix GetMatrix(int rows, int columns)
+ {
+ return MatrixBuilder.CreateMatrix(rows, columns, MatrixType.Dense);
+ }
+
+ public override Matrix GetMatrix(int order)
+ {
+ return MatrixBuilder.CreateMatrix(order, MatrixType.Dense);
+ }
+
+ public override void GetDecompositionClass(String file, out Matrix matrix, out Svd svd)
+ {
+ MatrixMarketReader mmr = new MatrixMarketReader();
+ matrix = mmr.ReadMatrix(file, MatrixType.Dense);
+ svd = new Svd(matrix, true);
+ }
+ }
+}
Modified: trunk/UnitTests/UnitTests.csproj
===================================================================
--- trunk/UnitTests/UnitTests.csproj 2006-01-20 11:10:32 UTC (rev 180)
+++ trunk/UnitTests/UnitTests.csproj 2006-01-20 15:52:02 UTC (rev 181)
@@ -52,6 +52,8 @@
<Compile Include="LinearAlgebra\Decomposition\AbstractLUTest.cs" />
<Compile Include="LinearAlgebra\Decomposition\AbstractCholeskyTest.cs" />
<Compile Include="LinearAlgebra\Decomposition\AbstractQRTest.cs" />
+ <Compile Include="LinearAlgebra\Decomposition\AbstractSvdTest.cs" />
+ <Compile Include="LinearAlgebra\Decomposition\DenseSvdTest.cs" />
<Compile Include="LinearAlgebra\Decomposition\DenseQRTest.cs" />
<Compile Include="LinearAlgebra\Decomposition\DenseCholeskyTest.cs" />
<Compile Include="LinearAlgebra\Decomposition\DenseLUTest.cs">
|
