search

mathnet / mathnet-numerics

Math.NET Numerics
http://numerics.mathdotnet.com
MIT License
3.48k stars 894 forks source link

QRFactor error in native providers #490

Open armwal opened 7 years ago

armwal commented 7 years ago

Hi, I noticed a discrepancy between the managed provider and the native providers in computing the QR factorization of a small test matrix.

According to Matlab, the 2x4 test matrix

testData =

  -2.185929268096052  -0.458914070137267                   0  -0.748405712906516
  -2.262514008467083   0.443380104476701                   0  -1.064633478923354

should be factored as:

>> [Q,R] = qr(testData,0)

Q =

  -0.694830212770919  -0.719173814470967
  -0.719173814470967   0.694830212770918

R =

   3.145990528216626   0.000000000000000                   0   1.285671420788610
                   0   0.638112874666953                   0  -0.201505715360515

If I reproduce this QR factorization in Math.NET with the Control.LinearAlgebraProvider.QRFactor method, the managed linear algebra provider produces essentially the same results as Matlab (except for sign changes), while both the MKL and the OpenBLAS native providers do not: For both of them, the entry R[1,0] equals about 0.42 instead of 0. This means, that for the native providers Q*R does not equal the test matrix, while this test works for the results of the managed provider.

I've attached a test class (sorry for the .txt format, GitHub wouldn't accept .cs) reproducing the problem. QRFactorizationTest.txt

I'm using the NuGet packages MathNet.Numerics V 3.19.0 MathNet.Numerics.MKL.Win-x64 V 2.2.0 MathNet.Numerics.OpenBLAS.Win V 0.2.0

The matrices I'm working on are all rectangular m x n with m < n and I'm aware that other entry points for the QR factorization in Math.NET such as the Double.Factorization.UserQR restrict the QR factorization to matrices with m >= n. If this restriction is in place because of problems like the one I noticed, do you by chance have a suggestion on how to achieve the QR factorization for matrices with m < n in Math.NET?

Kind regards Armin

cdrnet commented 7 years ago

I assume we do not support m < n in this case since the QR decomposition was mainly provided as a means to solve linear least squares problems, but with an under-determined system like this we no longer have a least squares problem. But I do see that the decomposition itself can still be useful in practice.

I just have a look at the way we call the LAPACK routines within the native providers, and they do indeed seem wrong for m < n. This code is shared between MKL and OpenBLAS, so the same (wrong) behavior is expected there. We should be able to fix this with a new provider release. Maybe we can then also drop the restriction entirely.

Thanks for reporting this, and for providing a repro!

tournierjc commented 4 years ago

Hello,

Is this issue still present upstream?

armwal commented 4 years ago

Hi,

I did a short test, and the most recent versions of MathNet.Numerics (V 4.12.0), MathNet.Numerics.MKL.Win.x64 (V 2.4.0) and MathNet.Numerics.OpenBLAS.Win (V 0.2.0) still have the same issue.