culaDeviceSgesv accuracy, for M
11 posts
• Page 1 of 1
culaDeviceSgesv accuracy, for M
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.
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
Re:culaDeviceSgesv accuracy, for M
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
- Administrator
- Posts: 587
- Joined: Thu Jul 23, 2009 2:31 pm
Re:culaDeviceSgesv accuracy, for M
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
I call it as:
culaDeviceSgels('N',
mA,
nA,
nB,
(culaDeviceFloat*) d_in,
lda,
(culaDeviceFloat*) d_out,
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
Re:culaDeviceSgesv accuracy, for M
I'll need some more details to continue. What are your mA, nA, nB and which Matlab function are you comparing to?
- john
- Administrator
- Posts: 587
- Joined: Thu Jul 23, 2009 2:31 pm
Re:culaDeviceSgesv accuracy, for M
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.
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
Re:culaDeviceSgesv accuracy, for M
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
Re:culaDeviceSgesv accuracy, for M
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
Thank you,
John
- john
- Administrator
- Posts: 587
- Joined: Thu Jul 23, 2009 2:31 pm
Re:culaDeviceSgesv accuracy, for M
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
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
Re:culaDeviceSgesv accuracy, for M
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
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
Re:culaDeviceSgesv accuracy, for M
I found the bug in my code. The GELS code works fine!!
Thanks
--Karthik
Thanks
--Karthik

Senior Developer, AccelerEyes
PhD Student, Dept. of ECE,
UC Santa Barbara
- jrk1982
- Posts: 26
- Joined: Thu Aug 13, 2009 8:48 pm
Re:culaDeviceSgesv accuracy, for M
Good to hear, thanks :)
- john
- Administrator
- Posts: 587
- Joined: Thu Jul 23, 2009 2:31 pm
11 posts
• Page 1 of 1
Who is online
Users browsing this forum: No registered users and 1 guest