Page **1** of **1**

### culaDeviceSgesv accuracy, for M

Posted:

**Tue Sep 08, 2009 11:54 am**
by **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.

### Re:culaDeviceSgesv accuracy, for M

Posted:

**Tue Sep 08, 2009 12:02 pm**
by **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.

### Re:culaDeviceSgesv accuracy, for M

Posted:

**Tue Sep 08, 2009 12:15 pm**
by **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

### Re:culaDeviceSgesv accuracy, for M

Posted:

**Tue Sep 08, 2009 12:46 pm**
by **john**

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

### Re:culaDeviceSgesv accuracy, for M

Posted:

**Tue Sep 08, 2009 1:10 pm**
by **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.

### Re:culaDeviceSgesv accuracy, for M

Posted:

**Thu Sep 10, 2009 6:57 am**
by **jrk1982**

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

### Re:culaDeviceSgesv accuracy, for M

Posted:

**Thu Sep 10, 2009 7:49 am**
by **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

### Re:culaDeviceSgesv accuracy, for M

Posted:

**Thu Sep 10, 2009 12:15 pm**
by **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

### Re:culaDeviceSgesv accuracy, for M

Posted:

**Thu Sep 10, 2009 12:21 pm**
by **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

### Re:culaDeviceSgesv accuracy, for M

Posted:

**Fri Sep 11, 2009 7:49 am**
by **jrk1982**

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

Thanks

--Karthik

### Re:culaDeviceSgesv accuracy, for M

Posted:

**Fri Sep 11, 2009 12:41 pm**
by **john**

Good to hear, thanks :)