forked from tsai/mathnet-numerics
Marcus Cuda
16 years ago
16 changed files with 4116 additions and 678 deletions
File diff suppressed because it is too large
@ -0,0 +1,180 @@ |
|||||
|
// <copyright file="DenseSvd.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.Factorization |
||||
|
{ |
||||
|
using System; |
||||
|
using Properties; |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// <para>A class which encapsulates the functionality of the singular value decomposition (SVD) for <see cref="DenseMatrix"/>.</para>
|
||||
|
/// <para>Suppose M is an m-by-n matrix whose entries are real numbers.
|
||||
|
/// Then there exists a factorization of the form M = UΣVT where:
|
||||
|
/// - U is an m-by-m unitary matrix;
|
||||
|
/// - Σ is m-by-n diagonal matrix with nonnegative real numbers on the diagonal;
|
||||
|
/// - VT denotes transpose of V, an n-by-n unitary matrix;
|
||||
|
/// Such a factorization is called a singular-value decomposition of M. A common convention is to order the diagonal
|
||||
|
/// entries Σ(i,i) in descending order. In this case, the diagonal matrix Σ is uniquely determined
|
||||
|
/// by M (though the matrices U and V are not). The diagonal entries of Σ are known as the singular values of M.</para>
|
||||
|
/// </summary>
|
||||
|
/// <remarks>
|
||||
|
/// The computation of the singular value decomposition is done at construction time.
|
||||
|
/// </remarks>
|
||||
|
public class DenseSvd : Svd |
||||
|
{ |
||||
|
/// <summary>
|
||||
|
/// Initializes a new instance of the <see cref="DenseSvd"/> class. This object will compute the
|
||||
|
/// the singular value decomposition when the constructor is called and cache it's decomposition.
|
||||
|
/// </summary>
|
||||
|
/// <param name="matrix">The matrix to factor.</param>
|
||||
|
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
|
||||
|
/// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <b>null</b>.</exception>
|
||||
|
/// <exception cref="ArgumentException">If SVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
|
||||
|
public DenseSvd(DenseMatrix matrix, bool computeVectors) |
||||
|
{ |
||||
|
if (matrix == null) |
||||
|
{ |
||||
|
throw new ArgumentNullException("matrix"); |
||||
|
} |
||||
|
|
||||
|
ComputeVectors = computeVectors; |
||||
|
var nm = Math.Min(matrix.RowCount, matrix.ColumnCount); |
||||
|
VectorS = new DenseVector(nm); |
||||
|
MatrixU = new DenseMatrix(matrix.RowCount); |
||||
|
MatrixVT = new DenseMatrix(matrix.ColumnCount); |
||||
|
Control.LinearAlgebraProvider.SingularValueDecomposition(computeVectors, ((DenseMatrix)matrix.Clone()).Data, matrix.RowCount, matrix.ColumnCount, ((DenseVector)VectorS).Data, ((DenseMatrix)MatrixU).Data, ((DenseMatrix)MatrixVT).Data); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Solves a system of linear equations, <b>AX = B</b>, with A SVD factorized.
|
||||
|
/// </summary>
|
||||
|
/// <param name="input">The right hand side <see cref="Matrix"/>, <b>B</b>.</param>
|
||||
|
/// <param name="result">The left hand side <see cref="Matrix"/>, <b>X</b>.</param>
|
||||
|
public override void Solve(Matrix input, Matrix result) |
||||
|
{ |
||||
|
// Check for proper arguments.
|
||||
|
if (input == null) |
||||
|
{ |
||||
|
throw new ArgumentNullException("input"); |
||||
|
} |
||||
|
|
||||
|
if (result == null) |
||||
|
{ |
||||
|
throw new ArgumentNullException("result"); |
||||
|
} |
||||
|
|
||||
|
if (!ComputeVectors) |
||||
|
{ |
||||
|
throw new InvalidOperationException(Resources.SingularVectorsNotComputed); |
||||
|
} |
||||
|
|
||||
|
// The solution X should have the same number of columns as B
|
||||
|
if (input.ColumnCount != result.ColumnCount) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); |
||||
|
} |
||||
|
|
||||
|
// The dimension compatibility conditions for X = A\B require the two matrices A and B to have the same number of rows
|
||||
|
if (MatrixU.RowCount != input.RowCount) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); |
||||
|
} |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
if (MatrixVT.ColumnCount != result.RowCount) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); |
||||
|
} |
||||
|
|
||||
|
var dinput = input as DenseMatrix; |
||||
|
if (dinput == null) |
||||
|
{ |
||||
|
throw new NotImplementedException("Can only do SVD factorization for dense matrices at the moment."); |
||||
|
} |
||||
|
|
||||
|
var dresult = result as DenseMatrix; |
||||
|
if (dresult == null) |
||||
|
{ |
||||
|
throw new NotImplementedException("Can only do SVD factorization for dense matrices at the moment."); |
||||
|
} |
||||
|
|
||||
|
Control.LinearAlgebraProvider.SvdSolveFactored(MatrixU.RowCount, MatrixVT.ColumnCount, ((DenseVector)VectorS).Data, ((DenseMatrix)MatrixU).Data, ((DenseMatrix)MatrixVT).Data, dinput.Data, input.ColumnCount, dresult.Data); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Solves a system of linear equations, <b>Ax = b</b>, with A SVD factorized.
|
||||
|
/// </summary>
|
||||
|
/// <param name="input">The right hand side vector, <b>b</b>.</param>
|
||||
|
/// <param name="result">The left hand side <see cref="Matrix"/>, <b>x</b>.</param>
|
||||
|
public override void Solve(Vector input, Vector result) |
||||
|
{ |
||||
|
if (input == null) |
||||
|
{ |
||||
|
throw new ArgumentNullException("input"); |
||||
|
} |
||||
|
|
||||
|
if (result == null) |
||||
|
{ |
||||
|
throw new ArgumentNullException("result"); |
||||
|
} |
||||
|
|
||||
|
if (!ComputeVectors) |
||||
|
{ |
||||
|
throw new InvalidOperationException(Resources.SingularVectorsNotComputed); |
||||
|
} |
||||
|
|
||||
|
// Ax=b where A is an m x n matrix
|
||||
|
// Check that b is a column vector with m entries
|
||||
|
if (MatrixU.RowCount != input.Count) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentVectorsSameLength); |
||||
|
} |
||||
|
|
||||
|
// Check that x is a column vector with n entries
|
||||
|
if (MatrixVT.ColumnCount != result.Count) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixDimensions); |
||||
|
} |
||||
|
|
||||
|
var dinput = input as DenseVector; |
||||
|
if (dinput == null) |
||||
|
{ |
||||
|
throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment."); |
||||
|
} |
||||
|
|
||||
|
var dresult = result as DenseVector; |
||||
|
if (dresult == null) |
||||
|
{ |
||||
|
throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment."); |
||||
|
} |
||||
|
|
||||
|
Control.LinearAlgebraProvider.SvdSolveFactored(MatrixU.RowCount, MatrixVT.ColumnCount, ((DenseVector)VectorS).Data, ((DenseMatrix)MatrixU).Data, ((DenseMatrix)MatrixVT).Data, dinput.Data, 1, dresult.Data); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,283 @@ |
|||||
|
// <copyright file="Svd.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.Factorization |
||||
|
{ |
||||
|
using System; |
||||
|
using Properties; |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// <para>A class which encapsulates the functionality of the singular value decomposition (SVD).</para>
|
||||
|
/// <para>Suppose M is an m-by-n matrix whose entries are real numbers.
|
||||
|
/// Then there exists a factorization of the form M = UΣVT where:
|
||||
|
/// - U is an m-by-m unitary matrix;
|
||||
|
/// - Σ is m-by-n diagonal matrix with nonnegative real numbers on the diagonal;
|
||||
|
/// - VT denotes transpose of V, an n-by-n unitary matrix;
|
||||
|
/// Such a factorization is called a singular-value decomposition of M. A common convention is to order the diagonal
|
||||
|
/// entries Σ(i,i) in descending order. In this case, the diagonal matrix Σ is uniquely determined
|
||||
|
/// by M (though the matrices U and V are not). The diagonal entries of Σ are known as the singular values of M.</para>
|
||||
|
/// </summary>
|
||||
|
/// <remarks>
|
||||
|
/// The computation of the singular value decomposition is done at construction time.
|
||||
|
/// </remarks>
|
||||
|
public abstract class Svd : ISolver |
||||
|
{ |
||||
|
/// <summary>
|
||||
|
/// Gets or sets a value indicating whether to compute U and VT matrices during SVD factorization or not
|
||||
|
/// </summary>
|
||||
|
protected bool ComputeVectors |
||||
|
{ |
||||
|
get; |
||||
|
set; |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Gets or sets the singular values (Σ) of matrix in ascending value.
|
||||
|
/// </summary>
|
||||
|
protected Vector VectorS |
||||
|
{ |
||||
|
get; |
||||
|
set; |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Gets or sets left singular vectors (U - m-by-m unitary matrix)
|
||||
|
/// </summary>
|
||||
|
protected Matrix MatrixU |
||||
|
{ |
||||
|
get; |
||||
|
set; |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Gets or sets transpose right singular vectors (transpose of V, an n-by-n unitary matrix
|
||||
|
/// </summary>
|
||||
|
protected Matrix MatrixVT |
||||
|
{ |
||||
|
get; |
||||
|
set; |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Gets the effective numerical matrix rank.
|
||||
|
/// </summary>
|
||||
|
/// <value>The number of non-negligible singular values.</value>
|
||||
|
public virtual int Rank |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
var eps = Math.Pow(2.0, -52.0); |
||||
|
var tol = Math.Max(MatrixU.RowCount, MatrixVT.ColumnCount) * VectorS[0] * eps; |
||||
|
var nm = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount); |
||||
|
var rank = 0; |
||||
|
for (var h = 0; h < nm; h++) |
||||
|
{ |
||||
|
if (VectorS[h] > tol) |
||||
|
{ |
||||
|
rank++; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
return rank; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Internal method which routes the call to perform the singular value decomposition to the appropriate class.
|
||||
|
/// </summary>
|
||||
|
/// <param name="matrix">The matrix to factor.</param>
|
||||
|
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
|
||||
|
/// <returns>An SVD object.</returns>
|
||||
|
internal static Svd Create(Matrix matrix, bool computeVectors) |
||||
|
{ |
||||
|
var dense = matrix as DenseMatrix; |
||||
|
if (dense != null) |
||||
|
{ |
||||
|
return new DenseSvd(dense, computeVectors); |
||||
|
} |
||||
|
|
||||
|
throw new NotImplementedException(); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Gets the two norm of the <see cref="Matrix"/>.
|
||||
|
/// </summary>
|
||||
|
/// <returns>The 2-norm of the <see cref="Matrix"/>.</returns>
|
||||
|
public virtual double Norm2 |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
return VectorS[0]; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Gets the condition number <b>max(S) / min(S)</b>
|
||||
|
/// </summary>
|
||||
|
/// <returns>The condition number.</returns>
|
||||
|
public virtual double ConditionNumber |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
var tmp = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount) - 1; |
||||
|
return VectorS[0] / VectorS[tmp]; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Gets the determinant of the square matrix for which the SVD was computed.
|
||||
|
/// </summary>
|
||||
|
public virtual double Determinant |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
if (MatrixU.RowCount != MatrixVT.ColumnCount) |
||||
|
{ |
||||
|
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
||||
|
} |
||||
|
|
||||
|
var det = 1.0; |
||||
|
for (var i = 0; i < VectorS.Count; i++) |
||||
|
{ |
||||
|
det *= VectorS[i]; |
||||
|
if (Math.Abs(VectorS[i]).AlmostEqualInDecimalPlaces(0.0, 15)) |
||||
|
{ |
||||
|
return 0; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
return Math.Abs(det); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>Returns the left singular vectors as a <see cref="Matrix"/>.</summary>
|
||||
|
/// <returns>The left singular vectors. The matrix will be <c>null</c>, if <b>computeVectors</b> in the constructor is set to <c>false</c>.</returns>
|
||||
|
public Matrix U() |
||||
|
{ |
||||
|
return ComputeVectors ? MatrixU.Clone() : null; |
||||
|
} |
||||
|
|
||||
|
/// <summary>Returns the right singular vectors as a <see cref="Matrix"/>.</summary>
|
||||
|
/// <returns>The right singular vectors. The matrix will be <c>null</c>, if <b>computeVectors</b> in the constructor is set to <c>false</c>.</returns>
|
||||
|
/// <remarks>This is the transpose of the V matrix.</remarks>
|
||||
|
public Matrix VT() |
||||
|
{ |
||||
|
return ComputeVectors ? MatrixVT.Clone() : null; |
||||
|
} |
||||
|
|
||||
|
/// <summary>Returns the singular values as a diagonal <see cref="Matrix"/>.</summary>
|
||||
|
/// <returns>The singular values as a diagonal <see cref="Matrix"/>.</returns>
|
||||
|
public Matrix W() |
||||
|
{ |
||||
|
var rows = MatrixU.RowCount; |
||||
|
var columns = MatrixVT.ColumnCount; |
||||
|
var result = MatrixU.CreateMatrix(rows, columns); |
||||
|
for (var i = 0; i < rows; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < columns; j++) |
||||
|
{ |
||||
|
if (i == j) |
||||
|
{ |
||||
|
result.At(i, i, VectorS[i]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
return result; |
||||
|
} |
||||
|
|
||||
|
/// <summary>Returns the singular values as a <see cref="Vector"/>.</summary>
|
||||
|
/// <returns>the singular values as a <see cref="Vector"/>.</returns>
|
||||
|
public Vector S() |
||||
|
{ |
||||
|
return VectorS.Clone(); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Solves a system of linear equations, <b>AX = B</b>, with A SVD factorized.
|
||||
|
/// </summary>
|
||||
|
/// <param name="input">The right hand side <see cref="Matrix"/>, <b>B</b>.</param>
|
||||
|
/// <returns>The left hand side <see cref="Matrix"/>, <b>X</b>.</returns>
|
||||
|
public virtual Matrix Solve(Matrix input) |
||||
|
{ |
||||
|
// Check for proper arguments.
|
||||
|
if (input == null) |
||||
|
{ |
||||
|
throw new ArgumentNullException("input"); |
||||
|
} |
||||
|
|
||||
|
if (!ComputeVectors) |
||||
|
{ |
||||
|
throw new InvalidOperationException(Resources.SingularVectorsNotComputed); |
||||
|
} |
||||
|
|
||||
|
Matrix result = MatrixU.CreateMatrix(MatrixVT.ColumnCount, input.ColumnCount); |
||||
|
Solve(input, result); |
||||
|
return result; |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Solves a system of linear equations, <b>AX = B</b>, with A SVD factorized.
|
||||
|
/// </summary>
|
||||
|
/// <param name="input">The right hand side <see cref="Matrix"/>, <b>B</b>.</param>
|
||||
|
/// <param name="result">The left hand side <see cref="Matrix"/>, <b>X</b>.</param>
|
||||
|
public abstract void Solve(Matrix input, Matrix result); |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Solves a system of linear equations, <b>Ax = b</b>, with A SVD factorized.
|
||||
|
/// </summary>
|
||||
|
/// <param name="input">The right hand side vector, <b>b</b>.</param>
|
||||
|
/// <returns>The left hand side <see cref="Vector"/>, <b>x</b>.</returns>
|
||||
|
public virtual Vector Solve(Vector input) |
||||
|
{ |
||||
|
// Check for proper arguments.
|
||||
|
if (input == null) |
||||
|
{ |
||||
|
throw new ArgumentNullException("input"); |
||||
|
} |
||||
|
|
||||
|
if (!ComputeVectors) |
||||
|
{ |
||||
|
throw new InvalidOperationException(Resources.SingularVectorsNotComputed); |
||||
|
} |
||||
|
|
||||
|
var x = MatrixU.CreateVector(MatrixVT.ColumnCount); |
||||
|
Solve(input, x); |
||||
|
return x; |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Solves a system of linear equations, <b>Ax = b</b>, with A SVD factorized.
|
||||
|
/// </summary>
|
||||
|
/// <param name="input">The right hand side vector, <b>b</b>.</param>
|
||||
|
/// <param name="result">The left hand side <see cref="Matrix"/>, <b>x</b>.</param>
|
||||
|
public abstract void Solve(Vector input, Vector result); |
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,370 @@ |
|||||
|
// <copyright file="SvdTests.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.Factorization |
||||
|
{ |
||||
|
using System; |
||||
|
using MbUnit.Framework; |
||||
|
using LinearAlgebra.Double; |
||||
|
using LinearAlgebra.Double.Factorization; |
||||
|
|
||||
|
public class SvdTests |
||||
|
{ |
||||
|
|
||||
|
[Test] |
||||
|
[ExpectedArgumentNullException] |
||||
|
public void ConstructorNull() |
||||
|
{ |
||||
|
new DenseSvd(null, true); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(10)] |
||||
|
[Row(100)] |
||||
|
public void CanFactorizeIdentity(int order) |
||||
|
{ |
||||
|
var I = DenseMatrix.Identity(order); |
||||
|
var factorSvd = I.Svd(true); |
||||
|
|
||||
|
Assert.AreEqual(I.RowCount, factorSvd.U().RowCount); |
||||
|
Assert.AreEqual(I.RowCount, factorSvd.U().ColumnCount); |
||||
|
|
||||
|
Assert.AreEqual(I.ColumnCount, factorSvd.VT().RowCount); |
||||
|
Assert.AreEqual(I.ColumnCount, factorSvd.VT().ColumnCount); |
||||
|
|
||||
|
Assert.AreEqual(I.RowCount, factorSvd.W().RowCount); |
||||
|
Assert.AreEqual(I.ColumnCount, factorSvd.W().ColumnCount); |
||||
|
|
||||
|
for (var i = 0; i < factorSvd.W().RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < factorSvd.W().ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(i == j ? 1.0 : 0.0, factorSvd.W()[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1,1)] |
||||
|
[Row(2,2)] |
||||
|
[Row(5,5)] |
||||
|
[Row(10,6)] |
||||
|
[Row(48,52)] |
||||
|
[Row(100,93)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanFactorizeRandomMatrix(int row, int column) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomMatrix(row, column); |
||||
|
var factorSvd = matrixA.Svd(true); |
||||
|
|
||||
|
// Make sure the U has the right dimensions.
|
||||
|
Assert.AreEqual(row, factorSvd.U().RowCount); |
||||
|
Assert.AreEqual(row, factorSvd.U().ColumnCount); |
||||
|
|
||||
|
// Make sure the VT has the right dimensions.
|
||||
|
Assert.AreEqual(column, factorSvd.VT().RowCount); |
||||
|
Assert.AreEqual(column, factorSvd.VT().ColumnCount); |
||||
|
|
||||
|
// Make sure the W has the right dimensions.
|
||||
|
Assert.AreEqual(row, factorSvd.W().RowCount); |
||||
|
Assert.AreEqual(column, factorSvd.W().ColumnCount); |
||||
|
|
||||
|
// Make sure the U*W*VT is the original matrix.
|
||||
|
var matrix = factorSvd.U() * factorSvd.W() * factorSvd.VT(); |
||||
|
for (var i = 0; i < matrix.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrix.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixA[i, j], matrix[i, j], 1.0e-11); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(10, 8)] |
||||
|
[Row(48, 52)] |
||||
|
[Row(100, 93)] |
||||
|
[MultipleAsserts] |
||||
|
public void CheckRankOfNonSquare(int row, int column) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomMatrix(row, column); |
||||
|
var factorSvd = matrixA.Svd(true); |
||||
|
|
||||
|
var mn = Math.Min(row, column); |
||||
|
Assert.AreEqual(factorSvd.Rank, mn); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1)] |
||||
|
[Row(2)] |
||||
|
[Row(5)] |
||||
|
[Row(9)] |
||||
|
[Row(50)] |
||||
|
[Row(90)] |
||||
|
[MultipleAsserts] |
||||
|
public void CheckRankSquare(int order) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); |
||||
|
var factorSvd = matrixA.Svd(true); |
||||
|
|
||||
|
if (factorSvd.Determinant != 0) |
||||
|
{ |
||||
|
Assert.AreEqual(factorSvd.Rank, order); |
||||
|
} |
||||
|
else |
||||
|
{ |
||||
|
Assert.AreEqual(factorSvd.Rank, order - 1); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(10)] |
||||
|
[Row(50)] |
||||
|
[Row(100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CheckRankOfSquareSingular(int order) |
||||
|
{ |
||||
|
var matrixA = new DenseMatrix(order, order); |
||||
|
matrixA[0, 0] = 1; |
||||
|
matrixA[order - 1, order - 1] = 1; |
||||
|
for (var i = 1; i < order - 1; i++) |
||||
|
{ |
||||
|
matrixA[i, i - 1] = 1; |
||||
|
matrixA[i, i + 1] = 1; |
||||
|
matrixA[i - 1, i] = 1; |
||||
|
matrixA[i + 1, i] = 1; |
||||
|
} |
||||
|
var factorSvd = matrixA.Svd(true); |
||||
|
|
||||
|
Assert.AreEqual(factorSvd.Determinant, 0); |
||||
|
Assert.AreEqual(factorSvd.Rank, order - 1); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[ExpectedException(typeof(InvalidOperationException))] |
||||
|
public void CannotSolveMatrixIfVectorsNotComputed() |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomMatrix(10, 10); |
||||
|
var factorSvd = matrixA.Svd(false); |
||||
|
|
||||
|
var matrixB = MatrixLoader.GenerateRandomMatrix(10, 10); |
||||
|
factorSvd.Solve(matrixB); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[ExpectedException(typeof(InvalidOperationException))] |
||||
|
public void CannotSolveVectorIfVectorsNotComputed() |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomMatrix(10, 10); |
||||
|
var factorSvd = matrixA.Svd(false); |
||||
|
|
||||
|
var vectorb = MatrixLoader.GenerateRandomVector(10); |
||||
|
factorSvd.Solve(vectorb); |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1, 1)] |
||||
|
[Row(2, 2)] |
||||
|
[Row(5, 5)] |
||||
|
[Row(9, 10)] |
||||
|
[Row(50, 50)] |
||||
|
[Row(90, 100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomVector(int row, int column) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomMatrix(row, column); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorSvd = matrixA.Svd(true); |
||||
|
|
||||
|
var vectorb = MatrixLoader.GenerateRandomVector(row); |
||||
|
var resultx = factorSvd.Solve(vectorb); |
||||
|
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, resultx.Count); |
||||
|
|
||||
|
var bReconstruct = matrixA * resultx; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1.0e-11); |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1, 1)] |
||||
|
[Row(4, 4)] |
||||
|
[Row(7, 8)] |
||||
|
[Row(10, 10)] |
||||
|
[Row(45, 50)] |
||||
|
[Row(80, 100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomMatrix(int row, int count) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomMatrix(row, count); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorSvd = matrixA.Svd(true); |
||||
|
|
||||
|
var matrixB = MatrixLoader.GenerateRandomMatrix(row, count); |
||||
|
var matrixX = factorSvd.Solve(matrixB); |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); |
||||
|
// The solution X has the same number of columns as B
|
||||
|
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); |
||||
|
|
||||
|
var matrixBReconstruct = matrixA * matrixX; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1, 1)] |
||||
|
[Row(2, 2)] |
||||
|
[Row(5, 5)] |
||||
|
[Row(9, 10)] |
||||
|
[Row(50, 50)] |
||||
|
[Row(90, 100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomVectorWhenResultVectorGiven(int row, int column) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomMatrix(row, column); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorSvd = matrixA.Svd(true); |
||||
|
var vectorb = MatrixLoader.GenerateRandomVector(row); |
||||
|
var vectorbCopy = vectorb.Clone(); |
||||
|
var resultx = new DenseVector(column); |
||||
|
factorSvd.Solve(vectorb,resultx); |
||||
|
|
||||
|
var bReconstruct = matrixA * resultx; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1.0e-11); |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure b didn't change.
|
||||
|
for (var i = 0; i < vectorb.Count; i++) |
||||
|
{ |
||||
|
Assert.AreEqual(vectorbCopy[i], vectorb[i]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
[Test] |
||||
|
[Row(1, 1)] |
||||
|
[Row(4, 4)] |
||||
|
[Row(7, 8)] |
||||
|
[Row(10, 10)] |
||||
|
[Row(45, 50)] |
||||
|
[Row(80, 100)] |
||||
|
[MultipleAsserts] |
||||
|
public void CanSolveForRandomMatrixWhenResultMatrixGiven(int row, int column) |
||||
|
{ |
||||
|
var matrixA = MatrixLoader.GenerateRandomMatrix(row, column); |
||||
|
var matrixACopy = matrixA.Clone(); |
||||
|
var factorSvd = matrixA.Svd(true); |
||||
|
|
||||
|
var matrixB = MatrixLoader.GenerateRandomMatrix(row, column); |
||||
|
var matrixBCopy = matrixB.Clone(); |
||||
|
|
||||
|
var matrixX = new DenseMatrix(column, column); |
||||
|
factorSvd.Solve(matrixB,matrixX); |
||||
|
|
||||
|
// The solution X row dimension is equal to the column dimension of A
|
||||
|
Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); |
||||
|
// The solution X has the same number of columns as B
|
||||
|
Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); |
||||
|
|
||||
|
var matrixBReconstruct = matrixA * matrixX; |
||||
|
|
||||
|
// Check the reconstruction.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure A didn't change.
|
||||
|
for (var i = 0; i < matrixA.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixA.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// Make sure B didn't change.
|
||||
|
for (var i = 0; i < matrixB.RowCount; i++) |
||||
|
{ |
||||
|
for (var j = 0; j < matrixB.ColumnCount; j++) |
||||
|
{ |
||||
|
Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
|
} |
||||
Loading…
Reference in new issue