From b22b5a20aa95565920b96ad742349d07e0980a64 Mon Sep 17 00:00:00 2001 From: Marcus Cuda Date: Fri, 14 Dec 2012 12:13:10 +0200 Subject: [PATCH] fixed QR thin solve bug native QR full now returns the full Q matrix --- src/NativeWrappers/MKL/lapack.cpp | 2 +- .../Mkl/MklLinearAlgebraProvider.Common.cs | 3 +- .../Complex/Factorization/DenseQR.cs | 5 ++- .../Complex/Factorization/UserQR.cs | 1 + .../Complex32/Factorization/DenseQR.cs | 5 ++- .../Complex32/Factorization/UserQR.cs | 1 + .../Double/Factorization/DenseQR.cs | 5 ++- .../Double/Factorization/UserQR.cs | 1 + .../LinearAlgebra/Generic/Factorization/QR.cs | 10 +++++ .../Single/Factorization/DenseQR.cs | 5 ++- .../Single/Factorization/UserQR.cs | 1 + .../Complex/LinearAlgebraProviderTests.cs | 1 - .../Complex/Factorization/QRTests.cs | 37 +++++++++++++++++-- .../Complex/MatrixTests.Arithmetic.cs | 8 +++- .../Complex32/Factorization/QRTests.cs | 11 ++++++ .../Complex32/MatrixTests.cs | 1 - .../Double/Factorization/QRTests.cs | 34 +++++++++++++++++ .../LinearAlgebraTests/Double/MatrixTests.cs | 14 +++---- .../Single/Factorization/QRTests.cs | 35 ++++++++++++++++++ 19 files changed, 155 insertions(+), 25 deletions(-) diff --git a/src/NativeWrappers/MKL/lapack.cpp b/src/NativeWrappers/MKL/lapack.cpp index 09bd04ce..e3610b86 100644 --- a/src/NativeWrappers/MKL/lapack.cpp +++ b/src/NativeWrappers/MKL/lapack.cpp @@ -180,7 +180,7 @@ inline MKL_INT qr_factor(MKL_INT m, MKL_INT n, T r[], T tau[], T q[], T work[], } else { - orgqr(&m, &n, &n, q, &m, tau, work, &len, &info); + orgqr(&m, &m, &n, q, &m, tau, work, &len, &info); } return info; diff --git a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Common.cs b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Common.cs index d495ca05..1be99a38 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Common.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Common.cs @@ -40,7 +40,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// public partial class MklLinearAlgebraProvider : ManagedLinearAlgebraProvider { - /* /// + /// /// Computes the requested of the matrix. /// /// The type of norm to compute. @@ -200,6 +200,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl return SafeNativeMethods.d_matrix_norm((byte)norm, rows, columns, matrix, work); } + /* BUG in MKL'S ZLANGE routine. Using managed code until it is fixed. /// /// Computes the requested of the matrix. /// diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs index 95e71f8b..512b631a 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs @@ -77,6 +77,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw Matrix.DimensionsDontMatch(matrix); } + QrMethod = method; Tau = new Complex[Math.Min(matrix.RowCount, matrix.ColumnCount)]; if (method == QRMethod.Full) @@ -143,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); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixR.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, input.ColumnCount, dresult.Values, QrMethod); } /// @@ -188,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); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixR.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, 1, dresult.Values, QrMethod); } } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/UserQR.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/UserQR.cs index 04a1fcbe..2eba8cd0 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/UserQR.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/UserQR.cs @@ -69,6 +69,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw Matrix.DimensionsDontMatch(matrix); } + QrMethod = method; var minmn = Math.Min(matrix.RowCount, matrix.ColumnCount); var u = new Complex[minmn][]; diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs index 02741957..78a2dd95 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs @@ -77,6 +77,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw Matrix.DimensionsDontMatch(matrix); } + QrMethod = method; Tau = new Complex32[Math.Min(matrix.RowCount, matrix.ColumnCount)]; if (method == QRMethod.Full) @@ -143,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); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixR.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, input.ColumnCount, dresult.Values, QrMethod); } /// @@ -188,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); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixR.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, 1, dresult.Values, QrMethod); } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserQR.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserQR.cs index 4cc9907a..d2677ebc 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserQR.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserQR.cs @@ -69,6 +69,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw Matrix.DimensionsDontMatch(matrix); } + QrMethod = method; var minmn = Math.Min(matrix.RowCount, matrix.ColumnCount); var u = new Complex32[minmn][]; diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs index 3be73a16..9294d033 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs @@ -76,6 +76,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw Matrix.DimensionsDontMatch(matrix); } + QrMethod = method; Tau = new double[Math.Min(matrix.RowCount, matrix.ColumnCount)]; if (method == QRMethod.Full) @@ -143,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); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixR.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, input.ColumnCount, dresult.Values, QrMethod); } /// @@ -188,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); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixR.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, 1, dresult.Values, QrMethod); } } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs b/src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs index eaac74ee..b2246ec4 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs @@ -68,6 +68,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw Matrix.DimensionsDontMatch(matrix); } + QrMethod = method; var minmn = Math.Min(matrix.RowCount, matrix.ColumnCount); var u = new double[minmn][]; diff --git a/src/Numerics/LinearAlgebra/Generic/Factorization/QR.cs b/src/Numerics/LinearAlgebra/Generic/Factorization/QR.cs index 0afbe870..9d4dee24 100644 --- a/src/Numerics/LinearAlgebra/Generic/Factorization/QR.cs +++ b/src/Numerics/LinearAlgebra/Generic/Factorization/QR.cs @@ -81,6 +81,15 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization set; } + /// + /// The QR factorization method. + /// + protected QRMethod QrMethod + { + get; + set; + } + /// /// Internal method which routes the call to perform the QR factorization to the appropriate class. /// @@ -89,6 +98,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic.Factorization /// A QR factorization object. internal static QR Create(Matrix matrix, QRMethod method = QRMethod.Full) { + if (typeof(T) == typeof(double)) { var dense = matrix as LinearAlgebra.Double.DenseMatrix; diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs index b8701d61..2fb6d3a3 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs @@ -76,6 +76,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw Matrix.DimensionsDontMatch(matrix); } + QrMethod = method; Tau = new float[Math.Min(matrix.RowCount, matrix.ColumnCount)]; if (method == QRMethod.Full) @@ -142,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); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixR.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, input.ColumnCount, dresult.Values, QrMethod); } /// @@ -187,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); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Values, ((DenseMatrix)MatrixR).Values, MatrixR.RowCount, MatrixR.ColumnCount, Tau, dinput.Values, 1, dresult.Values, QrMethod); } } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/UserQR.cs b/src/Numerics/LinearAlgebra/Single/Factorization/UserQR.cs index 047c4ff7..6aad3dcb 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/UserQR.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/UserQR.cs @@ -68,6 +68,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw Matrix.DimensionsDontMatch(matrix); } + QrMethod = method; var minmn = Math.Min(matrix.RowCount, matrix.ColumnCount); var u = new float[minmn][]; diff --git a/src/UnitTests/LinearAlgebraProviderTests/Complex/LinearAlgebraProviderTests.cs b/src/UnitTests/LinearAlgebraProviderTests/Complex/LinearAlgebraProviderTests.cs index 08ce02ed..cdd11c21 100644 --- a/src/UnitTests/LinearAlgebraProviderTests/Complex/LinearAlgebraProviderTests.cs +++ b/src/UnitTests/LinearAlgebraProviderTests/Complex/LinearAlgebraProviderTests.cs @@ -919,7 +919,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex var mx = new DenseMatrix(matrix.ColumnCount, 2, x); var mb = matrix * mx; - Console.WriteLine(mx); AssertHelpers.AlmostEqual(mb[0, 0], b[0], 14); AssertHelpers.AlmostEqual(mb[1, 0], b[1], 14); AssertHelpers.AlmostEqual(mb[2, 0], b[2], 14); diff --git a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/QRTests.cs b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/QRTests.cs index 0be6f6e6..00b682a9 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex/Factorization/QRTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex/Factorization/QRTests.cs @@ -181,6 +181,25 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization AssertHelpers.AlmostEqual(matrixA[i, j], matrixQfromR[i, j], 9); } } + + // Make sure the Q is unitary --> (Q*)x(Q) = I + var matrixQсtQ = q.ConjugateTranspose() * q; + for (var i = 0; i < matrixQсtQ.RowCount; i++) + { + for (var j = 0; j < matrixQсtQ.ColumnCount; j++) + { + if (i == j) + { + Assert.AreEqual(matrixQсtQ[i, j].Real, 1.0, 1e-3); + Assert.AreEqual(matrixQсtQ[i, j].Imaginary, 0.0, 1e-3); + } + else + { + Assert.AreEqual(matrixQсtQ[i, j].Real, 0.0, 1e-3); + Assert.AreEqual(matrixQсtQ[i, j].Imaginary, 0.0, 1e-3); + } + } + } } /// @@ -221,6 +240,16 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization } } + // Make sure the Q*R is the original matrix. + var matrixQfromR = q * r; + for (var i = 0; i < matrixQfromR.RowCount; i++) + { + for (var j = 0; j < matrixQfromR.ColumnCount; j++) + { + AssertHelpers.AlmostEqual(matrixA[i, j], matrixQfromR[i, j], 9); + } + } + // Make sure the Q is unitary --> (Q*)x(Q) = I var matrixQсtQ = q.ConjugateTranspose() * q; for (var i = 0; i < matrixQсtQ.RowCount; i++) @@ -229,13 +258,13 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Factorization { if (i == j) { - Assert.AreEqual(matrixQсtQ[i, j].Real, 1.0f, 1e-3f); - Assert.AreEqual(matrixQсtQ[i, j].Imaginary, 0.0f, 1e-3f); + Assert.AreEqual(matrixQсtQ[i, j].Real, 1.0, 1e-3); + Assert.AreEqual(matrixQсtQ[i, j].Imaginary, 0.0, 1e-3); } else { - Assert.AreEqual(matrixQсtQ[i, j].Real, 0.0f, 1e-3f); - Assert.AreEqual(matrixQсtQ[i, j].Imaginary, 0.0f, 1e-3f); + Assert.AreEqual(matrixQсtQ[i, j].Real, 0.0, 1e-3); + Assert.AreEqual(matrixQсtQ[i, j].Imaginary, 0.0, 1e-3); } } } diff --git a/src/UnitTests/LinearAlgebraTests/Complex/MatrixTests.Arithmetic.cs b/src/UnitTests/LinearAlgebraTests/Complex/MatrixTests.Arithmetic.cs index f7aefa6e..e1117706 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex/MatrixTests.Arithmetic.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex/MatrixTests.Arithmetic.cs @@ -1067,7 +1067,9 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex { for (var j = 0; j < data.ColumnCount; j++) { - Assert.AreEqual(data[i, j] * other[i, j], result[i, j]); + var value = data[i, j]*other[i, j]; + Assert.AreEqual(value.Real, result[i, j].Real, 1e-12); + Assert.AreEqual(value.Imaginary, result[i, j].Imaginary, 1e-12); } } @@ -1076,7 +1078,9 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex { for (var j = 0; j < data.ColumnCount; j++) { - Assert.AreEqual(data[i, j] * other[i, j], result[i, j]); + var value = data[i, j] * other[i, j]; + Assert.AreEqual(value.Real, result[i, j].Real, 1e-12); + Assert.AreEqual(value.Imaginary, result[i, j].Imaginary, 1e-12); } } } diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/QRTests.cs b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/QRTests.cs index 9509d97e..c4df5c6c 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/QRTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex32/Factorization/QRTests.cs @@ -242,6 +242,17 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Factorization } } + // Make sure the Q*R is the original matrix. + var matrixQfromR = q * r; + for (var i = 0; i < matrixQfromR.RowCount; i++) + { + for (var j = 0; j < matrixQfromR.ColumnCount; j++) + { + Assert.AreEqual(matrixA[i, j].Real, matrixQfromR[i, j].Real, 1e-3f); + Assert.AreEqual(matrixA[i, j].Imaginary, matrixQfromR[i, j].Imaginary, 1e-3f); + } + } + // Make sure the Q is unitary --> (Q*)x(Q) = I var matrixQсtQ = q.ConjugateTranspose() * q; for (var i = 0; i < matrixQсtQ.RowCount; i++) diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/MatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Complex32/MatrixTests.cs index 35b46a5c..36c463c7 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex32/MatrixTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex32/MatrixTests.cs @@ -27,7 +27,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32 { using NUnit.Framework; - using Complex32 = Numerics.Complex32; /// /// Abstract class with the common set of matrix tests diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/QRTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/QRTests.cs index a6578b8d..f75c7fe3 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Factorization/QRTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/QRTests.cs @@ -180,6 +180,23 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization Assert.AreEqual(matrixA[i, j], matrixQfromR[i, j], 1.0e-11); } } + + // Make sure the Q is unitary --> (Q*)x(Q) = + var matrixQtQ = q.Transpose() * q; + for (var i = 0; i < matrixQtQ.RowCount; i++) + { + for (var j = 0; j < matrixQtQ.ColumnCount; j++) + { + if (i == j) + { + Assert.AreEqual(matrixQtQ[i, j], 1.0, 1e-3); + } + else + { + Assert.AreEqual(matrixQtQ[i, j], 0.0, 1e-3); + } + } + } } /// @@ -229,6 +246,23 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization Assert.AreEqual(matrixA[i, j], matrixQfromR[i, j], 1.0e-11); } } + + // Make sure the Q is unitary --> (Q*)x(Q) = + var matrixQtQ = q.Transpose() * q; + for (var i = 0; i < matrixQtQ.RowCount; i++) + { + for (var j = 0; j < matrixQtQ.ColumnCount; j++) + { + if (i == j) + { + Assert.AreEqual(matrixQtQ[i, j], 1.0, 1e-3); + } + else + { + Assert.AreEqual(matrixQtQ[i, j], 0.0, 1e-3); + } + } + } } /// diff --git a/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.cs index 388ac1f8..94edad1f 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.cs @@ -66,13 +66,13 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double public virtual void CanComputeFrobeniusNorm() { var matrix = TestMatrices["Square3x3"]; - AssertHelpers.AlmostEqual(10.77775486824598, matrix.FrobeniusNorm(), 14); + AssertHelpers.AlmostEqual(10.77775486824598, matrix.FrobeniusNorm(), 7); matrix = TestMatrices["Wide2x3"]; - AssertHelpers.AlmostEqual(4.79478883789474, matrix.FrobeniusNorm(), 14); + AssertHelpers.AlmostEqual(4.79478883789474, matrix.FrobeniusNorm(), 7); matrix = TestMatrices["Tall3x2"]; - AssertHelpers.AlmostEqual(7.54122006044115, matrix.FrobeniusNorm(), 14); + AssertHelpers.AlmostEqual(7.54122006044115, matrix.FrobeniusNorm(), 7); } /// @@ -85,10 +85,10 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double Assert.AreEqual(16.5, matrix.InfinityNorm()); matrix = TestMatrices["Wide2x3"]; - Assert.AreEqual(6.6, matrix.InfinityNorm()); + Assert.AreEqual(6.6, matrix.InfinityNorm(), 1e-7); matrix = TestMatrices["Tall3x2"]; - Assert.AreEqual(9.9, matrix.InfinityNorm()); + Assert.AreEqual(9.9, matrix.InfinityNorm(), 1e-4); } /// @@ -98,13 +98,13 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double public virtual void CanComputeL1Norm() { var matrix = TestMatrices["Square3x3"]; - Assert.AreEqual(12.1, matrix.L1Norm()); + Assert.AreEqual(12.1, matrix.L1Norm(), 1e-4); matrix = TestMatrices["Wide2x3"]; Assert.AreEqual(5.5, matrix.L1Norm()); matrix = TestMatrices["Tall3x2"]; - Assert.AreEqual(8.8, matrix.L1Norm()); + Assert.AreEqual(8.8, matrix.L1Norm(), 1e-4); } /// diff --git a/src/UnitTests/LinearAlgebraTests/Single/Factorization/QRTests.cs b/src/UnitTests/LinearAlgebraTests/Single/Factorization/QRTests.cs index 93d06e9f..ee910936 100644 --- a/src/UnitTests/LinearAlgebraTests/Single/Factorization/QRTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Single/Factorization/QRTests.cs @@ -181,6 +181,24 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization Assert.AreEqual(matrixA[i, j], matrixQfromR[i, j], 1e-4); } } + + // Make sure the Q is unitary --> (Q*)x(Q) = I + var matrixQtQ = q.Transpose() * q; + + for (var i = 0; i < matrixQtQ.RowCount; i++) + { + for (var j = 0; j < matrixQtQ.ColumnCount; j++) + { + if (i == j) + { + Assert.AreEqual(matrixQtQ[i, j], 1.0f, 1e-3f); + } + else + { + Assert.AreEqual(matrixQtQ[i, j], 0.0f, 1e-3f); + } + } + } } /// @@ -230,6 +248,23 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization Assert.AreEqual(matrixA[i, j], matrixQfromR[i, j], 1.0e-4); } } + + // Make sure the Q is unitary --> (Q*)x(Q) = I + var matrixQtQ = q.Transpose() * q; + for (var i = 0; i < matrixQtQ.RowCount; i++) + { + for (var j = 0; j < matrixQtQ.ColumnCount; j++) + { + if (i == j) + { + Assert.AreEqual(matrixQtQ[i, j], 1.0f, 1e-3f); + } + else + { + Assert.AreEqual(matrixQtQ[i, j], 0.0f, 1e-3f); + } + } + } } ///