Lost of precision using sparse matrices?

I have two matrices J1 (sparse) = J2(full).

The dimension of the matrices are ~ 5200x2600

Then when I do:

hlm1 = (J1'*J1 + u*I)\g, where I = eye(n);

and

hlm2 = (J2'*J2 + u*I)\g, where I = eye(n);

i have after that: norm(hlm1 - hlm12, Inf) is 4.8625e-05 ...

That difference is my problem, is correct the way to use the matrice sparse ?.

Thx.

[Matt J]

Incidentally, you've given no information proving that the sparse version is less accurate. How do you know that the problem isn't with the full version?

[Matt J]

You should also be looking at percent error, rather than absolute error

 percentError = norm(hlm1 - hlm12, Inf)/norm(hlm2,inf)*100

Answer 1

[Matt J]

Precision should improve if you do

 [J;sqrt(u)*I]\g

but it's also possible that your J matrices are just very ill-conditioned. What does rcond(J2) say?

[Kei]

>> rcond(J2) Error using rcond Input must be a square matrix.

The J2 version "full" is a matrix that was not taken from a full matrix of a sparse as such, but J1 == J2.

[Matt J]

Then try rcond(J2'*J2 + u*I)

[Kei]

>> rcond(J2'*J2 + u*I) ans = 3.8290e-15

The value of u varies continuously (in each interaction), so in the above u = e-11 for example, but can be e+2 then the result is 0.046 and so...

[Matt J]

Well, it makes sense that you have poor precision when rcond=e-15. It means the matrix is almost singular.

[Kei]

But all this knowing that J1 == J2 ... makes one think that the result of hlm1 is equivalent to hlm2.

When I do hlm1 = (full(J1)'*full(J1) + u*Is)\-g; the result of norm(hlm1 - hlm2,Inf) = 6.78e-18 and for my is OK, but using 'full' .... I think it is something internal to the operation \ with sparse matrix.

[Matt J]

Well the algorithms used to solve linear equations are different for sparse and full matrices. If the solution is unstable, as it inevitably will be for ill-conditioned (J'*J+u*I), the result will be sensitive to the solution method used.