fixed Gram Schmidt bug for solving tall matrices

src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs

@@ -412,8 +412,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="q">The Q matrix obtained by QR factor. This is only used for the managed provider and can be
         /// <c>null</c> for the native provider. The native provider uses the Q portion stored in the R matrix.</param>
         /// <param name="r">The R matrix obtained by calling <see cref="QRFactor(T[],int,int,T[],T[])"/>. </param>
-        /// <param name="rowsR">The number of rows in the A matrix.</param>
-        /// <param name="columnsR">The number of columns in the A matrix.</param>
+        /// <param name="rowsA">The number of rows in the A matrix.</param>
+        /// <param name="columnsA">The number of columns in the A matrix.</param>
         /// <param name="tau">Contains additional information on Q. Only used for the native solver
         /// and can be <c>null</c> for the managed provider.</param>
         /// <param name="b">On entry the B matrix; on exit the X matrix.</param>
@@ -421,7 +421,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="x">On exit, the solution matrix.</param>
         /// <remarks>Rows must be greater or equal to columns.</remarks>
         /// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
-        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);
 
         /// <summary>
         /// Solves A*X=B for X using a previously QR factored matrix.
@@ -429,8 +429,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="q">The Q matrix obtained by QR factor. This is only used for the managed provider and can be
         /// <c>null</c> for the native provider. The native provider uses the Q portion stored in the R matrix.</param>
         /// <param name="r">The R matrix obtained by calling <see cref="QRFactor(T[],int,int,T[],T[])"/>. </param>
-        /// <param name="rowsR">The number of rows in the A matrix.</param>
-        /// <param name="columnsR">The number of columns in the A matrix.</param>
+        /// <param name="rowsA">The number of rows in the A matrix.</param>
+        /// <param name="columnsA">The number of columns in the A matrix.</param>
         /// <param name="tau">Contains additional information on Q. Only used for the native solver
         /// and can be <c>null</c> for the managed provider.</param>
         /// <param name="b">On entry the B matrix; on exit the X matrix.</param>
@@ -441,7 +441,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// work size value.</param>
         /// <remarks>Rows must be greater or equal to columns.</remarks>
         /// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
-        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);
         
         /// <summary>
         /// Computes the singular value decomposition of A.

src/Numerics/LinearAlgebra/Complex/Factorization/DenseGramSchmidt.cs

@@ -28,6 +28,8 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
+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);
         }
 
         /// <summary>
@@ -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);
         }
     }
 }

src/Numerics/LinearAlgebra/Complex32/Factorization/DenseGramSchmidt.cs

@@ -28,6 +28,8 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
+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);
         }
 
         /// <summary>
@@ -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);
         }
     }
 }

src/Numerics/LinearAlgebra/Double/Factorization/DenseGramSchmidt.cs

@@ -28,6 +28,8 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
+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);
         }
 
         /// <summary>
@@ -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);
         }
     }
 }

src/Numerics/LinearAlgebra/Single/Factorization/DenseGramSchmidt.cs

@@ -28,6 +28,8 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
+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);
         }
 
         /// <summary>
@@ -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);
         }
     }
 }

src/UnitTests/LinearAlgebraTests/Complex/Factorization/GramSchmidtTests.cs

@@ -372,5 +372,77 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization
                 }
             }
         }
+
+        /// <summary>
+        /// Can solve when using a tall matrix.
+        /// </summary>
+        [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]);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Can solve when using a tall matrix.
+        /// </summary>
+        [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]);
+                }
+            }
+        }
     }
 }

src/UnitTests/LinearAlgebraTests/Complex32/Factorization/GramSchmidtTests.cs

@@ -379,5 +379,77 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization
                 }
             }
         }
+
+        /// <summary>
+        /// Can solve when using a tall matrix.
+        /// </summary>
+        [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]);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Can solve when using a tall matrix.
+        /// </summary>
+        [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]);
+                }
+            }
+        }
     }
 }

src/UnitTests/LinearAlgebraTests/Double/Factorization/GramSchmidtTests.cs

@@ -354,5 +354,77 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
                 }
             }
         }
+
+        /// <summary>
+        /// Can solve when using a tall matrix.
+        /// </summary>
+        [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]);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Can solve when using a tall matrix.
+        /// </summary>
+        [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]);
+                }
+            }
+        }
     }
 }

src/UnitTests/LinearAlgebraTests/Single/Factorization/GramSchmidtTests.cs

@@ -354,5 +354,77 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization
                 }
             }
         }
+
+        /// <summary>
+        /// Can solve when using a tall matrix.
+        /// </summary>
+        [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]);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Can solve when using a tall matrix.
+        /// </summary>
+        [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]);
+                }
+            }
+        }
     }
 }