General CULA Dense (LAPACK & BLAS) support and troubleshooting. Use this forum if you are having a general problem or have encountered a bug.

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.

Senior Developer, AccelerEyes
PhD Student, Dept. of ECE,
UC Santa Barbara
jrk1982

Posts: 26
Joined: Thu Aug 13, 2009 8:48 pm

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.
john

Posts: 587
Joined: Thu Jul 23, 2009 2:31 pm

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

I call it as:
mA,
nA,
nB,
lda,
ldb);

Regards
--Karthik

Senior Developer, AccelerEyes
PhD Student, Dept. of ECE,
UC Santa Barbara
jrk1982

Posts: 26
Joined: Thu Aug 13, 2009 8:48 pm

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

Posts: 587
Joined: Thu Jul 23, 2009 2:31 pm

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.

Senior Developer, AccelerEyes
PhD Student, Dept. of ECE,
UC Santa Barbara
jrk1982

Posts: 26
Joined: Thu Aug 13, 2009 8:48 pm

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

Senior Developer, AccelerEyes
PhD Student, Dept. of ECE,
UC Santa Barbara
jrk1982

Posts: 26
Joined: Thu Aug 13, 2009 8:48 pm

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
john

Posts: 587
Joined: Thu Jul 23, 2009 2:31 pm

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

Senior Developer, AccelerEyes
PhD Student, Dept. of ECE,
UC Santa Barbara
jrk1982

Posts: 26
Joined: Thu Aug 13, 2009 8:48 pm

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

Senior Developer, AccelerEyes
PhD Student, Dept. of ECE,
UC Santa Barbara
jrk1982

Posts: 26
Joined: Thu Aug 13, 2009 8:48 pm

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

Thanks
--Karthik

Senior Developer, AccelerEyes
PhD Student, Dept. of ECE,
UC Santa Barbara
jrk1982

Posts: 26
Joined: Thu Aug 13, 2009 8:48 pm

Good to hear, thanks :)
john