Asked by Mathieu Walsh
on 14 Jun 2019

I have a very large sparse matrix, Amat (192611x192611), that has a density of 6.7190e-05 . When I am trying to use the backslash command with this matrix (cnew =Amat/rhs, where rhs is a column vector (192611 x 1) with small values, around 1e-13) on matlab 2019a, I get an error that says ''Warning: Matrix is singular to working precision." and the output array is full of NaN. However, when I run the same command on matlab 2011b (the code was made in 2011) or matlab 2014a, I get the same error, but I don't get these Nan values. Instead, I get small values, around 1e-14. My questions are the following :

1) Has the backslash command changed between these versions?

2) Are the values that I'm getting with the other versions reliable?

Answer by Walter Roberson
on 14 Jun 2019

Accepted Answer

1) Yes, the \ command has changed. Between those versions, the underlying LAPACK libraries were upgraded -- twice, I think.

2) No, LAPACK has gotten more reliable and less likely to give results for matrices that are effectively singular.

Mathieu Walsh
on 14 Jun 2019

Thank you. Do you know where can I get the changes that they've made on the backslash ?

Walter Roberson
on 14 Jun 2019

Sign in to comment.

Answer by Christine Tobler
on 14 Jun 2019

For more details about what backslash does for sparse matrices, use

spparms('spumoni', 1)

cnew = Amat \ rhs;

This will display information about which solver is being used. If you run this command in both MATLAB versions and copy the displayed information here, I may be able to tell what changed between versions (if it's anything else then just the LAPACK update).

Mathieu Walsh
on 15 Jun 2019

Christine Tobler
on 18 Jun 2019

Sorry for the late reply, I was unsure of the best way to do this. If you send me a direct message through MATLAB Answers (by clicking on my profile), we can exchange email addresses and you could send your file as an email attachment.

Alternatively, I can set up a secure file transfer to do this - but I would also need your email address for that.

Christine Tobler
on 18 Jun 2019

Thank you for sending me the .mat file. I have found the reason for this change - it was an intentional change in behavior, introduced in R2018b.

Before R2018b, mldivide was treating singular matrices differently if they were sparse symmetric indefinite. It skipped over the division by zero which is needed there.

>> A = sparse([1:4 1 1 3 4], [1:4 4 3 1 1], 1, 6, 6);

>> full(A)

ans =

1 0 1 1 0 0

0 1 0 0 0 0

1 0 1 0 0 0

1 0 0 1 0 0

0 0 0 0 0 0

0 0 0 0 0 0

>> x = A \ ones(6, 1)

Warning: Matrix is singular to working precision.

x =

1

1

0

0

2

2

% Compare to full matrix case

>> x = full(A) \ ones(6, 1)

Warning: Matrix is singular to working precision.

x =

NaN

NaN

NaN

NaN

NaN

Inf

% Compare to non-symmetric sparse case

>> A = sparse([1:4 1 1 3], [1:4 4 3 1], 1, 6, 6);

>> x = A \ ones(6, 1)

Warning: Matrix is singular to working precision.

x =

Inf

1

-Inf

1

Inf

Inf

% Compare if zero diagonal entries are replaced by eps

>> A = sparse([1:4 1 1 3 4 5 6], [1:4 4 3 1 1 5 6], [ones(1, 8) eps eps], 6, 6);

>> x = A \ ones(6, 1)

x =

1.0e+15 *

0.0000

0.0000

0

0

4.5036

4.5036

This inconsistency was removed for R2018b, by having the sparse symmetric indefinite case return a vector of all-NaNs in backslash if the input matrix does not have full rank.

Workaround: If a sparse matrix is singular to working precision, this is most commonly because there are some rows and/or columns which are exactly zero. For a symmetric matrix, these will be the same. By replacing the diagonals for those rows and columns with a finite number, you can get the previous behavior back.

>> A = sparse([1:4 1 1 3 4], [1:4 4 3 1 1], 1, 6, 6);

>> ind = ~any(A~=0, 2)

ind =

(5,1) 1

(6,1) 1

>> A(ind, ind) = 0.5*speye(nnz(ind));

>> x = A \ ones(6, 1)

x =

1

1

0

0

2

2

I also tried this with your original matrix and got behavior close to the previous one

>> ind = ~any(Amat~=0, 2)

ind =

(2258,1) 1

(32919,1) 1

(33002,1) 1

(146136,1) 1

(147100,1) 1

>> Amat(ind, ind) = 0.5*speye(nnz(ind));

>> x = Amat \ rhs;

>> max(abs(x))

ans =

1.2954e-13

>> all(~isnan(x))

ans =

1

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.