diff --git a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs index a06f1af9..a6005e41 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs @@ -412,8 +412,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be /// null for the native provider. The native provider uses the Q portion stored in the R matrix. /// The R matrix obtained by calling . - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver /// and can be null for the managed provider. /// On entry the B matrix; on exit the X matrix. @@ -421,7 +421,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// On exit, the solution matrix. /// Rows must be greater or equal to columns. /// The type of QR factorization to perform. - void QRSolveFactored(T[] q, T[] r, int rowsR, int columnsR, T[] tau, T[] b, int columnsB, T[] x, QRMethod method = QRMethod.Full); + void QRSolveFactored(T[] q, T[] r, int rowsA, int columnsA, T[] tau, T[] b, int columnsB, T[] x, QRMethod method = QRMethod.Full); /// /// Solves A*X=B for X using a previously QR factored matrix. @@ -429,8 +429,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be /// null for the native provider. The native provider uses the Q portion stored in the R matrix. /// The R matrix obtained by calling . - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver /// and can be null for the managed provider. /// On entry the B matrix; on exit the X matrix. @@ -441,7 +441,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// work size value. /// Rows must be greater or equal to columns. /// The type of QR factorization to perform. - void QRSolveFactored(T[] q, T[] r, int rowsR, int columnsR, T[] tau, T[] b, int columnsB, T[] x, T[] work, QRMethod method = QRMethod.Full); + void QRSolveFactored(T[] q, T[] r, int rowsA, int columnsA, T[] tau, T[] b, int columnsB, T[] x, T[] work, QRMethod method = QRMethod.Full); /// /// Computes the singular value decomposition of A. diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseGramSchmidt.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseGramSchmidt.cs index df7aa4b5..318a3ea3 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseGramSchmidt.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseGramSchmidt.cs @@ -28,6 +28,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using MathNet.Numerics.LinearAlgebra.Generic.Factorization; + namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization { using System; @@ -172,7 +174,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment."); } - _provider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixQ.ColumnCount, null, dinput.Values, input.ColumnCount, dresult.Values); + _provider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixR.ColumnCount, null, dinput.Values, input.ColumnCount, dresult.Values, QRMethod.Thin); } /// @@ -217,7 +219,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment."); } - _provider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixQ.ColumnCount, null, dinput.Values, 1, dresult.Values); + _provider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixR.ColumnCount, null, dinput.Values, 1, dresult.Values, QRMethod.Thin); } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseGramSchmidt.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseGramSchmidt.cs index 2ee13acd..258f037d 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseGramSchmidt.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseGramSchmidt.cs @@ -28,6 +28,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using MathNet.Numerics.LinearAlgebra.Generic.Factorization; + namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { using System; @@ -172,7 +174,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment."); } - _provider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixQ.ColumnCount, null, dinput.Values, input.ColumnCount, dresult.Values); + _provider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixR.ColumnCount, null, dinput.Values, input.ColumnCount, dresult.Values, QRMethod.Thin); } /// @@ -217,7 +219,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment."); } - _provider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixQ.ColumnCount, null, dinput.Values, 1, dresult.Values); + _provider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixR.ColumnCount, null, dinput.Values, 1, dresult.Values, QRMethod.Thin); } } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseGramSchmidt.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseGramSchmidt.cs index 80ff1d93..f1b04b82 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseGramSchmidt.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseGramSchmidt.cs @@ -28,6 +28,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using MathNet.Numerics.LinearAlgebra.Generic.Factorization; + namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { using System; @@ -171,7 +173,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment."); } - _provider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixR.RowCount, MatrixR.ColumnCount, null, dinput.Values, input.ColumnCount, dresult.Values); + _provider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixR.ColumnCount, null, dinput.Values, input.ColumnCount, dresult.Values, QRMethod.Thin); } /// @@ -216,7 +218,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment."); } - _provider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixQ.ColumnCount, null, dinput.Values, 1, dresult.Values); + _provider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixR.ColumnCount, null, dinput.Values, 1, dresult.Values, QRMethod.Thin); } } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseGramSchmidt.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseGramSchmidt.cs index d0f82cd0..04ce0aad 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseGramSchmidt.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseGramSchmidt.cs @@ -28,6 +28,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using MathNet.Numerics.LinearAlgebra.Generic.Factorization; + namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { using System; @@ -171,7 +173,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment."); } - _provider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixQ.ColumnCount, null, dinput.Values, input.ColumnCount, dresult.Values); + _provider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixR.ColumnCount, null, dinput.Values, input.ColumnCount, dresult.Values, QRMethod.Thin); } /// @@ -216,7 +218,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment."); } - _provider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixQ.ColumnCount, null, dinput.Values, 1, dresult.Values); + _provider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixQ.RowCount, MatrixR.ColumnCount, null, dinput.Values, 1, dresult.Values, QRMethod.Thin); } } } diff --git a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/GramSchmidtTests.cs b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/GramSchmidtTests.cs index 3337c886..abad3099 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/GramSchmidtTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/GramSchmidtTests.cs @@ -372,5 +372,77 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization } } } + + /// + /// Can solve when using a tall matrix. + /// + [Test] + public void CanSolveForMatrixWithTallRandomMatrix() + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(20, 10); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.GramSchmidt(); + + 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. + /// + [Test] + public void CanSolveForVectorWithTallRandomMatrix() + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(20, 10); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.GramSchmidt(); + + 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/GramSchmidtTests.cs b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/GramSchmidtTests.cs index fca0d180..d48c4122 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/GramSchmidtTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/GramSchmidtTests.cs @@ -379,5 +379,77 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization } } } + + /// + /// Can solve when using a tall matrix. + /// + [Test] + public void CanSolveForMatrixWithTallRandomMatrix() + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(20, 10); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.GramSchmidt(); + + 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. + /// + [Test] + public void CanSolveForVectorWithTallRandomMatrix() + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(20, 10); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.GramSchmidt(); + + 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/GramSchmidtTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/GramSchmidtTests.cs index 657cc48d..f66d9e20 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Factorization/GramSchmidtTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/GramSchmidtTests.cs @@ -354,5 +354,77 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization } } } + + /// + /// Can solve when using a tall matrix. + /// + [Test] + public void CanSolveForMatrixWithTallRandomMatrix() + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(20, 10); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.GramSchmidt(); + + 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. + /// + [Test] + public void CanSolveForVectorWithTallRandomMatrix() + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(20, 10); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.GramSchmidt(); + + 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/GramSchmidtTests.cs b/src/UnitTests/LinearAlgebraTests/Single/Factorization/GramSchmidtTests.cs index 5f299fda..a114940d 100644 --- a/src/UnitTests/LinearAlgebraTests/Single/Factorization/GramSchmidtTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Single/Factorization/GramSchmidtTests.cs @@ -354,5 +354,77 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization } } } + + /// + /// Can solve when using a tall matrix. + /// + [Test] + public void CanSolveForMatrixWithTallRandomMatrix() + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(20, 10); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.GramSchmidt(); + + 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. + /// + [Test] + public void CanSolveForVectorWithTallRandomMatrix() + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(20, 10); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.GramSchmidt(); + + 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]); + } + } + } } }