diff --git a/src/MathNet.Numerics.5.0.ReSharper b/src/MathNet.Numerics.5.0.ReSharper index 20bd5dbf..995c83e8 100644 --- a/src/MathNet.Numerics.5.0.ReSharper +++ b/src/MathNet.Numerics.5.0.ReSharper @@ -863,6 +863,7 @@ Cholesky + diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs index 54c8a990..d7e3d74d 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs @@ -1361,7 +1361,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// Perform calculation of Q or R /// - /// Work arrat + /// Work array /// Index of colunn in work array /// Q or R matrices /// The number of rows @@ -3053,16 +3053,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } /// - /// Сonstruct givens plane rotation + /// Given the Cartesian coordinates (da, db) of a point p, these fucntion return the parameters da, db, c, and s + /// associated with the Givens rotation that zeros the y-coordinate of the point. /// - /// - /// - /// - /// + /// Provides the x-coordinate of the point p. On exit contains the parameter r associated with the Givens rotation + /// Provides the y-coordinate of the point p. On exit contains the parameter z associated with the Givens rotation + /// Contains the parameter c associated with the Givens rotation + /// Contains the parameter s associated with the Givens rotation + /// This is equivalent to the DROTG LAPACK routine. private static void Drotg(ref double da, ref double db, ref double c, ref double s) { - // Сonstruct givens plane rotation. - // jack dongarra, linpack, 3/11/78. double r, z; var roe = db; @@ -3107,7 +3107,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra da = r; db = z; - return; } /// diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs b/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs index 800e528f..140ffb07 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs @@ -56,16 +56,27 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization return new DenseCholesky(dense); } - throw new NotImplementedException(); + return new UserCholesky(matrix); } /// - /// Gets or sets the lower triangular form of the Cholesky matrix. + /// Gets or sets the lower triangular form of the Cholesky matrix /// - public virtual Matrix Factor + protected Matrix CholeskyFactor { get; - protected set; + set; + } + + /// + /// Gets the lower triangular form of the Cholesky matrix. + /// + public virtual Matrix Factor + { + get + { + return CholeskyFactor.Clone(); + } } /// @@ -76,9 +87,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization get { var det = 1.0; - for (var j = 0; j < Factor.RowCount; j++) + for (var j = 0; j < CholeskyFactor.RowCount; j++) { - det *= Factor[j, j] * Factor[j, j]; + det *= CholeskyFactor[j, j] * CholeskyFactor[j, j]; } return det; @@ -93,9 +104,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization get { var det = 0.0; - for (var j = 0; j < Factor.RowCount; j++) + for (var j = 0; j < CholeskyFactor.RowCount; j++) { - det += 2.0 * Math.Log(Factor[j, j]); + det += 2.0 * Math.Log(CholeskyFactor[j, j]); } return det; diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs index 0b1c0881..8b7785ad 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs @@ -67,7 +67,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting). var factor = (DenseMatrix)matrix.Clone(); Control.LinearAlgebraProvider.CholeskyFactor(factor.Data, factor.RowCount); - Factor = factor; + CholeskyFactor = factor; } /// @@ -99,7 +99,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); } - if (input.RowCount != Factor.RowCount) + if (input.RowCount != CholeskyFactor.RowCount) { throw new ArgumentException(Resources.ArgumentMatrixDimensions); } @@ -120,7 +120,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization Buffer.BlockCopy(dinput.Data, 0, dresult.Data, 0, dinput.Data.Length * Constants.SizeOfDouble); // Cholesky solve by overwriting result. - var dfactor = (DenseMatrix)Factor; + var dfactor = (DenseMatrix)CholeskyFactor; Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.RowCount, dresult.ColumnCount); } @@ -148,7 +148,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - if (input.Count != Factor.RowCount) + if (input.Count != CholeskyFactor.RowCount) { throw new ArgumentException(Resources.ArgumentMatrixDimensions); } @@ -169,7 +169,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization Buffer.BlockCopy(dinput.Data, 0, dresult.Data, 0, dinput.Data.Length * Constants.SizeOfDouble); // Cholesky solve by overwriting result. - var dfactor = (DenseMatrix)Factor; + var dfactor = (DenseMatrix)CholeskyFactor; Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.Count, 1); } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs b/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs index 2b0c35e5..03bdce8c 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs @@ -71,7 +71,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization return new DenseLU(dense); } - throw new NotImplementedException(); + return new UserLU(matrix); } /// diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/QR.cs b/src/Numerics/LinearAlgebra/Double/Factorization/QR.cs index e52326fe..c13ac636 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/QR.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/QR.cs @@ -71,7 +71,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization return new DenseQR(dense); } - throw new NotImplementedException(); + return new UserQR(matrix); } /// diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/Svd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/Svd.cs index bb5f9962..7440a5a6 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/Svd.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/Svd.cs @@ -123,7 +123,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization return new DenseSvd(dense, computeVectors); } - throw new NotImplementedException(); + return new UserSvd(matrix, computeVectors); } /// diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/UserCholesky.cs b/src/Numerics/LinearAlgebra/Double/Factorization/UserCholesky.cs new file mode 100644 index 00000000..88962b83 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Double/Factorization/UserCholesky.cs @@ -0,0 +1,222 @@ +// +// 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. +// + +namespace MathNet.Numerics.LinearAlgebra.Double.Factorization +{ + using System; + using Properties; + + /// + /// A class which encapsulates the functionality of a Cholesky factorization for user matrices. + /// For a symmetric, positive definite matrix A, the Cholesky factorization + /// is an lower triangular matrix L so that A = L*L'. + /// + /// + /// The computation of the Cholesky factorization is done at construction time. If the matrix is not symmetric + /// or positive definite, the constructor will throw an exception. + /// + public class UserCholesky : Cholesky + { + /// + /// Initializes a new instance of the class. This object will compute the + /// Cholesky factorization when the constructor is called and cache it's factorization. + /// + /// The matrix to factor. + /// If is null. + /// If is not a square matrix. + /// If is not positive definite. + public UserCholesky(Matrix 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). + CholeskyFactor = matrix.Clone(); + for (var j = 0; j < CholeskyFactor.RowCount; j++) + { + var d = 0.0; + for (var k = 0; k < j; k++) + { + var s = 0.0; + for (var i = 0; i < k; i++) + { + s += CholeskyFactor.At(k, i) * CholeskyFactor.At(j, i); + } + + s = (matrix.At(j, k) - s) / CholeskyFactor.At(k, k); + CholeskyFactor.At(j, k, s); + d += s * s; + } + + d = matrix.At(j, j) - d; + if (d <= 0.0) + { + throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite); + } + + CholeskyFactor.At(j, j, Math.Sqrt(d)); + for (var k = j + 1; k < CholeskyFactor.RowCount; k++) + { + CholeskyFactor.At(j, k, 0.0); + } + } + } + + /// + /// Solves a system of linear equations, AX = B, with A Cholesky factorized. + /// + /// The right hand side , B. + /// The left hand side , X. + public override void Solve(Matrix input, Matrix result) + { + if (input == null) + { + throw new ArgumentNullException("input"); + } + + if (result == null) + { + throw new ArgumentNullException("result"); + } + + // Check for proper dimensions. + if (result.RowCount != input.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); + } + + if (result.ColumnCount != input.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); + } + + if (input.RowCount != CholeskyFactor.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + input.CopyTo(result); + var order = CholeskyFactor.RowCount; + + for (var c = 0; c < result.ColumnCount; c++) + { + // Solve L*Y = B; + double sum; + for (var i = 0; i < order; i++) + { + sum = result.At(i, c); + for (var k = i - 1; k >= 0; k--) + { + sum -= CholeskyFactor.At(i, k) * result.At(k, c); + } + + result.At(i, c, sum / CholeskyFactor.At(i, i)); + } + + // Solve L'*X = Y; + for (var i = order - 1; i >= 0; i--) + { + sum = result.At(i, c); + for (var k = i + 1; k < order; k++) + { + sum -= CholeskyFactor.At(k, i) * result.At(k, c); + } + + result.At(i, c, sum / CholeskyFactor.At(i, i)); + } + } + } + + /// + /// Solves a system of linear equations, Ax = b, with A Cholesky factorized. + /// + /// The right hand side vector, b. + /// The left hand side , x. + public override void Solve(Vector input, Vector result) + { + // Check for proper arguments. + if (input == null) + { + throw new ArgumentNullException("input"); + } + + if (result == null) + { + throw new ArgumentNullException("result"); + } + + // Check for proper dimensions. + if (input.Count != result.Count) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + if (input.Count != CholeskyFactor.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + input.CopyTo(result); + var order = CholeskyFactor.RowCount; + + // Solve L*Y = B; + double sum; + for (var i = 0; i < order; i++) + { + sum = result[i]; + for (var k = i - 1; k >= 0; k--) + { + sum -= CholeskyFactor.At(i, k) * result[k]; + } + + result[i] = sum / CholeskyFactor.At(i, i); + } + + // Solve L'*X = Y; + for (var i = order - 1; i >= 0; i--) + { + sum = result[i]; + for (var k = i + 1; k < order; k++) + { + sum -= CholeskyFactor.At(k, i) * result[k]; + } + + result[i] = sum / CholeskyFactor.At(i, i); + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/UserLU.cs b/src/Numerics/LinearAlgebra/Double/Factorization/UserLU.cs new file mode 100644 index 00000000..de41dcd5 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Double/Factorization/UserLU.cs @@ -0,0 +1,300 @@ +// +// 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. +// + +namespace MathNet.Numerics.LinearAlgebra.Double.Factorization +{ + using System; + using Properties; + + /// + /// A class which encapsulates the functionality of an LU factorization. + /// For a matrix A, the LU factorization is a pair of lower triangular matrix L and + /// upper triangular matrix U so that A = L*U. + /// + /// + /// The computation of the LU factorization is done at construction time. + /// + public class UserLU : LU + { + /// + /// Initializes a new instance of the class. This object will compute the + /// LU factorization when the constructor is called and cache it's factorization. + /// + /// The matrix to factor. + /// If is null. + /// If is not a square matrix. + public UserLU(Matrix matrix) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (matrix.RowCount != matrix.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + // Create an array for the pivot indices. + var order = matrix.RowCount; + Factors = matrix.Clone(); + Pivots = new int[order]; + + // Initialize the pivot matrix to the identity permutation. + for (var i = 0; i < order; i++) + { + Pivots[i] = i; + } + + var vectorLUcolj = new double[order]; + for (var j = 0; j < order; j++) + { + // Make a copy of the j-th column to localize references. + for (var i = 0; i < order; i++) + { + vectorLUcolj[i] = Factors.At(i, j); + } + + // Apply previous transformations. + for (var i = 0; i < order; i++) + { + var kmax = Math.Min(i, j); + var s = 0.0; + for (var k = 0; k < kmax; k++) + { + s += Factors.At(i, k) * vectorLUcolj[k]; + } + + vectorLUcolj[i] -= s; + Factors.At(i, j, vectorLUcolj[i]); + } + + // Find pivot and exchange if necessary. + var p = j; + for (var i = j + 1; i < order; i++) + { + if (Math.Abs(vectorLUcolj[i]) > Math.Abs(vectorLUcolj[p])) + { + p = i; + } + } + + if (p != j) + { + for (var k = 0; k < order; k++) + { + var temp = Factors.At(p, k); + Factors.At(p, k, Factors.At(j, k)); + Factors.At(j, k, temp); + } + + Pivots[j] = p; + } + + // Compute multipliers. + if (j < order & Factors.At(j, j) != 0.0) + { + for (var i = j + 1; i < order; i++) + { + Factors.At(i, j, (Factors.At(i, j) / Factors.At(j, j))); + } + } + } + } + + /// + /// Solves a system of linear equations, AX = B, with A LU factorized. + /// + /// The right hand side , B. + /// The left hand side , X. + 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"); + } + + // Check for proper dimensions. + if (result.RowCount != input.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); + } + + if (result.ColumnCount != input.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); + } + + if (input.RowCount != Factors.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + // Copy the contents of input to result. + input.CopyTo(result); + for (var i = 0; i < Pivots.Length; i++) + { + if (Pivots[i] == i) + { + continue; + } + + var p = Pivots[i]; + for (var j = 0; j < result.ColumnCount; j++) + { + var temp = result.At(p, j); + result.At(p, j, result.At(i, j)); + result.At(i, j, temp); + } + } + + var order = Factors.RowCount; + + // Solve L*Y = P*B + for (var k = 0; k < order; k++) + { + for (var i = k + 1; i < order; i++) + { + for (var j = 0; j < result.ColumnCount; j++) + { + var temp = result.At(k, j) * Factors.At(i, k); + result.At(i, j, result.At(i, j) - temp); + } + } + } + + // Solve U*X = Y; + for (var k = order - 1; k >= 0; k--) + { + for (var j = 0; j < result.ColumnCount; j++) + { + result.At(k, j, (result.At(k, j) / Factors.At(k, k))); + } + + for (var i = 0; i < k; i++) + { + for (var j = 0; j < result.ColumnCount; j++) + { + var temp = result.At(k, j) * Factors.At(i, k); + result.At(i, j, result.At(i, j) - temp); + } + } + } + } + + /// + /// Solves a system of linear equations, Ax = b, with A LU factorized. + /// + /// The right hand side vector, b. + /// The left hand side , x. + public override void Solve(Vector input, Vector result) + { + // Check for proper arguments. + if (input == null) + { + throw new ArgumentNullException("input"); + } + + if (result == null) + { + throw new ArgumentNullException("result"); + } + + // Check for proper dimensions. + if (input.Count != result.Count) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + if (input.Count != Factors.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + // Copy the contents of input to result. + input.CopyTo(result); + for (var i = 0; i < Pivots.Length; i++) + { + if (Pivots[i] == i) + { + continue; + } + + var p = Pivots[i]; + var temp = result[p]; + result[p] = result[i]; + result[i] = temp; + } + + var order = Factors.RowCount; + + // Solve L*Y = P*B + for (var k = 0; k < order; k++) + { + for (var i = k + 1; i < order; i++) + { + result[i] -= result[k] * Factors.At(i, k); + } + } + + // Solve U*X = Y; + for (var k = order - 1; k >= 0; k--) + { + result[k] /= Factors.At(k, k); + for (var i = 0; i < k; i++) + { + result[i] -= result[k] * Factors.At(i, k); + } + } + } + + /// + /// Returns the inverse of this matrix. The inverse is calculated using LU decomposition. + /// + /// The inverse of this matrix. + public override Matrix Inverse() + { + var order = Factors.RowCount; + var inverse = Factors.CreateMatrix(order, order); + for (var i = 0; i < order; i++) + { + inverse.At(i, i, 1.0); + } + + return Solve(inverse); + } + } +} diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs b/src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs new file mode 100644 index 00000000..2595fadb --- /dev/null +++ b/src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs @@ -0,0 +1,332 @@ +// +// 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. +// + +namespace MathNet.Numerics.LinearAlgebra.Double.Factorization +{ + using System; + using System.Linq; + using Properties; + + /// + /// A class which encapsulates the functionality of the QR decomposition. + /// 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). + /// + /// + /// The computation of the QR decomposition is done at construction time by Householder transformation. + /// + public class UserQR : QR + { + /// + /// Initializes a new instance of the class. This object will compute the + /// QR factorization when the constructor is called and cache it's factorization. + /// + /// The matrix to factor. + /// If is null. + public UserQR(Matrix matrix) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (matrix.RowCount < matrix.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + MatrixR = matrix.Clone(); + MatrixQ = matrix.CreateMatrix(matrix.RowCount, matrix.RowCount); + + for (var i = 0; i < matrix.RowCount; i++) + { + MatrixQ.At(i, i, 1.0); + } + + var minmn = Math.Min(matrix.RowCount, matrix.ColumnCount); + var u = new double[minmn][]; + for (var i = 0; i < minmn; i++) + { + u[i] = GenerateColumn(MatrixR, i, matrix.RowCount - 1, i); + ComputeQR(u[i], MatrixR, i, matrix.RowCount - 1, i + 1, matrix.ColumnCount - 1); + } + + for (var i = minmn - 1; i >= 0; i--) + { + ComputeQR(u[i], MatrixQ, i, matrix.RowCount - 1, i, matrix.RowCount - 1); + } + } + + /// + /// Generate column from initial matrix to work array + /// + /// Initial matrix + /// The firts row + /// The last row + /// Column index + /// Generated vector + private static double[] GenerateColumn(Matrix a, int rowStart, int rowEnd, int column) + { + var ru = rowEnd - rowStart + 1; + var u = new double[ru]; + + for (var i = rowStart; i <= rowEnd; i++) + { + u[i - rowStart] = a.At(i, rowStart); + a.At(i, rowStart, 0.0); + } + + var norm = u.Sum(t => t * t); + norm = Math.Sqrt(norm); + + if (rowStart == rowEnd || norm == 0) + { + a.At(rowStart, column, -u[0]); + u[0] = Math.Sqrt(2.0); + return u; + } + + var scale = 1.0 / norm; + if (u[0] < 0.0) + { + scale *= -1.0; + } + + a.At(rowStart, column, -1.0 / scale); + + for (var i = 0; i < ru; i++) + { + u[i] *= scale; + } + + u[0] += 1.0; + var s = Math.Sqrt(1.0 / u[0]); + + for (var i = 0; i < ru; i++) + { + u[i] *= s; + } + + return u; + } + + /// + /// Perform calculation of Q or R + /// + /// Work array + /// Q or R matrices + /// The first row + /// The last row + /// The first column + /// The last column + private static void ComputeQR(double[] u, Matrix a, int rowStart, int rowEnd, int columnStart, int columnEnd) + { + if (rowEnd < rowStart || columnEnd < columnStart) + { + return; + } + + var v = new double[columnEnd - columnStart + 1]; + for (var j = columnStart; j <= columnEnd; j++) + { + v[j - columnStart] = 0.0; + } + + for (var i = rowStart; i <= rowEnd; i++) + { + for (var j = columnStart; j <= columnEnd; j++) + { + v[j - columnStart] = v[j - columnStart] + (u[i - rowStart] * a.At(i, j)); + } + } + + for (var i = rowStart; i <= rowEnd; i++) + { + for (var j = columnStart; j <= columnEnd; j++) + { + a.At(i, j, a.At(i, j) - (u[i - rowStart] * v[j - columnStart])); + } + } + } + + /// + /// Solves a system of linear equations, AX = B, with A QR factorized. + /// + /// The right hand side , B. + /// The left hand side , X. + 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 inputCopy = input.Clone(); + + // Compute Y = transpose(Q)*B + var bn = inputCopy.ColumnCount; + var column = new double[MatrixR.RowCount]; + for (var j = 0; j < bn; j++) + { + for (var k = 0; k < MatrixR.RowCount; k++) + { + column[k] = inputCopy.At(k, j); + } + + for (var i = 0; i < MatrixR.RowCount; i++) + { + double s = 0; + for (var k = 0; k < MatrixR.RowCount; k++) + { + s += MatrixQ.At(k, i) * column[k]; + } + + inputCopy.At(i, j, s); + } + } + + // Solve R*X = Y; + for (var k = MatrixR.ColumnCount - 1; k >= 0; k--) + { + for (var j = 0; j < bn; j++) + { + inputCopy.At(k, j, inputCopy.At(k, j) / MatrixR.At(k, k)); + } + + for (var i = 0; i < k; i++) + { + for (var j = 0; j < bn; j++) + { + inputCopy.At(i, j, inputCopy.At(i, j) - (inputCopy.At(k, j) * MatrixR.At(i, k))); + } + } + } + + for (var i = 0; i < MatrixR.ColumnCount; i++) + { + for (var j = 0; j < inputCopy.ColumnCount; j++) + { + result.At(i, j, inputCopy.At(i, j)); + } + } + } + + /// + /// Solves a system of linear equations, Ax = b, with A QR factorized. + /// + /// The right hand side vector, b. + /// The left hand side , x. + 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 inputCopy = input.Clone(); + + // Compute Y = transpose(Q)*B + var column = new double[MatrixR.RowCount]; + for (var k = 0; k < MatrixR.RowCount; k++) + { + column[k] = inputCopy[k]; + } + + for (var i = 0; i < MatrixR.RowCount; i++) + { + double s = 0; + for (var k = 0; k < MatrixR.RowCount; k++) + { + s += MatrixQ.At(k, i) * column[k]; + } + + inputCopy[i] = s; + } + + // Solve R*X = Y; + for (var k = MatrixR.ColumnCount - 1; k >= 0; k--) + { + inputCopy[k] /= MatrixR.At(k, k); + for (var i = 0; i < k; i++) + { + inputCopy[i] -= inputCopy[k] * MatrixR.At(i, k); + } + } + + for (var i = 0; i < MatrixR.ColumnCount; i++) + { + result[i] = inputCopy[i]; + } + } + } +} diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/UserSvd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/UserSvd.cs new file mode 100644 index 00000000..71b5762d --- /dev/null +++ b/src/Numerics/LinearAlgebra/Double/Factorization/UserSvd.cs @@ -0,0 +1,923 @@ +// +// 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. +// +namespace MathNet.Numerics.LinearAlgebra.Double.Factorization +{ + using System; + using Properties; + + /// + /// A class which encapsulates the functionality of the singular value decomposition (SVD) for . + /// 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. + /// + /// + /// The computation of the singular value decomposition is done at construction time. + /// + public class UserSvd : Svd + { + /// + /// Initializes a new instance of the class. This object will compute the + /// the singular value decomposition when the constructor is called and cache it's decomposition. + /// + /// The matrix to factor. + /// Compute the singular U and VT vectors or not. + /// If is null. + /// If SVD algorithm failed to converge with matrix . + public UserSvd(Matrix matrix, bool computeVectors) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + ComputeVectors = computeVectors; + var nm = Math.Min(matrix.RowCount + 1, matrix.ColumnCount); + var matrixCopy = matrix.Clone(); + + VectorS = matrixCopy.CreateVector(nm); + MatrixU = matrixCopy.CreateMatrix(matrixCopy.RowCount, matrixCopy.RowCount); + MatrixVT = matrixCopy.CreateMatrix(matrixCopy.ColumnCount, matrixCopy.ColumnCount); + + const int Maxiter = 1000; + var e = new double[matrixCopy.ColumnCount]; + var work = new double[matrixCopy.RowCount]; + + int i, j; + int l, lp1; + var cs = 0.0; + var sn = 0.0; + double t; + + var ncu = matrixCopy.RowCount; + + // Reduce matrixCopy to bidiagonal form, storing the diagonal elements + // In s and the super-diagonal elements in e. + var nct = Math.Min(matrixCopy.RowCount - 1, matrixCopy.ColumnCount); + var nrt = Math.Max(0, Math.Min(matrixCopy.ColumnCount - 2, matrixCopy.RowCount)); + var lu = Math.Max(nct, nrt); + for (l = 0; l < lu; l++) + { + lp1 = l + 1; + if (l < nct) + { + // Compute the transformation for the l-th column and place the l-th diagonal in VectorS[l]. + var xnorm = Dnrm2Column(matrixCopy, matrixCopy.RowCount, l, l); + VectorS[l] = xnorm; + if (VectorS[l] != 0.0) + { + if (matrixCopy.At(l, l) != 0.0) + { + VectorS[l] = Dsign(VectorS[l], matrixCopy.At(l, l)); + } + + DscalColumn(matrixCopy, matrixCopy.RowCount, l, l, 1.0 / VectorS[l]); + matrixCopy.At(l, l, (1.0 + matrixCopy.At(l, l))); + } + + VectorS[l] = -VectorS[l]; + } + + for (j = lp1; j < matrixCopy.ColumnCount; j++) + { + if (l < nct) + { + if (VectorS[l] != 0.0) + { + // Apply the transformation. + t = -Ddot(matrixCopy, matrixCopy.RowCount, l, j, l) / matrixCopy.At(l, l); + for (var ii = l; ii < matrixCopy.RowCount; ii++) + { + matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (t * matrixCopy.At(ii, l))); + } + } + } + + // Place the l-th row of matrixCopy into e for the + // Subsequent calculation of the row transformation. + e[j] = matrixCopy.At(l, j); + } + + if (ComputeVectors && l < nct) + { + // Place the transformation in u for subsequent back multiplication. + for (i = l; i < matrixCopy.RowCount; i++) + { + MatrixU.At(i, l, matrixCopy.At(i, l)); + } + } + + if (l >= nrt) + { + continue; + } + + // Compute the l-th row transformation and place the l-th super-diagonal in e(l). + var enorm = Dnrm2Vector(e, lp1); + e[l] = enorm; + if (e[l] != 0.0) + { + if (e[lp1] != 0.0) + { + e[l] = Dsign(e[l], e[lp1]); + } + + DscalVector(e, lp1, 1.0 / e[l]); + e[lp1] = 1.0 + e[lp1]; + } + + e[l] = -e[l]; + if (lp1 < matrixCopy.RowCount && e[l] != 0.0) + { + // Apply the transformation. + for (i = lp1; i < matrixCopy.RowCount; i++) + { + work[i] = 0.0; + } + + for (j = lp1; j < matrixCopy.ColumnCount; j++) + { + for (var ii = lp1; ii < matrixCopy.RowCount; ii++) + { + work[ii] += e[j] * matrixCopy.At(ii, j); + } + } + + for (j = lp1; j < matrixCopy.ColumnCount; j++) + { + var ww = -e[j] / e[lp1]; + for (var ii = lp1; ii < matrixCopy.RowCount; ii++) + { + matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (ww * work[ii])); + } + } + } + + if (ComputeVectors) + { + // Place the transformation in v for subsequent back multiplication. + for (i = lp1; i < matrixCopy.ColumnCount; i++) + { + MatrixVT.At(i, l, e[i]); + } + } + } + + // Set up the final bidiagonal matrixCopy or order m. + var m = Math.Min(matrixCopy.ColumnCount, matrixCopy.RowCount + 1); + var nctp1 = nct + 1; + var nrtp1 = nrt + 1; + if (nct < matrixCopy.ColumnCount) + { + VectorS[nctp1 - 1] = matrixCopy.At((nctp1 - 1), (nctp1 - 1)); + } + + if (matrixCopy.RowCount < m) + { + VectorS[m - 1] = 0.0; + } + + if (nrtp1 < m) + { + e[nrtp1 - 1] = matrixCopy.At((nrtp1 - 1), (m - 1)); + } + + e[m - 1] = 0.0; + + // If required, generate u. + if (ComputeVectors) + { + for (j = nctp1 - 1; j < ncu; j++) + { + for (i = 0; i < matrixCopy.RowCount; i++) + { + MatrixU.At(i, j, 0.0); + } + + MatrixU.At(j, j, 1.0); + } + + for (l = nct - 1; l >= 0; l--) + { + if (VectorS[l] != 0.0) + { + for (j = l + 1; j < ncu; j++) + { + t = -Ddot(MatrixU, matrixCopy.RowCount, l, j, l) / MatrixU.At(l, l); + for (var ii = l; ii < matrixCopy.RowCount; ii++) + { + MatrixU.At(ii, j, MatrixU.At(ii, j) + (t * MatrixU.At(ii, l))); + } + } + + DscalColumn(MatrixU, matrixCopy.RowCount, l, l, -1.0); + MatrixU.At(l, l, 1.0 + MatrixU.At(l, l)); + for (i = 0; i < l; i++) + { + MatrixU.At(i, l, 0.0); + } + } + else + { + for (i = 0; i < matrixCopy.RowCount; i++) + { + MatrixU.At(i, l, 0.0); + } + + MatrixU.At(l, l, 1.0); + } + } + } + + // If it is required, generate v. + if (ComputeVectors) + { + for (l = matrixCopy.ColumnCount - 1; l >= 0; l--) + { + lp1 = l + 1; + if (l < nrt) + { + if (e[l] != 0.0) + { + for (j = lp1; j < matrixCopy.ColumnCount; j++) + { + t = -Ddot(MatrixVT, matrixCopy.ColumnCount, l, j, lp1) / MatrixVT.At(lp1, l); + for (var ii = l; ii < matrixCopy.ColumnCount; ii++) + { + MatrixVT.At(ii, j, MatrixVT.At(ii, j) + (t * MatrixVT.At(ii, l))); + } + } + } + } + + for (i = 0; i < matrixCopy.ColumnCount; i++) + { + MatrixVT.At(i, l, 0.0); + } + + MatrixVT.At(l, l, 1.0); + } + } + + // Transform s and e so that they are double . + for (i = 0; i < m; i++) + { + double r; + if (VectorS[i] != 0.0) + { + t = VectorS[i]; + r = VectorS[i] / t; + VectorS[i] = t; + if (i < m - 1) + { + e[i] = e[i] / r; + } + + if (ComputeVectors) + { + DscalColumn(MatrixU, matrixCopy.RowCount, i, 0, r); + } + } + + // Exit + if (i == m - 1) + { + break; + } + + if (e[i] != 0.0) + { + t = e[i]; + r = t / e[i]; + e[i] = t; + VectorS[i + 1] = VectorS[i + 1] * r; + if (ComputeVectors) + { + DscalColumn(MatrixVT, matrixCopy.ColumnCount, i + 1, 0, r); + } + } + } + + // Main iteration loop for the singular values. + var mn = m; + var iter = 0; + + while (m > 0) + { + // Quit if all the singular values have been found. If too many iterations have been performed, + // throw exception that Convergence Failed + if (iter >= Maxiter) + { + throw new ArgumentException(Resources.ConvergenceFailed); + } + + // This section of the program inspects for negligible elements in the s and e arrays. On + // completion the variables kase and l are set as follows. + // Kase = 1 if VectorS[m] and e[l-1] are negligible and l < m + // Kase = 2 if VectorS[l] is negligible and l < m + // Kase = 3 if e[l-1] is negligible, l < m, and VectorS[l, ..., VectorS[m] are not negligible (qr step). + // Лase = 4 if e[m-1] is negligible (convergence). + double ztest; + double test; + for (l = m - 2; l >= 0; l--) + { + test = Math.Abs(VectorS[l]) + Math.Abs(VectorS[l + 1]); + ztest = test + Math.Abs(e[l]); + if (ztest.AlmostEqualInDecimalPlaces(test, 15)) + { + e[l] = 0.0; + break; + } + } + + int kase; + if (l == m - 2) + { + kase = 4; + } + else + { + int ls; + for (ls = m - 1; ls > l; ls--) + { + test = 0.0; + if (ls != m - 1) + { + test = test + Math.Abs(e[ls]); + } + + if (ls != l + 1) + { + test = test + Math.Abs(e[ls - 1]); + } + + ztest = test + Math.Abs(VectorS[ls]); + if (ztest.AlmostEqualInDecimalPlaces(test, 15)) + { + VectorS[ls] = 0.0; + break; + } + } + + if (ls == l) + { + kase = 3; + } + else if (ls == m - 1) + { + kase = 1; + } + else + { + kase = 2; + l = ls; + } + } + + l = l + 1; + + // Perform the task indicated by kase. + int k; + double f; + switch (kase) + { + // Deflate negligible VectorS[m]. + case 1: + f = e[m - 2]; + e[m - 2] = 0.0; + double t1; + for (var kk = l; kk < m - 1; kk++) + { + k = m - 2 - kk + l; + t1 = VectorS[k]; + Drotg(ref t1, ref f, ref cs, ref sn); + VectorS[k] = t1; + if (k != l) + { + f = -sn * e[k - 1]; + e[k - 1] = cs * e[k - 1]; + } + + if (ComputeVectors) + { + Drot(MatrixVT, matrixCopy.ColumnCount, k, m - 1, cs, sn); + } + } + + break; + + // Split at negligible VectorS[l]. + case 2: + f = e[l - 1]; + e[l - 1] = 0.0; + for (k = l; k < m; k++) + { + t1 = VectorS[k]; + Drotg(ref t1, ref f, ref cs, ref sn); + VectorS[k] = t1; + f = -sn * e[k]; + e[k] = cs * e[k]; + if (ComputeVectors) + { + Drot(MatrixU, matrixCopy.RowCount, k, l - 1, cs, sn); + } + } + + break; + + // Perform one qr step. + case 3: + // Calculate the shift. + var scale = 0.0; + scale = Math.Max(scale, Math.Abs(VectorS[m - 1])); + scale = Math.Max(scale, Math.Abs(VectorS[m - 2])); + scale = Math.Max(scale, Math.Abs(e[m - 2])); + scale = Math.Max(scale, Math.Abs(VectorS[l])); + scale = Math.Max(scale, Math.Abs(e[l])); + var sm = VectorS[m - 1] / scale; + var smm1 = VectorS[m - 2] / scale; + var emm1 = e[m - 2] / scale; + var sl = VectorS[l] / scale; + var el = e[l] / scale; + var b = (((smm1 + sm) * (smm1 - sm)) + (emm1 * emm1)) / 2.0; + var c = (sm * emm1) * (sm * emm1); + var shift = 0.0; + if (b != 0.0 || c != 0.0) + { + shift = Math.Sqrt((b * b) + c); + if (b < 0.0) + { + shift = -shift; + } + + shift = c / (b + shift); + } + + f = ((sl + sm) * (sl - sm)) + shift; + var g = sl * el; + + // Chase zeros. + for (k = l; k < m - 1; k++) + { + Drotg(ref f, ref g, ref cs, ref sn); + if (k != l) + { + e[k - 1] = f; + } + + f = (cs * VectorS[k]) + (sn * e[k]); + e[k] = (cs * e[k]) - (sn * VectorS[k]); + g = sn * VectorS[k + 1]; + VectorS[k + 1] = cs * VectorS[k + 1]; + if (ComputeVectors) + { + Drot(MatrixVT, matrixCopy.ColumnCount, k, k + 1, cs, sn); + } + + Drotg(ref f, ref g, ref cs, ref sn); + VectorS[k] = f; + f = (cs * e[k]) + (sn * VectorS[k + 1]); + VectorS[k + 1] = (-sn * e[k]) + (cs * VectorS[k + 1]); + g = sn * e[k + 1]; + e[k + 1] = cs * e[k + 1]; + if (ComputeVectors && k < matrixCopy.RowCount) + { + Drot(MatrixU, matrixCopy.RowCount, k, k + 1, cs, sn); + } + } + + e[m - 2] = f; + iter = iter + 1; + break; + + // Convergence. + case 4: + // Make the singular value positive + if (VectorS[l] < 0.0) + { + VectorS[l] = -VectorS[l]; + if (ComputeVectors) + { + DscalColumn(MatrixVT, matrixCopy.ColumnCount, l, 0, -1.0); + } + } + + // Order the singular value. + while (l != mn - 1) + { + if (VectorS[l] >= VectorS[l + 1]) + { + break; + } + + t = VectorS[l]; + VectorS[l] = VectorS[l + 1]; + VectorS[l + 1] = t; + if (ComputeVectors && l < matrixCopy.ColumnCount) + { + Dswap(MatrixVT, matrixCopy.ColumnCount, l, l + 1); + } + + if (ComputeVectors && l < matrixCopy.RowCount) + { + Dswap(MatrixU, matrixCopy.RowCount, l, l + 1); + } + + l = l + 1; + } + + iter = 0; + m = m - 1; + break; + } + } + + if (ComputeVectors) + { + MatrixVT = MatrixVT.Transpose(); + } + + // Adjust the size of s if rows < columns. We are using ported copy of linpack's svd code and it uses + // a singular vector of length mRows+1 when mRows < mColumns. The last element is not used and needs to be removed. + // we should port lapack's svd routine to remove this problem. + if (matrixCopy.RowCount < matrixCopy.ColumnCount) + { + nm--; + var tmp = matrixCopy.CreateVector(nm); + for (i = 0; i < nm; i++) + { + tmp[i] = VectorS[i]; + } + + VectorS = tmp; + } + } + + /// + /// Calculates absolute value of multiplied on signum function of + /// + /// Double value z1 + /// Double value z2 + /// Result multiplication of signum function and absolute value + private static double Dsign(double z1, double z2) + { + return Math.Abs(z1) * (z2 / Math.Abs(z2)); + } + + /// + /// Swap column and + /// + /// Source matrix + /// The number of rows in + /// Column A index to swap + /// Column B index to swap + private static void Dswap(Matrix a, int rowCount, int columnA, int columnB) + { + for (var i = 0; i < rowCount; i++) + { + var z = a.At(i, columnA); + a.At(i, columnA, a.At(i, columnB)); + a.At(i, columnB, z); + } + } + + /// + /// Scale column by starting from row + /// + /// Source matrix + /// The number of rows in + /// Column to scale + /// Row to scale from + /// Scale value + private static void DscalColumn(Matrix a, int rowCount, int column, int rowStart, double z) + { + for (var i = rowStart; i < rowCount; i++) + { + a.At(i, column, a.At(i, column) * z); + } + } + + /// + /// Scale vector by starting from index + /// + /// Source vector + /// Row to scale from + /// Scale value + private static void DscalVector(double[] a, int start, double z) + { + for (var i = start; i < a.Length; i++) + { + a[i] = a[i] * z; + } + } + + /// + /// Given the Cartesian coordinates (da, db) of a point p, these fucntion return the parameters da, db, c, and s + /// associated with the Givens rotation that zeros the y-coordinate of the point. + /// + /// Provides the x-coordinate of the point p. On exit contains the parameter r associated with the Givens rotation + /// Provides the y-coordinate of the point p. On exit contains the parameter z associated with the Givens rotation + /// Contains the parameter c associated with the Givens rotation + /// Contains the parameter s associated with the Givens rotation + /// This is equivalent to the DROTG LAPACK routine. + private static void Drotg(ref double da, ref double db, ref double c, ref double s) + { + double r, z; + + var roe = db; + var absda = Math.Abs(da); + var absdb = Math.Abs(db); + if (absda > absdb) + { + roe = da; + } + + var scale = absda + absdb; + if (scale == 0.0) + { + c = 1.0; + s = 0.0; + r = 0.0; + z = 0.0; + } + else + { + var sda = da / scale; + var sdb = db / scale; + r = scale * Math.Sqrt((sda * sda) + (sdb * sdb)); + if (roe < 0.0) + { + r = -r; + } + + c = da / r; + s = db / r; + z = 1.0; + if (absda > absdb) + { + z = s; + } + + if (absdb >= absda && c != 0.0) + { + z = 1.0 / c; + } + } + + da = r; + db = z; + } + + /// dded + /// Calculate Norm 2 of the column in matrix starting from row + /// + /// Source matrix + /// The number of rows in + /// Column index + /// Start row index + /// Norm2 (Euclidean norm) of trhe column + private static double Dnrm2Column(Matrix a, int rowCount, int column, int rowStart) + { + double s = 0; + for (var i = rowStart; i < rowCount; i++) + { + s += a.At(i, column) * a.At(i, column); + } + + return Math.Sqrt(s); + } + + /// + /// Calculate Norm 2 of the vector starting from index + /// + /// Source vector + /// Start index + /// Norm2 (Euclidean norm) of the vector + private static double Dnrm2Vector(double[] a, int rowStart) + { + double s = 0; + for (var i = rowStart; i < a.Length; i++) + { + s += a[i] * a[i]; + } + + return Math.Sqrt(s); + } + + /// + /// Calculate dot product of and + /// + /// Source matrix + /// The number of rows in + /// Index of column A + /// Index of column B + /// Starting row index + /// Dot product value + private static double Ddot(Matrix a, int rowCount, int columnA, int columnB, int rowStart) + { + var z = 0.0; + for (var i = rowStart; i < rowCount; i++) + { + z += a.At(i, columnB) * a.At(i, columnA); + } + + return z; + } + + /// + /// Performs rotation of points in the plane. Given two vectors x and y , + /// each vector element of these vectors is replaced as follows: x(i) = c*x(i) + s*y(i); y(i) = c*y(i) - s*x(i) + /// + /// Source matrix + /// The number of rows in + /// Index of column A + /// Index of column B + /// Scalar "c" value + /// Scalar "s" value + private static void Drot(Matrix a, int rowCount, int columnA, int columnB, double c, double s) + { + for (var i = 0; i < rowCount; i++) + { + var z = (c * a.At(i, columnA)) + (s * a.At(i, columnB)); + var tmp = (c * a.At(i, columnB)) - (s * a.At(i, columnA)); + a.At(i, columnB, tmp); + a.At(i, columnA, z); + } + } + + /// + /// Solves a system of linear equations, AX = B, with A SVD factorized. + /// + /// The right hand side , B. + /// The left hand side , X. + 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 mn = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount); + var bn = input.ColumnCount; + + var tmp = new double[MatrixVT.ColumnCount]; + + for (var k = 0; k < bn; k++) + { + for (var j = 0; j < MatrixVT.ColumnCount; j++) + { + double value = 0; + if (j < mn) + { + for (var i = 0; i < MatrixU.RowCount; i++) + { + value += MatrixU.At(i, j) * input.At(i, k); + } + + value /= VectorS[j]; + } + + tmp[j] = value; + } + + for (var j = 0; j < MatrixVT.ColumnCount; j++) + { + double value = 0; + for (var i = 0; i < MatrixVT.ColumnCount; i++) + { + value += MatrixVT.At(i, j) * tmp[i]; + } + + result[j, k] = value; + } + } + } + + /// + /// Solves a system of linear equations, Ax = b, with A SVD factorized. + /// + /// The right hand side vector, b. + /// The left hand side , x. + 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 mn = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount); + var tmp = new double[MatrixVT.ColumnCount]; + double value; + for (var j = 0; j < MatrixVT.ColumnCount; j++) + { + value = 0; + if (j < mn) + { + for (var i = 0; i < MatrixU.RowCount; i++) + { + value += MatrixU.At(i, j) * input[i]; + } + + value /= VectorS[j]; + } + + tmp[j] = value; + } + + for (var j = 0; j < MatrixVT.ColumnCount; j++) + { + value = 0; + for (int i = 0; i < MatrixVT.ColumnCount; i++) + { + value += MatrixVT.At(i, j) * tmp[i]; + } + + result[j] = value; + } + } + } +} diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 7979e567..84eec40a 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -103,6 +103,10 @@ + + + + diff --git a/src/Silverlight/Silverlight.csproj b/src/Silverlight/Silverlight.csproj index 369daafb..7de94595 100644 --- a/src/Silverlight/Silverlight.csproj +++ b/src/Silverlight/Silverlight.csproj @@ -248,6 +248,18 @@ LinearAlgebra\Double\Factorization\Svd.cs + + LinearAlgebra\Double\Factorization\UserCholesky.cs + + + LinearAlgebra\Double\Factorization\UserLU.cs + + + LinearAlgebra\Double\Factorization\UserQR.cs + + + LinearAlgebra\Double\Factorization\UserSvd.cs + LinearAlgebra\Double\ISolver.cs diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/CholeskyTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/CholeskyTests.cs index 99cc3bc1..f7d98a4f 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Factorization/CholeskyTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/CholeskyTests.cs @@ -30,7 +30,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization { - using System.Collections.Generic; using MbUnit.Framework; using LinearAlgebra.Double; using LinearAlgebra.Double.Factorization; @@ -44,23 +43,16 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization public void CanFactorizeIdentity(int order) { var I = DenseMatrix.Identity(order); - var C = I.Cholesky(); + var factorC = I.Cholesky(); - Assert.AreEqual(I.RowCount, C.Factor.RowCount); - Assert.AreEqual(I.ColumnCount, C.Factor.ColumnCount); + Assert.AreEqual(I.RowCount, factorC.Factor.RowCount); + Assert.AreEqual(I.ColumnCount, factorC.Factor.ColumnCount); - for (var i = 0; i < C.Factor.RowCount; i++) + for (var i = 0; i < factorC.Factor.RowCount; i++) { - for (var j = 0; j < C.Factor.ColumnCount; j++) + for (var j = 0; j < factorC.Factor.ColumnCount; j++) { - if (i == j) - { - Assert.AreEqual(1.0, C.Factor[i, j]); - } - else - { - Assert.AreEqual(0.0, C.Factor[i, j]); - } + Assert.AreEqual(i == j ? 1.0 : 0.0, factorC.Factor[i, j]); } } } @@ -71,7 +63,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization { var I = DenseMatrix.Identity(10); I[3, 3] = -4.0; - var C = I.Cholesky(); + I.Cholesky(); } [Test] @@ -81,7 +73,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization public void CholeskyFailsWithNonSquareMatrix(int row, int col) { var I = new DenseMatrix(row, col); - var C = I.Cholesky(); + I.Cholesky(); } [Test] @@ -91,9 +83,9 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization 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); + var factorC = I.Cholesky(); + Assert.AreEqual(1.0, factorC.Determinant); + Assert.AreEqual(0.0, factorC.DeterminantLn); } [Test] @@ -106,30 +98,30 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanFactorizeRandomMatrix(int order) { - var X = MatrixLoader.GenerateRandomPositiveDefiniteMatrix(order); - var chol = X.Cholesky(); - var C = chol.Factor; + var matrixX = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order); + var chol = matrixX.Cholesky(); + var factorC = chol.Factor; // Make sure the Cholesky factor has the right dimensions. - Assert.AreEqual(order, C.RowCount); - Assert.AreEqual(order, C.ColumnCount); + Assert.AreEqual(order, factorC.RowCount); + Assert.AreEqual(order, factorC.ColumnCount); // Make sure the Cholesky factor is lower triangular. - for (int i = 0; i < C.RowCount; i++) + for (var i = 0; i < factorC.RowCount; i++) { - for (int j = i+1; j < C.ColumnCount; j++) + for (var j = i+1; j < factorC.ColumnCount; j++) { - Assert.AreEqual(0.0, C[i, j]); + Assert.AreEqual(0.0, factorC[i, j]); } } // Make sure the cholesky factor times it's transpose is the original matrix. - var XfromC = C * C.Transpose(); - for (int i = 0; i < XfromC.RowCount; i++) + var matrixXfromC = factorC * factorC.Transpose(); + for (var i = 0; i < matrixXfromC.RowCount; i++) { - for (int j = 0; j < XfromC.ColumnCount; j++) + for (var j = 0; j < matrixXfromC.ColumnCount; j++) { - Assert.AreApproximatelyEqual(X[i,j], XfromC[i, j], 1.0e-11); + Assert.AreApproximatelyEqual(matrixX[i,j], matrixXfromC[i, j], 1.0e-11); } } } @@ -144,28 +136,28 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanSolveForRandomVector(int order) { - var A = MatrixLoader.GenerateRandomPositiveDefiniteMatrix(order); - var ACopy = A.Clone(); - var chol = A.Cholesky(); - var b = MatrixLoader.GenerateRandomVector(order); + var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order); + var matrixACopy = matrixA.Clone(); + var chol = matrixA.Cholesky(); + var b = MatrixLoader.GenerateRandomDenseVector(order); var x = chol.Solve(b); Assert.AreEqual(b.Count, x.Count); - var bReconstruct = A * x; + var bReconstruct = matrixA * x; // Check the reconstruction. - for (int i = 0; i < order; i++) + for (var i = 0; i < order; i++) { Assert.AreApproximatelyEqual(b[i], bReconstruct[i], 1.0e-11); } // Make sure A didn't change. - for (int i = 0; i < A.RowCount; i++) + for (var i = 0; i < matrixA.RowCount; i++) { - for (int j = 0; j < A.ColumnCount; j++) + for (var j = 0; j < matrixA.ColumnCount; j++) { - Assert.AreEqual(ACopy[i, j], A[i, j]); + Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); } } } @@ -180,32 +172,32 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanSolveForRandomMatrix(int row, int col) { - var A = MatrixLoader.GenerateRandomPositiveDefiniteMatrix(row); - var ACopy = A.Clone(); - var chol = A.Cholesky(); - var B = MatrixLoader.GenerateRandomMatrix(row, col); - var X = chol.Solve(B); + var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(row); + var matrixACopy = matrixA.Clone(); + var chol = matrixA.Cholesky(); + var matrixB = MatrixLoader.GenerateRandomDenseMatrix(row, col); + var matrixX = chol.Solve(matrixB); - Assert.AreEqual(B.RowCount, X.RowCount); - Assert.AreEqual(B.ColumnCount, X.ColumnCount); + Assert.AreEqual(matrixB.RowCount, matrixX.RowCount); + Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); - var BReconstruct = A * X; + var matrixBReconstruct = matrixA * matrixX; // Check the reconstruction. - for (int i = 0; i < B.RowCount; i++) + for (var i = 0; i < matrixB.RowCount; i++) { - for (int j = 0; j < B.ColumnCount; j++) + for (var j = 0; j < matrixB.ColumnCount; j++) { - Assert.AreApproximatelyEqual(B[i, j], BReconstruct[i, j], 1.0e-11); + Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11); } } // Make sure A didn't change. - for (int i = 0; i < A.RowCount; i++) + for (var i = 0; i < matrixA.RowCount; i++) { - for (int j = 0; j < A.ColumnCount; j++) + for (var j = 0; j < matrixA.ColumnCount; j++) { - Assert.AreEqual(ACopy[i, j], A[i, j]); + Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); } } } @@ -220,35 +212,35 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanSolveForRandomVectorWhenResultVectorGiven(int order) { - var A = MatrixLoader.GenerateRandomPositiveDefiniteMatrix(order); - var ACopy = A.Clone(); - var chol = A.Cholesky(); - var b = MatrixLoader.GenerateRandomVector(order); + var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(order); + var matrixACopy = matrixA.Clone(); + var chol = matrixA.Cholesky(); + var b = MatrixLoader.GenerateRandomDenseVector(order); var bCopy = b.Clone(); var x = new DenseVector(order); chol.Solve(b, x); Assert.AreEqual(b.Count, x.Count); - var bReconstruct = A * x; + var bReconstruct = matrixA * x; // Check the reconstruction. - for (int i = 0; i < order; i++) + for (var i = 0; i < order; i++) { Assert.AreApproximatelyEqual(b[i], bReconstruct[i], 1.0e-11); } // Make sure A didn't change. - for (int i = 0; i < A.RowCount; i++) + for (var i = 0; i < matrixA.RowCount; i++) { - for (int j = 0; j < A.ColumnCount; j++) + for (var j = 0; j < matrixA.ColumnCount; j++) { - Assert.AreEqual(ACopy[i, j], A[i, j]); + Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); } } // Make sure b didn't change. - for (int i = 0; i < order; i++) + for (var i = 0; i < order; i++) { Assert.AreEqual(bCopy[i], b[i]); } @@ -264,43 +256,43 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanSolveForRandomMatrixWhenResultMatrixGiven(int row, int col) { - var A = MatrixLoader.GenerateRandomPositiveDefiniteMatrix(row); - var ACopy = A.Clone(); - var chol = A.Cholesky(); - var B = MatrixLoader.GenerateRandomMatrix(row, col); - var BCopy = B.Clone(); - var X = new DenseMatrix(row, col); - chol.Solve(B, X); + var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteDenseMatrix(row); + var matrixACopy = matrixA.Clone(); + var chol = matrixA.Cholesky(); + var matrixB = MatrixLoader.GenerateRandomDenseMatrix(row, col); + var matrixBCopy = matrixB.Clone(); + var matrixX = new DenseMatrix(row, col); + chol.Solve(matrixB, matrixX); - Assert.AreEqual(B.RowCount, X.RowCount); - Assert.AreEqual(B.ColumnCount, X.ColumnCount); + Assert.AreEqual(matrixB.RowCount, matrixX.RowCount); + Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); - var BReconstruct = A * X; + var matrixBReconstruct = matrixA * matrixX; // Check the reconstruction. - for (int i = 0; i < B.RowCount; i++) + for (var i = 0; i < matrixB.RowCount; i++) { - for (int j = 0; j < B.ColumnCount; j++) + for (var j = 0; j < matrixB.ColumnCount; j++) { - Assert.AreApproximatelyEqual(B[i, j], BReconstruct[i, j], 1.0e-11); + Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11); } } // Make sure A didn't change. - for (int i = 0; i < A.RowCount; i++) + for (var i = 0; i < matrixA.RowCount; i++) { - for (int j = 0; j < A.ColumnCount; j++) + for (var j = 0; j < matrixA.ColumnCount; j++) { - Assert.AreEqual(ACopy[i, j], A[i, j]); + Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); } } // Make sure B didn't change. - for (int i = 0; i < B.RowCount; i++) + for (var i = 0; i < matrixB.RowCount; i++) { - for (int j = 0; j < B.ColumnCount; j++) + for (var j = 0; j < matrixB.ColumnCount; j++) { - Assert.AreEqual(BCopy[i, j], B[i, j]); + Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]); } } } diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/LUTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/LUTests.cs index 2fb150ca..3c8da2af 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Factorization/LUTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/LUTests.cs @@ -101,7 +101,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanFactorizeRandomMatrix(int order) { - var matrixX = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixX = MatrixLoader.GenerateRandomDenseMatrix(order, order); var factorLU = matrixX.LU(); var matrixL = factorLU.L; var matrixU = factorLU.U; @@ -154,11 +154,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanSolveForRandomVector(int order) { - var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); var matrixACopy = matrixA.Clone(); var factorLU = matrixA.LU(); - var vectorb = MatrixLoader.GenerateRandomVector(order); + var vectorb = MatrixLoader.GenerateRandomDenseVector(order); var resultx = factorLU.Solve(vectorb); Assert.AreEqual(matrixA.ColumnCount, resultx.Count); @@ -191,11 +191,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanSolveForRandomMatrix(int order) { - var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); var matrixACopy = matrixA.Clone(); var factorLU = matrixA.LU(); - var matrixB = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order); var matrixX = factorLU.Solve(matrixB); // The solution X row dimension is equal to the column dimension of A @@ -234,10 +234,10 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanSolveForRandomVectorWhenResultVectorGiven(int order) { - var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); var matrixACopy = matrixA.Clone(); var factorLU = matrixA.LU(); - var vectorb = MatrixLoader.GenerateRandomVector(order); + var vectorb = MatrixLoader.GenerateRandomDenseVector(order); var vectorbCopy = vectorb.Clone(); var resultx = new DenseVector(order); factorLU.Solve(vectorb, resultx); @@ -278,11 +278,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanSolveForRandomMatrixWhenResultMatrixGiven(int order) { - var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); var matrixACopy = matrixA.Clone(); var factorLU = matrixA.LU(); - var matrixB = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order); var matrixBCopy = matrixB.Clone(); var matrixX = new DenseMatrix(order, order); @@ -333,7 +333,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanInverse(int order) { - var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); var matrixACopy = matrixA.Clone(); var factorLU = matrixA.LU(); diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/QRTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/QRTests.cs index 51088763..ac08952d 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Factorization/QRTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/QRTests.cs @@ -101,7 +101,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanFactorizeRandomMatrix(int row, int column) { - var matrixA = MatrixLoader.GenerateRandomMatrix(row, column); + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(row, column); var factorQR = matrixA.QR(); // Make sure the R has the right dimensions. @@ -145,11 +145,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanSolveForRandomVector(int order) { - var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); var matrixACopy = matrixA.Clone(); var factorQR = matrixA.QR(); - var vectorb = MatrixLoader.GenerateRandomVector(order); + var vectorb = MatrixLoader.GenerateRandomDenseVector(order); var resultx = factorQR.Solve(vectorb); Assert.AreEqual(matrixA.ColumnCount, resultx.Count); @@ -182,11 +182,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanSolveForRandomMatrix(int order) { - var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); var matrixACopy = matrixA.Clone(); var factorQR = matrixA.QR(); - var matrixB = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order); var matrixX = factorQR.Solve(matrixB); // The solution X row dimension is equal to the column dimension of A @@ -225,10 +225,10 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanSolveForRandomVectorWhenResultVectorGiven(int order) { - var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); var matrixACopy = matrixA.Clone(); var factorQR = matrixA.QR(); - var vectorb = MatrixLoader.GenerateRandomVector(order); + var vectorb = MatrixLoader.GenerateRandomDenseVector(order); var vectorbCopy = vectorb.Clone(); var resultx = new DenseVector(order); factorQR.Solve(vectorb,resultx); @@ -269,11 +269,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanSolveForRandomMatrixWhenResultMatrixGiven(int order) { - var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); var matrixACopy = matrixA.Clone(); var factorQR = matrixA.QR(); - var matrixB = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order); var matrixBCopy = matrixB.Clone(); var matrixX = new DenseMatrix(order, order); diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/SvdTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/SvdTests.cs index fdc9bb12..fa35945d 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Factorization/SvdTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/SvdTests.cs @@ -82,7 +82,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanFactorizeRandomMatrix(int row, int column) { - var matrixA = MatrixLoader.GenerateRandomMatrix(row, column); + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(row, column); var factorSvd = matrixA.Svd(true); // Make sure the U has the right dimensions. @@ -115,7 +115,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CheckRankOfNonSquare(int row, int column) { - var matrixA = MatrixLoader.GenerateRandomMatrix(row, column); + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(row, column); var factorSvd = matrixA.Svd(true); var mn = Math.Min(row, column); @@ -132,7 +132,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CheckRankSquare(int order) { - var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); var factorSvd = matrixA.Svd(true); if (factorSvd.Determinant != 0) @@ -172,10 +172,10 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [ExpectedException(typeof(InvalidOperationException))] public void CannotSolveMatrixIfVectorsNotComputed() { - var matrixA = MatrixLoader.GenerateRandomMatrix(10, 10); + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(10, 10); var factorSvd = matrixA.Svd(false); - var matrixB = MatrixLoader.GenerateRandomMatrix(10, 10); + var matrixB = MatrixLoader.GenerateRandomDenseMatrix(10, 10); factorSvd.Solve(matrixB); } @@ -183,10 +183,10 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [ExpectedException(typeof(InvalidOperationException))] public void CannotSolveVectorIfVectorsNotComputed() { - var matrixA = MatrixLoader.GenerateRandomMatrix(10, 10); + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(10, 10); var factorSvd = matrixA.Svd(false); - var vectorb = MatrixLoader.GenerateRandomVector(10); + var vectorb = MatrixLoader.GenerateRandomDenseVector(10); factorSvd.Solve(vectorb); } @@ -200,11 +200,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanSolveForRandomVector(int row, int column) { - var matrixA = MatrixLoader.GenerateRandomMatrix(row, column); + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(row, column); var matrixACopy = matrixA.Clone(); var factorSvd = matrixA.Svd(true); - var vectorb = MatrixLoader.GenerateRandomVector(row); + var vectorb = MatrixLoader.GenerateRandomDenseVector(row); var resultx = factorSvd.Solve(vectorb); Assert.AreEqual(matrixA.ColumnCount, resultx.Count); @@ -237,11 +237,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanSolveForRandomMatrix(int row, int count) { - var matrixA = MatrixLoader.GenerateRandomMatrix(row, count); + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(row, count); var matrixACopy = matrixA.Clone(); var factorSvd = matrixA.Svd(true); - var matrixB = MatrixLoader.GenerateRandomMatrix(row, count); + var matrixB = MatrixLoader.GenerateRandomDenseMatrix(row, count); var matrixX = factorSvd.Solve(matrixB); // The solution X row dimension is equal to the column dimension of A @@ -280,10 +280,10 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanSolveForRandomVectorWhenResultVectorGiven(int row, int column) { - var matrixA = MatrixLoader.GenerateRandomMatrix(row, column); + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(row, column); var matrixACopy = matrixA.Clone(); var factorSvd = matrixA.Svd(true); - var vectorb = MatrixLoader.GenerateRandomVector(row); + var vectorb = MatrixLoader.GenerateRandomDenseVector(row); var vectorbCopy = vectorb.Clone(); var resultx = new DenseVector(column); factorSvd.Solve(vectorb,resultx); @@ -322,11 +322,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization [MultipleAsserts] public void CanSolveForRandomMatrixWhenResultMatrixGiven(int row, int column) { - var matrixA = MatrixLoader.GenerateRandomMatrix(row, column); + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(row, column); var matrixACopy = matrixA.Clone(); var factorSvd = matrixA.Svd(true); - var matrixB = MatrixLoader.GenerateRandomMatrix(row, column); + var matrixB = MatrixLoader.GenerateRandomDenseMatrix(row, column); var matrixBCopy = matrixB.Clone(); var matrixX = new DenseMatrix(column, column); diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserCholeskyTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserCholeskyTests.cs new file mode 100644 index 00000000..9840babb --- /dev/null +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserCholeskyTests.cs @@ -0,0 +1,299 @@ +// +// 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. +// + +namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization +{ + using MbUnit.Framework; + using LinearAlgebra.Double.Factorization; + + public class UserCholeskyTests + { + [Test] + [Row(1)] + [Row(10)] + [Row(100)] + public void CanFactorizeIdentity(int order) + { + var I = UserDefinedMatrix.Identity(order); + var factorC = I.Cholesky(); + + Assert.AreEqual(I.RowCount, factorC.Factor.RowCount); + Assert.AreEqual(I.ColumnCount, factorC.Factor.ColumnCount); + + for (var i = 0; i < factorC.Factor.RowCount; i++) + { + for (var j = 0; j < factorC.Factor.ColumnCount; j++) + { + Assert.AreEqual(i == j ? 1.0 : 0.0, factorC.Factor[i, j]); + } + } + } + + [Test] + [ExpectedArgumentException] + public void CholeskyFailsWithDiagonalNonPositiveDefiniteMatrix() + { + var I = UserDefinedMatrix.Identity(10); + I[3, 3] = -4.0; + I.Cholesky(); + } + + [Test] + [Row(3,5)] + [Row(5,3)] + [ExpectedArgumentException] + public void CholeskyFailsWithNonSquareMatrix(int row, int col) + { + var I = new UserDefinedMatrix(row, col); + I.Cholesky(); + } + + [Test] + [Row(1)] + [Row(10)] + [Row(100)] + public void IdentityDeterminantIsOne(int order) + { + var I = UserDefinedMatrix.Identity(order); + var factorC = I.Cholesky(); + Assert.AreEqual(1.0, factorC.Determinant); + Assert.AreEqual(0.0, factorC.DeterminantLn); + } + + [Test] + [Row(1)] + [Row(2)] + [Row(5)] + [Row(10)] + [Row(50)] + [Row(100)] + [MultipleAsserts] + public void CanFactorizeRandomMatrix(int order) + { + var matrixX = MatrixLoader.GenerateRandomPositiveDefiniteUserDefinedMatrix(order); + var chol = matrixX.Cholesky(); + var factorC = chol.Factor; + + // Make sure the Cholesky factor has the right dimensions. + Assert.AreEqual(order, factorC.RowCount); + Assert.AreEqual(order, factorC.ColumnCount); + + // Make sure the Cholesky factor is lower triangular. + for (var i = 0; i < factorC.RowCount; i++) + { + for (var j = i+1; j < factorC.ColumnCount; j++) + { + Assert.AreEqual(0.0, factorC[i, j]); + } + } + + // Make sure the cholesky factor times it's transpose is the original matrix. + var matrixXfromC = factorC * factorC.Transpose(); + for (var i = 0; i < matrixXfromC.RowCount; i++) + { + for (var j = 0; j < matrixXfromC.ColumnCount; j++) + { + Assert.AreApproximatelyEqual(matrixX[i,j], matrixXfromC[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.GenerateRandomPositiveDefiniteUserDefinedMatrix(order); + var matrixACopy = matrixA.Clone(); + var chol = matrixA.Cholesky(); + var b = MatrixLoader.GenerateRandomUserDefinedVector(order); + var x = chol.Solve(b); + + Assert.AreEqual(b.Count, x.Count); + + var bReconstruct = matrixA * x; + + // Check the reconstruction. + for (var i = 0; i < order; i++) + { + Assert.AreApproximatelyEqual(b[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(2,4)] + [Row(5,8)] + [Row(10,3)] + [Row(50,10)] + [Row(100,100)] + [MultipleAsserts] + public void CanSolveForRandomMatrix(int row, int col) + { + var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteUserDefinedMatrix(row); + var matrixACopy = matrixA.Clone(); + var chol = matrixA.Cholesky(); + var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(row, col); + var matrixX = chol.Solve(matrixB); + + Assert.AreEqual(matrixB.RowCount, matrixX.RowCount); + 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.GenerateRandomPositiveDefiniteUserDefinedMatrix(order); + var matrixACopy = matrixA.Clone(); + var chol = matrixA.Cholesky(); + var b = MatrixLoader.GenerateRandomUserDefinedVector(order); + var bCopy = b.Clone(); + var x = new UserDefinedVector(order); + chol.Solve(b, x); + + Assert.AreEqual(b.Count, x.Count); + + var bReconstruct = matrixA * x; + + // Check the reconstruction. + for (var i = 0; i < order; i++) + { + Assert.AreApproximatelyEqual(b[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 < order; i++) + { + Assert.AreEqual(bCopy[i], b[i]); + } + } + + [Test] + [Row(1, 1)] + [Row(2, 4)] + [Row(5, 8)] + [Row(10, 3)] + [Row(50, 10)] + [Row(100, 100)] + [MultipleAsserts] + public void CanSolveForRandomMatrixWhenResultMatrixGiven(int row, int col) + { + var matrixA = MatrixLoader.GenerateRandomPositiveDefiniteUserDefinedMatrix(row); + var matrixACopy = matrixA.Clone(); + var chol = matrixA.Cholesky(); + var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(row, col); + var matrixBCopy = matrixB.Clone(); + var matrixX = new UserDefinedMatrix(row, col); + chol.Solve(matrixB, matrixX); + + Assert.AreEqual(matrixB.RowCount, matrixX.RowCount); + 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]); + } + } + } + } +} diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserLUTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserLUTests.cs new file mode 100644 index 00000000..21d42a5e --- /dev/null +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserLUTests.cs @@ -0,0 +1,363 @@ +// +// 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. +// + +namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization +{ + using MbUnit.Framework; + using LinearAlgebra.Double.Factorization; + + public class UserLUTests + { + [Test] + [Row(1)] + [Row(10)] + [Row(100)] + public void CanFactorizeIdentity(int order) + { + var matrixI = UserDefinedMatrix.Identity(order); + var factorLU = matrixI.LU(); + + // Check lower triangular part. + var matrixL = factorLU.L; + Assert.AreEqual(matrixI.RowCount, matrixL.RowCount); + Assert.AreEqual(matrixI.ColumnCount, matrixL.ColumnCount); + for (var i = 0; i < matrixL.RowCount; i++) + { + for (var j = 0; j < matrixL.ColumnCount; j++) + { + Assert.AreEqual(i == j ? 1.0 : 0.0, matrixL[i, j]); + } + } + + // Check upper triangular part. + var matrixU = factorLU.U; + Assert.AreEqual(matrixI.RowCount, matrixU.RowCount); + Assert.AreEqual(matrixI.ColumnCount, matrixU.ColumnCount); + for (var i = 0; i < matrixU.RowCount; i++) + { + for (var j = 0; j < matrixU.ColumnCount; j++) + { + Assert.AreEqual(i == j ? 1.0 : 0.0, matrixU[i, j]); + } + } + } + + [Test] + [Row(3,5)] + [Row(5,3)] + [ExpectedArgumentException] + public void LUFailsWithNonSquareMatrix(int row, int col) + { + var I = new UserDefinedMatrix(row, col); + I.LU(); + } + + [Test] + [Row(1)] + [Row(10)] + [Row(100)] + public void IdentityDeterminantIsOne(int order) + { + var I = UserDefinedMatrix.Identity(order); + var lu = I.LU(); + Assert.AreEqual(1.0, lu.Determinant); + } + + [Test] + [Row(1)] + [Row(2)] + [Row(5)] + [Row(10)] + [Row(50)] + [Row(100)] + [MultipleAsserts] + public void CanFactorizeRandomMatrix(int order) + { + var matrixX = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); + var factorLU = matrixX.LU(); + var matrixL = factorLU.L; + var matrixU = factorLU.U; + + // Make sure the factors have the right dimensions. + Assert.AreEqual(order, matrixL.RowCount); + Assert.AreEqual(order, matrixL.ColumnCount); + Assert.AreEqual(order, matrixU.RowCount); + Assert.AreEqual(order, matrixU.ColumnCount); + + // Make sure the L factor is lower triangular. + for (var i = 0; i < matrixL.RowCount; i++) + { + Assert.AreEqual(1.0, matrixL[i, i]); + for (var j = i+1; j < matrixL.ColumnCount; j++) + { + Assert.AreEqual(0.0, matrixL[i, j]); + } + } + + // Make sure the U factor is upper triangular. + for (var i = 0; i < matrixL.RowCount; i++) + { + for (var j = 0; j < i; j++) + { + Assert.AreEqual(0.0, matrixU[i, j]); + } + } + + // Make sure the LU factor times it's transpose is the original matrix. + var matrixXfromLU = matrixL * matrixU; + var permutationInverse = factorLU.P.Inverse(); + matrixXfromLU.PermuteRows(permutationInverse); + for (var i = 0; i < matrixXfromLU.RowCount; i++) + { + for (var j = 0; j < matrixXfromLU.ColumnCount; j++) + { + Assert.AreApproximatelyEqual(matrixX[i, j], matrixXfromLU[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.GenerateRandomUserDefinedMatrix(order, order); + var matrixACopy = matrixA.Clone(); + var factorLU = matrixA.LU(); + + var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(order); + var resultx = factorLU.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.GenerateRandomUserDefinedMatrix(order, order); + var matrixACopy = matrixA.Clone(); + var factorLU = matrixA.LU(); + + var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); + var matrixX = factorLU.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.GenerateRandomUserDefinedMatrix(order, order); + var matrixACopy = matrixA.Clone(); + var factorLU = matrixA.LU(); + var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(order); + var vectorbCopy = vectorb.Clone(); + var resultx = new UserDefinedVector(order); + factorLU.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.GenerateRandomUserDefinedMatrix(order, order); + var matrixACopy = matrixA.Clone(); + var factorLU = matrixA.LU(); + + var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); + var matrixBCopy = matrixB.Clone(); + + var matrixX = new UserDefinedMatrix(order, order); + factorLU.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]); + } + } + } + + [Test] + [Row(1)] + [Row(4)] + [Row(8)] + [Row(10)] + [Row(50)] + [Row(100)] + [MultipleAsserts] + public void CanInverse(int order) + { + var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); + var matrixACopy = matrixA.Clone(); + var factorLU = matrixA.LU(); + + var matrixAInverse = factorLU.Inverse(); + + // The inverse dimension is equal A + Assert.AreEqual(matrixAInverse.RowCount, matrixAInverse.RowCount); + Assert.AreEqual(matrixAInverse.ColumnCount, matrixAInverse.ColumnCount); + + var matrixIdentity = matrixA * matrixAInverse; + + // 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]); + } + } + + // Check if multiplication of A and AI produced identity matrix. + for (var i = 0; i < matrixIdentity.RowCount; i++) + { + Assert.AreApproximatelyEqual(matrixIdentity[i, i], 1.0, 1.0e-11); + } + } + } +} diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserQRTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserQRTests.cs new file mode 100644 index 00000000..b6a39f30 --- /dev/null +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserQRTests.cs @@ -0,0 +1,316 @@ +// +// 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. +// + +namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization +{ + using MbUnit.Framework; + using LinearAlgebra.Double.Factorization; + + public class UserQRTests + { + + [Test] + [ExpectedArgumentNullException] + public void ConstructorNull() + { + new UserQR(null); + } + + [Test] + [ExpectedArgumentException] + public void WideMatrixThrowsInvalidMatrixOperationException() + { + new UserQR(new UserDefinedMatrix(3, 4)); + } + + [Test] + [Row(1)] + [Row(10)] + [Row(100)] + public void CanFactorizeIdentity(int order) + { + var I = UserDefinedMatrix.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 = UserDefinedMatrix.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.GenerateRandomUserDefinedMatrix(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.GenerateRandomUserDefinedMatrix(order, order); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.QR(); + + var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(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.GenerateRandomUserDefinedMatrix(order, order); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.QR(); + + var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(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.GenerateRandomUserDefinedMatrix(order, order); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.QR(); + var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(order); + var vectorbCopy = vectorb.Clone(); + var resultx = new UserDefinedVector(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.GenerateRandomUserDefinedMatrix(order, order); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.QR(); + + var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(order, order); + var matrixBCopy = matrixB.Clone(); + + var matrixX = new UserDefinedMatrix(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]); + } + } + } + } +} diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserSvdTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserSvdTests.cs new file mode 100644 index 00000000..607ca312 --- /dev/null +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/UserSvdTests.cs @@ -0,0 +1,369 @@ +// +// 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. +// + +namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization +{ + using System; + using MbUnit.Framework; + using LinearAlgebra.Double.Factorization; + + public class UserSvdTests + { + + [Test] + [ExpectedArgumentNullException] + public void ConstructorNull() + { + new UserSvd(null, true); + } + + [Test] + [Row(1)] + [Row(10)] + [Row(100)] + public void CanFactorizeIdentity(int order) + { + var I = UserDefinedMatrix.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.GenerateRandomUserDefinedMatrix(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.GenerateRandomUserDefinedMatrix(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.GenerateRandomUserDefinedMatrix(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 UserDefinedMatrix(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.GenerateRandomUserDefinedMatrix(10, 10); + var factorSvd = matrixA.Svd(false); + + var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(10, 10); + factorSvd.Solve(matrixB); + } + + [Test] + [ExpectedException(typeof(InvalidOperationException))] + public void CannotSolveVectorIfVectorsNotComputed() + { + var matrixA = MatrixLoader.GenerateRandomUserDefinedMatrix(10, 10); + var factorSvd = matrixA.Svd(false); + + var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(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.GenerateRandomUserDefinedMatrix(row, column); + var matrixACopy = matrixA.Clone(); + var factorSvd = matrixA.Svd(true); + + var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(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.GenerateRandomUserDefinedMatrix(row, count); + var matrixACopy = matrixA.Clone(); + var factorSvd = matrixA.Svd(true); + + var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(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.GenerateRandomUserDefinedMatrix(row, column); + var matrixACopy = matrixA.Clone(); + var factorSvd = matrixA.Svd(true); + var vectorb = MatrixLoader.GenerateRandomUserDefinedVector(row); + var vectorbCopy = vectorb.Clone(); + var resultx = new UserDefinedVector(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.GenerateRandomUserDefinedMatrix(row, column); + var matrixACopy = matrixA.Clone(); + var factorSvd = matrixA.Svd(true); + + var matrixB = MatrixLoader.GenerateRandomUserDefinedMatrix(row, column); + var matrixBCopy = matrixB.Clone(); + + var matrixX = new UserDefinedMatrix(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]); + } + } + } + } +} diff --git a/src/UnitTests/LinearAlgebraTests/Double/MatrixLoader.cs b/src/UnitTests/LinearAlgebraTests/Double/MatrixLoader.cs index 56884526..776bf053 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/MatrixLoader.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/MatrixLoader.cs @@ -63,7 +63,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double } } - public static Matrix GenerateRandomMatrix(int row, int col) + public static Matrix GenerateRandomDenseMatrix(int row, int col) { // Fill a matrix with standard random numbers. var normal = new Distributions.Normal(); @@ -81,7 +81,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double return A; } - public static Matrix GenerateRandomPositiveDefiniteMatrix(int order) + public static Matrix GenerateRandomPositiveDefiniteDenseMatrix(int order) { // Fill a matrix with standard random numbers. var normal = new Distributions.Normal(); @@ -99,7 +99,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double return A.Transpose() * A; } - public static Vector GenerateRandomVector(int order) + public static Vector GenerateRandomDenseVector(int order) { // Fill a matrix with standard random numbers. var normal = new Distributions.Normal(); @@ -113,5 +113,56 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double // Generate a matrix which is positive definite. return v; } + + public static Matrix GenerateRandomUserDefinedMatrix(int row, int col) + { + // Fill a matrix with standard random numbers. + var normal = new Distributions.Normal(); + normal.RandomSource = new Random.MersenneTwister(1); + var A = new UserDefinedMatrix(row, col); + for (int i = 0; i < row; i++) + { + for (int j = 0; j < col; j++) + { + A[i, j] = normal.Sample(); + } + } + + // Generate a matrix which is positive definite. + return A; + } + + public static Matrix GenerateRandomPositiveDefiniteUserDefinedMatrix(int order) + { + // Fill a matrix with standard random numbers. + var normal = new Distributions.Normal(); + normal.RandomSource = new Random.MersenneTwister(1); + var A = new UserDefinedMatrix(order); + for (int i = 0; i < order; i++) + { + for (int j = 0; j < order; j++) + { + A[i, j] = normal.Sample(); + } + } + + // Generate a matrix which is positive definite. + return A.Transpose() * A; + } + + public static Vector GenerateRandomUserDefinedVector(int order) + { + // Fill a matrix with standard random numbers. + var normal = new Distributions.Normal(); + normal.RandomSource = new Random.MersenneTwister(1); + var v = new UserDefinedVector(order); + for (int i = 0; i < order; i++) + { + v[i] = normal.Sample(); + } + + // Generate a matrix which is positive definite. + return v; + } } } diff --git a/src/UnitTests/LinearAlgebraTests/Double/UserDefinedMatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Double/UserDefinedMatrixTests.cs index b7ecdae3..0a3b98ea 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/UserDefinedMatrixTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/UserDefinedMatrixTests.cs @@ -36,6 +36,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double { private readonly double[,] _data; + public UserDefinedMatrix(int order): base(order, order) + { + _data = new double[order, order]; + } + public UserDefinedMatrix(int rows, int columns) : base(rows, columns) { _data = new double[rows, columns]; @@ -65,6 +70,17 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double { return new UserDefinedVector(size); } + + public static UserDefinedMatrix Identity(int order) + { + var m = new UserDefinedMatrix(order, order); + for (var i = 0; i < order; i++) + { + m[i, i] = 1.0; + } + + return m; + } } public class UserDefinedMatrixTests : MatrixTests diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 528eedac..7b2e5698 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -90,6 +90,10 @@ + + + +