forked from tsai/mathnet-numerics
Marcus Cuda
16 years ago
13 changed files with 1623 additions and 648 deletions
File diff suppressed because it is too large
@ -0,0 +1,183 @@ |
|||
// <copyright file="DenseQR.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; |
|||
using Threading; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of the QR decomposition.</para>
|
|||
/// <para>Any real square matrix A may be decomposed as A = QR where Q is an orthogonal matrix
|
|||
/// (its columns are orthogonal unit vectors meaning QTQ = I) and R is an upper triangular matrix
|
|||
/// (also called right triangular matrix).</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the QR decomposition is done at construction time by Householder transformation.
|
|||
/// </remarks>
|
|||
public class DenseQR : QR |
|||
{ |
|||
/// <summary>
|
|||
/// Initializes a new instance of the <see cref="DenseQR"/> class. This object will compute the
|
|||
/// QR factorization 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 <c>null</c>.</exception>
|
|||
public DenseQR(DenseMatrix matrix) |
|||
{ |
|||
if (matrix == null) |
|||
{ |
|||
throw new ArgumentNullException("matrix"); |
|||
} |
|||
|
|||
if (matrix.RowCount < matrix.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixDimensions); |
|||
} |
|||
|
|||
MatrixR = matrix.Clone(); |
|||
MatrixQ = new DenseMatrix(matrix.RowCount); |
|||
Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, ((DenseMatrix)MatrixQ).Data); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Solves a system of linear equations, <b>AX = B</b>, with A QR 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"); |
|||
} |
|||
|
|||
// 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 (MatrixR.RowCount != input.RowCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); |
|||
} |
|||
|
|||
// The solution X row dimension is equal to the column dimension of A
|
|||
if (MatrixR.ColumnCount != result.RowCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); |
|||
} |
|||
|
|||
var dinput = input as DenseMatrix; |
|||
if (dinput == null) |
|||
{ |
|||
throw new NotImplementedException("Can only do QR factorization for dense matrices at the moment."); |
|||
} |
|||
|
|||
var dresult = result as DenseMatrix; |
|||
if (dresult == null) |
|||
{ |
|||
throw new NotImplementedException("Can only do QR factorization for dense matrices at the moment."); |
|||
} |
|||
|
|||
var solution = new double[dinput.Data.Length]; |
|||
Control.LinearAlgebraProvider.QRSolveFactored(input.ColumnCount, ((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, dinput.Data, solution); |
|||
CommonParallel.For( |
|||
0, |
|||
dresult.RowCount, |
|||
row => |
|||
{ |
|||
for (var col = 0; col < dresult.ColumnCount; col++) |
|||
{ |
|||
dresult[row, col] = solution[row + (col * dinput.RowCount)]; |
|||
} |
|||
}); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Solves a system of linear equations, <b>Ax = b</b>, with A QR 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"); |
|||
} |
|||
|
|||
// Ax=b where A is an m x n matrix
|
|||
// Check that b is a column vector with m entries
|
|||
if (MatrixR.RowCount != input.Count) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentVectorsSameLength); |
|||
} |
|||
|
|||
// Check that x is a column vector with n entries
|
|||
if (MatrixR.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."); |
|||
} |
|||
|
|||
var solution = new double[dinput.Data.Length]; |
|||
Control.LinearAlgebraProvider.QRSolveFactored(1, ((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, dinput.Data, solution); |
|||
CommonParallel.For( |
|||
0, |
|||
dresult.Count, |
|||
index => { dresult[index] = solution[index]; }); |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,195 @@ |
|||
// <copyright file="QR.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 QR decomposition.</para>
|
|||
/// <para>Any real square matrix A (m x n) may be decomposed as A = QR where Q is an orthogonal matrix (m x m)
|
|||
/// (its columns are orthogonal unit vectors meaning QTQ = I) and R (m x n) is an upper triangular matrix
|
|||
/// (also called right triangular matrix).</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the QR decomposition is done at construction time by Householder transformation.
|
|||
/// </remarks>
|
|||
public abstract class QR : ISolver |
|||
{ |
|||
/// <summary>
|
|||
/// Gets or sets orthogonal Q matrix
|
|||
/// </summary>
|
|||
protected virtual Matrix MatrixQ |
|||
{ |
|||
get; |
|||
set; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets or sets upper triangular factor R
|
|||
/// </summary>
|
|||
protected virtual Matrix MatrixR |
|||
{ |
|||
get; |
|||
set; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Internal method which routes the call to perform the QR factorization to the appropriate class.
|
|||
/// </summary>
|
|||
/// <param name="matrix">The matrix to factor.</param>
|
|||
/// <returns>A QR factorization object.</returns>
|
|||
internal static QR Create(Matrix matrix) |
|||
{ |
|||
var dense = matrix as DenseMatrix; |
|||
if (dense != null) |
|||
{ |
|||
return new DenseQR(dense); |
|||
} |
|||
|
|||
throw new NotImplementedException(); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets orthogonal Q matrix
|
|||
/// </summary>
|
|||
public virtual Matrix Q |
|||
{ |
|||
get |
|||
{ |
|||
return MatrixQ; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the upper triangular factor R.
|
|||
/// </summary>
|
|||
public virtual Matrix R |
|||
{ |
|||
get |
|||
{ |
|||
return MatrixR.UpperTriangle(); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the determinant of the matrix for which the QR matrix was computed.
|
|||
/// </summary>
|
|||
public virtual double Determinant |
|||
{ |
|||
get |
|||
{ |
|||
if (MatrixR.RowCount != MatrixR.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
|||
} |
|||
|
|||
var det = 1.0; |
|||
for (var i = 0; i < MatrixR.ColumnCount; i++) |
|||
{ |
|||
det *= MatrixR.At(i, i); |
|||
if (Math.Abs(MatrixR.At(i, i)).AlmostEqualInDecimalPlaces(0.0, 15)) |
|||
{ |
|||
return 0; |
|||
} |
|||
} |
|||
|
|||
return Math.Abs(det); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets a value indicating whether the matrix is full rank or not.
|
|||
/// </summary>
|
|||
/// <value><c>true</c> if the matrix is full rank; otherwise <c>false</c>.</value>
|
|||
public virtual bool IsFullRank |
|||
{ |
|||
get |
|||
{ |
|||
for (var i = 0; i < MatrixR.ColumnCount; i++) |
|||
{ |
|||
if (Math.Abs(MatrixR.At(i, i)).AlmostEqualInDecimalPlaces(0.0, 15)) |
|||
{ |
|||
return false; |
|||
} |
|||
} |
|||
|
|||
return true; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Solves a system of linear equations, <b>AX = B</b>, with A QR 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"); |
|||
} |
|||
|
|||
var matrixX = input.CreateMatrix(MatrixR.ColumnCount, input.ColumnCount); |
|||
Solve(input, matrixX); |
|||
return matrixX; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Solves a system of linear equations, <b>AX = B</b>, with A QR 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 QR 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"); |
|||
} |
|||
|
|||
var x = input.CreateVector(MatrixR.ColumnCount); |
|||
Solve(input, x); |
|||
return x; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Solves a system of linear equations, <b>Ax = b</b>, with A QR 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,64 @@ |
|||
// <copyright file="ISolver.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 |
|||
{ |
|||
/// <summary>
|
|||
/// Classes that solves a system of linear equations, <c>AX = B</c>.
|
|||
/// </summary>
|
|||
public interface ISolver |
|||
{ |
|||
/// <summary>
|
|||
/// Solves a system of linear equations, <c>AX = B</c>.
|
|||
/// </summary>
|
|||
/// <param name="input">The right hand side <see cref="Matrix"/>, <c>B</c>.</param>
|
|||
/// <returns>The left hand side <see cref="Matrix"/>, <c>X</c>.</returns>
|
|||
Matrix Solve(Matrix input); |
|||
|
|||
/// <summary>
|
|||
/// Solves a system of linear equations, <c>AX = B</c>.
|
|||
/// </summary>
|
|||
/// <param name="input">The right hand side <see cref="Matrix"/>, <c>B</c>.</param>
|
|||
/// <param name="result">The left hand side <see cref="Matrix"/>, <c>X</c>.</param>
|
|||
void Solve(Matrix input, Matrix result); |
|||
|
|||
/// <summary>
|
|||
/// Solves a system of linear equations, <c>Ax = b</c>
|
|||
/// </summary>
|
|||
/// <param name="input">The right hand side vector, <c>b</c>.</param>
|
|||
/// <returns>The left hand side <see cref="Vector"/>, <c>x</c>.</returns>
|
|||
Vector Solve(Vector input); |
|||
|
|||
/// <summary>
|
|||
/// Solves a system of linear equations, <c>Ax = b</c>.
|
|||
/// </summary>
|
|||
/// <param name="input">The right hand side vector, <c>b</c>.</param>
|
|||
/// <param name="result">The left hand side <see cref="Matrix"/>, <c>x</c>.</param>
|
|||
void Solve(Vector input, Vector result); |
|||
} |
|||
} |
|||
@ -0,0 +1,317 @@ |
|||
// <copyright file="QRTests.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 MbUnit.Framework; |
|||
using LinearAlgebra.Double; |
|||
using LinearAlgebra.Double.Factorization; |
|||
|
|||
public class QRTests |
|||
{ |
|||
|
|||
[Test] |
|||
[ExpectedArgumentNullException] |
|||
public void ConstructorNull() |
|||
{ |
|||
new DenseQR(null); |
|||
} |
|||
|
|||
[Test] |
|||
[ExpectedArgumentException] |
|||
public void WideMatrixThrowsInvalidMatrixOperationException() |
|||
{ |
|||
new DenseQR(new DenseMatrix(3, 4)); |
|||
} |
|||
|
|||
[Test] |
|||
[Row(1)] |
|||
[Row(10)] |
|||
[Row(100)] |
|||
public void CanFactorizeIdentity(int order) |
|||
{ |
|||
var I = DenseMatrix.Identity(order); |
|||
var factorQR = I.QR(); |
|||
|
|||
Assert.AreEqual(I.RowCount, factorQR.R.RowCount); |
|||
Assert.AreEqual(I.ColumnCount, factorQR.R.ColumnCount); |
|||
|
|||
for (var i = 0; i < factorQR.R.RowCount; i++) |
|||
{ |
|||
for (var j = 0; j < factorQR.R.ColumnCount; j++) |
|||
{ |
|||
if (i == j) |
|||
{ |
|||
Assert.AreEqual(-1.0, factorQR.R[i, j]); |
|||
} |
|||
else |
|||
{ |
|||
Assert.AreEqual(0.0, factorQR.R[i, j]); |
|||
} |
|||
} |
|||
} |
|||
} |
|||
|
|||
|
|||
[Test] |
|||
[Row(1)] |
|||
[Row(10)] |
|||
[Row(100)] |
|||
public void IdentityDeterminantIsOne(int order) |
|||
{ |
|||
var I = DenseMatrix.Identity(order); |
|||
var factorQR = I.QR(); |
|||
Assert.AreEqual(1.0, factorQR.Determinant); |
|||
} |
|||
|
|||
[Test] |
|||
[Row(1,1)] |
|||
[Row(2,2)] |
|||
[Row(5,5)] |
|||
[Row(10,6)] |
|||
[Row(50,48)] |
|||
[Row(100,98)] |
|||
[MultipleAsserts] |
|||
public void CanFactorizeRandomMatrix(int row, int column) |
|||
{ |
|||
var matrixA = MatrixLoader.GenerateRandomMatrix(row, column); |
|||
var factorQR = matrixA.QR(); |
|||
|
|||
// Make sure the R has the right dimensions.
|
|||
Assert.AreEqual(row, factorQR.R.RowCount); |
|||
Assert.AreEqual(column, factorQR.R.ColumnCount); |
|||
|
|||
// Make sure the Q has the right dimensions.
|
|||
Assert.AreEqual(row, factorQR.Q.RowCount); |
|||
Assert.AreEqual(row, factorQR.Q.ColumnCount); |
|||
|
|||
// Make sure the R factor is upper triangular.
|
|||
for (var i = 0; i < factorQR.R.RowCount; i++) |
|||
{ |
|||
for (var j = 0; j < factorQR.R.ColumnCount; j++) |
|||
{ |
|||
if (i > j) |
|||
{ |
|||
Assert.AreEqual(0.0, factorQR.R[i, j]); |
|||
} |
|||
} |
|||
} |
|||
|
|||
// Make sure the Q*R is the original matrix.
|
|||
var matrixQfromR = factorQR.Q * factorQR.R; |
|||
for (int i = 0; i < matrixQfromR.RowCount; i++) |
|||
{ |
|||
for (int j = 0; j < matrixQfromR.ColumnCount; j++) |
|||
{ |
|||
Assert.AreApproximatelyEqual(matrixA[i, j], matrixQfromR[i, j], 1.0e-11); |
|||
} |
|||
} |
|||
} |
|||
|
|||
[Test] |
|||
[Row(1)] |
|||
[Row(2)] |
|||
[Row(5)] |
|||
[Row(10)] |
|||
[Row(50)] |
|||
[Row(100)] |
|||
[MultipleAsserts] |
|||
public void CanSolveForRandomVector(int order) |
|||
{ |
|||
var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); |
|||
var matrixACopy = matrixA.Clone(); |
|||
var factorQR = matrixA.QR(); |
|||
|
|||
var vectorb = MatrixLoader.GenerateRandomVector(order); |
|||
var resultx = factorQR.Solve(vectorb); |
|||
|
|||
Assert.AreEqual(matrixA.ColumnCount, resultx.Count); |
|||
|
|||
var bReconstruct = matrixA * resultx; |
|||
|
|||
// Check the reconstruction.
|
|||
for (var i = 0; i < order; 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)] |
|||
[Row(4)] |
|||
[Row(8)] |
|||
[Row(10)] |
|||
[Row(50)] |
|||
[Row(100)] |
|||
[MultipleAsserts] |
|||
public void CanSolveForRandomMatrix(int order) |
|||
{ |
|||
var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); |
|||
var matrixACopy = matrixA.Clone(); |
|||
var factorQR = matrixA.QR(); |
|||
|
|||
var matrixB = MatrixLoader.GenerateRandomMatrix(order, order); |
|||
var matrixX = factorQR.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)] |
|||
[Row(2)] |
|||
[Row(5)] |
|||
[Row(10)] |
|||
[Row(50)] |
|||
[Row(100)] |
|||
[MultipleAsserts] |
|||
public void CanSolveForRandomVectorWhenResultVectorGiven(int order) |
|||
{ |
|||
var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); |
|||
var matrixACopy = matrixA.Clone(); |
|||
var factorQR = matrixA.QR(); |
|||
var vectorb = MatrixLoader.GenerateRandomVector(order); |
|||
var vectorbCopy = vectorb.Clone(); |
|||
var resultx = new DenseVector(order); |
|||
factorQR.Solve(vectorb,resultx); |
|||
|
|||
Assert.AreEqual(vectorb.Count, 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]); |
|||
} |
|||
} |
|||
|
|||
// Make sure b didn't change.
|
|||
for (var i = 0; i < vectorb.Count; i++) |
|||
{ |
|||
Assert.AreEqual(vectorbCopy[i], vectorb[i]); |
|||
} |
|||
} |
|||
|
|||
[Test] |
|||
[Row(1)] |
|||
[Row(4)] |
|||
[Row(8)] |
|||
[Row(10)] |
|||
[Row(50)] |
|||
[Row(100)] |
|||
[MultipleAsserts] |
|||
public void CanSolveForRandomMatrixWhenResultMatrixGiven(int order) |
|||
{ |
|||
var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); |
|||
var matrixACopy = matrixA.Clone(); |
|||
var factorQR = matrixA.QR(); |
|||
|
|||
var matrixB = MatrixLoader.GenerateRandomMatrix(order, order); |
|||
var matrixBCopy = matrixB.Clone(); |
|||
|
|||
var matrixX = new DenseMatrix(order, order); |
|||
factorQR.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