culaDeviceSgesv accuracy, for M

[jrk1982]

culaDeviceSgesv's results are very accurate for square and skinny systems (M > N), but
for 'fat' systems (M < N), the results seem inaccurate... I am not able to get even a single digit of precision.

Could you please verify this? Again, am using random matrices, with N = 1000.

[john]

Hello, can you send me the relevant calling code so I can test on my own? gesv is only valid for square matrices (ie it doesn't even have a parameter "M"), although it can accept multiple right-hand-sides.

[jrk1982]

I apologise. I mean the function Sgels, which returns the minimum norm solution for M<N.

I call it as:
culaDeviceSgels('N',
mA,
nA,
nB,
(culaDeviceFloat*) d_in,
lda,
(culaDeviceFloat*) d_out,
ldb);

Regards
--Karthik

[john]

I'll need some more details to continue. What are your mA, nA, nB and which Matlab function are you comparing to?

[jrk1982]

mA = 1000, nA = 1010, mB = nB = 1000,
lda = max(1,mA);
ldb = max(mA, nA);

For matlab, I compare it with the following:
[q r] = qr(A'); s = q*(r' \ B );

This is because I think Matlab has a bug in its minimum norm solution in Linux.

[jrk1982]

Was the information I gave sufficient? Please let me know if you need to know anything else... Is the issue reproducible by you?

[john]

We have investigated GELS thoroughly and have been unable to produce any incorrect results for these test conditions. Because least squares solutions are non-unique, it is not safe to assume that a given method (whether Matlab's left divide or your QR method) will produce identical results. What matters is that norm(A*x-b) is minimal - that is the guarantee given by GELS. It's hard to determine if and when Matlab calls this function under the hood. If you can demonstrate that GELS is not producing a minimum norm solution then we will continue investigating. For any future reports, it would be very helpful to see how you are comparing the two solutions - ie norm(us-them) etc, as well as all setup and/call code.

Thank you,
John

[jrk1982]

Hi John

For the case under consideration (M < N , under constrained), full row-rank A, in a system
Ax=b, the unique, minimum euclidean norm solution(min||x||) is given by pinv(A)*b, which we can prove to be equal to q*(r'\b), where A' = qr. Hence, one would expect a unique minimum norm solution, across methods.

The solution we are looking for, is not the one for min(||Ax-b||), which is typically the case (if M > N, an over constrained system).


In my test, I'm looking at the maximum absolute value of the error, between the qr solution and the cula sgels solution.

--Karthik

[jrk1982]

The other observation is that this happens for multiple RHS case. The results are accurate for a single RHS vector.
A =

0.1209 0.8746 0.5964 0.6757
0.6432 0.3369 0.4949 0.4142

>> x

x =

0.0892 0.5664
0.5626 0.1655

sgesv solution from cula:

b =

0.7732 -0.2930
-0.1973 0.6134
0.2211 0.2504
0.0540 0.3765

matlab qr solution:
b1 =

0.7732 -0.2153
-0.1973 0.3942
0.2211 0.1490
0.0540 0.2350


norm(b(2,: ))

ans =

0.6444

>> norm(b1(2,: ))

ans =

0.4409

[jrk1982]

I found the bug in my code. The GELS code works fine!!

Thanks
--Karthik

[john]

Good to hear, thanks :)