From 53679353bf791efd4b63564379d6aad0efcedba5 Mon Sep 17 00:00:00 2001 From: Marcus Cuda Date: Wed, 13 Feb 2013 14:52:55 +0200 Subject: [PATCH] fixed dimension check with tall matrices and thin qr --- .../Complex/Factorization/DenseQR.cs | 8 +- .../Complex32/Factorization/DenseQR.cs | 8 +- .../Double/Factorization/DenseQR.cs | 8 +- .../Single/Factorization/DenseQR.cs | 8 +- .../Complex/Factorization/QRTests.cs | 76 +++++++++++++++++++ .../Complex32/Factorization/QRTests.cs | 75 ++++++++++++++++++ .../Double/Factorization/QRTests.cs | 76 +++++++++++++++++++ .../Single/Factorization/QRTests.cs | 76 +++++++++++++++++++ 8 files changed, 319 insertions(+), 16 deletions(-) diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs index 512b631a..e88d62ee 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs @@ -121,7 +121,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization } // 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) + if (MatrixQ.RowCount != input.RowCount) { throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); } @@ -144,7 +144,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixR.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, input.ColumnCount, dresult.Values, QrMethod); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, input.ColumnCount, dresult.Values, QrMethod); } /// @@ -166,7 +166,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization // 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) + if (MatrixQ.RowCount != input.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); } @@ -189,7 +189,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixR.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, 1, dresult.Values, QrMethod); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, 1, dresult.Values, QrMethod); } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs index 78a2dd95..1eb56b12 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs @@ -121,7 +121,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization } // 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) + if (MatrixQ.RowCount != input.RowCount) { throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); } @@ -144,7 +144,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixR.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, input.ColumnCount, dresult.Values, QrMethod); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, input.ColumnCount, dresult.Values, QrMethod); } /// @@ -166,7 +166,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization // 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) + if (MatrixQ.RowCount != input.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); } @@ -189,7 +189,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixR.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, 1, dresult.Values, QrMethod); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, 1, dresult.Values, QrMethod); } } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs index 9294d033..c92b1ac5 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs @@ -121,7 +121,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization } // 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) + if (MatrixQ.RowCount != input.RowCount) { throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); } @@ -144,7 +144,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixR.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, input.ColumnCount, dresult.Values, QrMethod); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, input.ColumnCount, dresult.Values, QrMethod); } /// @@ -166,7 +166,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization // 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) + if (MatrixQ.RowCount != input.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); } @@ -189,7 +189,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixR.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, 1, dresult.Values, QrMethod); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, 1, dresult.Values, QrMethod); } } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs index 2fb6d3a3..474e1e5a 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs @@ -120,7 +120,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization } // 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) + if (MatrixQ.RowCount != input.RowCount) { throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); } @@ -143,7 +143,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixR.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, input.ColumnCount, dresult.Values, QrMethod); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, input.ColumnCount, dresult.Values, QrMethod); } /// @@ -165,7 +165,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization // 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) + if (MatrixQ.RowCount != input.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); } @@ -188,7 +188,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixR.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, 1, dresult.Values, QrMethod); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, 1, dresult.Values, QrMethod); } } } diff --git a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/QRTests.cs b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/QRTests.cs index 00b682a9..505f1dce 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/QRTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/QRTests.cs @@ -647,5 +647,81 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization } } } + + /// + /// Can solve when using a tall matrix. + /// + /// The QR decomp method to use. + [TestCase(QRMethod.Full)] + [TestCase(QRMethod.Thin)] + public void CanSolveForMatrixWithTallRandomMatrix(QRMethod method) + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(20, 10); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.QR(method); + + var matrixB = MatrixLoader.GenerateRandomDenseMatrix(20, 5); + 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 test = (matrixA.ConjugateTranspose() * matrixA).Inverse() * matrixA.ConjugateTranspose() * matrixB; + + for (var i = 0; i < matrixX.RowCount; i++) + { + for (var j = 0; j < matrixX.ColumnCount; j++) + { + AssertHelpers.AlmostEqual(test[i, j], matrixX[i, j], 9); + } + } + + // 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]); + } + } + } + + /// + /// Can solve when using a tall matrix. + /// + /// The QR decomp method to use. + [TestCase(QRMethod.Full)] + [TestCase(QRMethod.Thin)] + public void CanSolveForVectorWithTallRandomMatrix(QRMethod method) + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(20, 10); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.QR(method); + + var vectorB = MatrixLoader.GenerateRandomDenseVector(20); + var vectorX = factorQR.Solve(vectorB); + + // The solution x dimension is equal to the column dimension of A + Assert.AreEqual(matrixA.ColumnCount, vectorX.Count); + + var test = (matrixA.ConjugateTranspose() * matrixA).Inverse() * matrixA.ConjugateTranspose()* vectorB; + + for (var i = 0; i < vectorX.Count; i++) + { + AssertHelpers.AlmostEqual(test[i], vectorX[i], 9); + } + + // 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]); + } + } + } } } diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/QRTests.cs b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/QRTests.cs index 04987dff..5387d77d 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/QRTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/QRTests.cs @@ -656,5 +656,80 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization } } } + /// + /// Can solve when using a tall matrix. + /// + /// The QR decomp method to use. + [TestCase(QRMethod.Full)] + [TestCase(QRMethod.Thin)] + public void CanSolveForMatrixWithTallRandomMatrix(QRMethod method) + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(20, 10); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.QR(method); + + var matrixB = MatrixLoader.GenerateRandomDenseMatrix(20, 5); + 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 test = (matrixA.ConjugateTranspose() * matrixA).Inverse() * matrixA.ConjugateTranspose() * matrixB; + + for (var i = 0; i < matrixX.RowCount; i++) + { + for (var j = 0; j < matrixX.ColumnCount; j++) + { + AssertHelpers.AlmostEqual(test[i, j], matrixX[i, j], 4); + } + } + + // 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]); + } + } + } + + /// + /// Can solve when using a tall matrix. + /// + /// The QR decomp method to use. + [TestCase(QRMethod.Full)] + [TestCase(QRMethod.Thin)] + public void CanSolveForVectorWithTallRandomMatrix(QRMethod method) + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(20, 10); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.QR(method); + + var vectorB = MatrixLoader.GenerateRandomDenseVector(20); + var vectorX = factorQR.Solve(vectorB); + + // The solution x dimension is equal to the column dimension of A + Assert.AreEqual(matrixA.ColumnCount, vectorX.Count); + + var test = (matrixA.ConjugateTranspose() * matrixA).Inverse() * matrixA.ConjugateTranspose() * vectorB; + + for (var i = 0; i < vectorX.Count; i++) + { + AssertHelpers.AlmostEqual(test[i], vectorX[i], 4); + } + + // 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]); + } + } + } } } diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/QRTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/QRTests.cs index f75c7fe3..cb44e96c 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Factorization/QRTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/QRTests.cs @@ -642,5 +642,81 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization } } } + + /// + /// Can solve when using a tall matrix. + /// + /// The QR decomp method to use. + [TestCase(QRMethod.Full)] + [TestCase(QRMethod.Thin)] + public void CanSolveForMatrixWithTallRandomMatrix(QRMethod method) + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(20, 10); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.QR(method); + + var matrixB = MatrixLoader.GenerateRandomDenseMatrix(20, 5); + 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 test = (matrixA.Transpose() * matrixA).Inverse() * matrixA.Transpose() * matrixB; + + for (var i = 0; i < matrixX.RowCount; i++) + { + for (var j = 0; j < matrixX.ColumnCount; j++) + { + AssertHelpers.AlmostEqual(test[i, j], matrixX[i, j], 9); + } + } + + // 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]); + } + } + } + + /// + /// Can solve when using a tall matrix. + /// + /// The QR decomp method to use. + [TestCase(QRMethod.Full)] + [TestCase(QRMethod.Thin)] + public void CanSolveForVectorWithTallRandomMatrix(QRMethod method) + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(20, 10); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.QR(method); + + var vectorB = MatrixLoader.GenerateRandomDenseVector(20); + var vectorX = factorQR.Solve(vectorB); + + // The solution x dimension is equal to the column dimension of A + Assert.AreEqual(matrixA.ColumnCount, vectorX.Count); + + var test = (matrixA.Transpose() * matrixA).Inverse() * matrixA.Transpose() * vectorB; + + for (var i = 0; i < vectorX.Count; i++) + { + AssertHelpers.AlmostEqual(test[i], vectorX[i], 9); + } + + // 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]); + } + } + } } } diff --git a/src/UnitTests/LinearAlgebraTests/Single/Factorization/QRTests.cs b/src/UnitTests/LinearAlgebraTests/Single/Factorization/QRTests.cs index 2ad97326..67dd3448 100644 --- a/src/UnitTests/LinearAlgebraTests/Single/Factorization/QRTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Single/Factorization/QRTests.cs @@ -643,5 +643,81 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization } } } + + /// + /// Can solve when using a tall matrix. + /// + /// The QR decomp method to use. + [TestCase(QRMethod.Full)] + [TestCase(QRMethod.Thin)] + public void CanSolveForMatrixWithTallRandomMatrix(QRMethod method) + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(20, 10); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.QR(method); + + var matrixB = MatrixLoader.GenerateRandomDenseMatrix(20, 5); + 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 test = (matrixA.Transpose() * matrixA).Inverse() * matrixA.Transpose() * matrixB; + + for (var i = 0; i < matrixX.RowCount; i++) + { + for (var j = 0; j < matrixX.ColumnCount; j++) + { + AssertHelpers.AlmostEqual(test[i, j], matrixX[i, j], 4); + } + } + + // 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]); + } + } + } + + /// + /// Can solve when using a tall matrix. + /// + /// The QR decomp method to use. + [TestCase(QRMethod.Full)] + [TestCase(QRMethod.Thin)] + public void CanSolveForVectorWithTallRandomMatrix(QRMethod method) + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(20, 10); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.QR(method); + + var vectorB = MatrixLoader.GenerateRandomDenseVector(20); + var vectorX = factorQR.Solve(vectorB); + + // The solution x dimension is equal to the column dimension of A + Assert.AreEqual(matrixA.ColumnCount, vectorX.Count); + + var test = (matrixA.Transpose() * matrixA).Inverse() * matrixA.Transpose() * vectorB; + + for (var i = 0; i < vectorX.Count; i++) + { + AssertHelpers.AlmostEqual(test[i], vectorX[i], 4); + } + + // 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]); + } + } + } } }