Backslash error ''Warning: Matrix is singular to working precision."

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 1

[Walter Roberson]

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]

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

[Walter Roberson]

If I recall correctly, they added a few more specialized detections to be more efficient. However it is unlikely you are encountering that: You are more likely encountering the fact that the LAPACK library was brought up to date.

http://www.netlib.org/lapack/improvement.html

Answer 2

[Christine Tobler]

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]

I used the command spparms('spumoni', 2)

For matlab 2019a (64 bit) :

cnew=diffusion2solve(dt,AmatAll,BmatAll,oc,os,[c;cg],im.Mesh.nodevol);

sp\: bandwidth = 186372+1+186372.

sp\: is A diagonal? no.

sp\: is band density (6.72603e-05) > bandden (0.5) to try banded solver? no.

sp\: is A triangular? no.

sp\: is A morally triangular? no.

sp\: is A a candidate for Cholesky (symmetric, real positive or negative diagonal)? no.

sp\: But A is real symmetric; try MA57.

MA57 Control Parameters

Ordering chosen: Automatic

Block Size for Level 3 BLAS: 16

Merge assembly tree nodes if number of eliminations

is less than: 16

Scaling (via symmetrized MC64): YES

Maximum iterative refinement steps: 10

Pivot Threshold: 1.000000e-02

Zero Pivot Tolerance: 1.000000e-20

MA57 Symbolic Analysis

Ordering used: METIS

Forecast number of Real entries in factors: 30868603

Forecast number of Integer entries in factors: 1281306

Forecast maximum frontal size: 2595

Number of nodes in Assembly tree: 40066

Forecast length of FACT array (real)

without numerical pivot: 39886389

Forecast length of ifact array (integer)

without numerical pivot 4990988

Length of fact required for a successful

completion of the numerical phase allowing

data compression (without numerical pivoting) 38556283

Length of ifact required for a successful

completion of the numerical phase allowing

data compression (without numerical pivoting) 4990988

Number of data compresses 6

Forecast number of floating point additions

for the assembly processes 1.367403e+08

Forecast number of floating point operations

to perform the elimination operations

counting multiply-add pair as two operations 2.821416e+10

*** Return code from MA57AD: 1

MA57 Numeric Factorization

Number of entries in factors: 30868603

Storage for real data in factors: 30868603

Storage for integer data in factors: 1281306

Minimum length of fact required for a successful

completion of the numerical phase: 38556283

Minimum length of ifact required for a successful

completion of the numerical phase: 4990988

Minimum length of fact required for a successful

completion of the numerical phase without

data compression: 39886389

Minimum length of ifact required for a successful

completion of the numerical phase without

data compression: 4990988

Order of the largest frontal matrix: 2595

Number of 2x2 numerical pivots: 0

Total number of fully-summed variables passed to

the father node because of numerical pivot: 4

Number of negative eigenvalues: 0

Rank of factorization of the matrix: 192606

Pivot step where static pivot commences and

is set to zero if no modification performed: 0

Number of data compresses on real

data structures: 0

Number of data compresses on integer

data structures: 0

Number of block pivots in factors 40066

Number of zeros in triangle of factors: 201318

Number of zeros in rectangle of factors: 1792324

Number of zero columns in rectangle of factors: 0

Number of static pivots chosen: 0

Number of floating point additions

for the assembly processes: 1.367403e+08

Number of floating point operations

to perform the elimination operations

counting multiply-add pair as two operations: 2.821416e+10

Minimum value of the scaling factor (MC64): 1.000000e+00

Maximum value of the scaling factor (MC64): 1.000000e+00

Maximum modulus of matrix entry: 8.144467e+04

*** Return code from MA57BD: 4

Number of passes to properly factor the matrix: 1

Warning: Matrix is singular to working precision.

For matlab 2011b (32 bit) :

cnew=diffusion2solve(dt,AmatAll,BmatAll,oc,os,[c;cg],im.Mesh.nodevol);

sp\: bandwidth = 186372+1+186372.

sp\: is A diagonal? no.

sp\: is band density (0.00) > bandden (0.50) to try banded solver? no.

sp\: is A triangular? no.

sp\: is A morally triangular? no.

sp\: is A a candidate for Cholesky (symmetric, real positive diagonal)? no.

sp\: But A is real symmetric; try MA57.

MA57 Control Parameters

Ordering used: AMD via MC47

Block Size for Level 3 BLAS: 16

Merge assembly tree nodes if number of eliminations

is less than: 16

Scaling (via symmetrized MC64): YES

Maximum iterative refinement steps: 10

Pivot Threshold: 1.000000e-002

Zero Pivot Tolerance: 1.000000e-020

MA57 Symbolic Analysis

Forecast number of Real entries in factors: 56843947

Forecast number of Integer entries in factors: 1346263

Forecast maximum frontal size: 4491

Number of nodes in Assembly tree: 40355

Forecast length of FACT array (real)

without numerical pivot: 83158500

Forecast length of ifact array (integer)

without numerical pivot 4990988

Length of fact required for a successful

completion of the numerical phase allowing

data compression (without numerical pivoting) 81832219

Length of ifact required for a successful

completion of the numerical phase allowing

data compression (without numerical pivoting) 4990988

Number of data compresses 0

Forecast number of floating point additions

for the assembly processes 1.555713e+008

Forecast number of floating point operations

to perform the elimination operations

counting multiply-add pair as two operations 1.176029e+011

*** Return code from MA57AD: 1

MA57 Numeric Factorization

Number of entries in factors: 56843947

Storage for real data in factors: 56843947

Storage for integer data in factors: 1346263

Minimum length of fact required for a successful

completion of the numerical phase: 81832219

Minimum length of ifact required for a successful

completion of the numerical phase: 4990988

Minimum length of fact required for a successful

completion of the numerical phase without

data compression: 83158500

Minimum length of ifact required for a successful

completion of the numerical phase without

data compression: 4990988

Order of the largest frontal matrix: 4491

Number of 2x2 numerical pivots: 0

Total number of fully-summed variables passed to

the father node because of numerical pivot: 5

Number of negative eigenvalues: 0

Rank of factorization of the matrix: 192606

Pivot step where static pivot commences and

is set to zero if no modification performed: 0

Number of data compresses on real

data structures: 0

Number of data compresses on integer

data structures: 0

Number of block pivots in factors 40355

Number of zeros in triangle of factors: 123876

Number of zeros in rectangle of factors: 1851432

Number of zero columns in rectangle of factors: 0

Number of static pivots chosen: 0

Number of floating point additions

for the assembly processes: 1.555713e+008

Number of floating point operations

to perform the elimination operations

counting multiply-add pair as two operations: 5.747331e+010

Minimum value of the scaling factor (MC64): 3.504037e-003

Maximum value of the scaling factor (MC64): 1.000000e+000

Maximum modulus of matrix entry: 8.144467e+004

*** Return code from MA57BD: 4

Number of passes to properly factor the matrix: 1

Warning: Matrix is singular to working precision.

[Walter Roberson]

I do see some differences, some of which are pure cosmetic and some not cosmetic, but it looks to me as if the same algorithm was used.

Below, lines prefixed with '<' are the R2019a output, and lines prefixed with '>' are the corresponding R2011a output

$ diff r2019a r2011b

4c4

< sp\: is band density (6.72603e-05) > bandden (0.5) to try banded solver? no.

---

> sp\: is band density (0.00) > bandden (0.50) to try banded solver? no.

7c7

< sp\: is A a candidate for Cholesky (symmetric, real positive or negative diagonal)? no.

---

> sp\: is A a candidate for Cholesky (symmetric, real positive diagonal)? no.

10c10

< Ordering chosen: Automatic

---

> Ordering used: AMD via MC47

16,17c16,17

< Pivot Threshold: 1.000000e-02

< Zero Pivot Tolerance: 1.000000e-20

---

> Pivot Threshold: 1.000000e-002

> Zero Pivot Tolerance: 1.000000e-020

19,23c19,22

< Ordering used: METIS

< Forecast number of Real entries in factors: 30868603

< Forecast number of Integer entries in factors: 1281306

< Forecast maximum frontal size: 2595

< Number of nodes in Assembly tree: 40066

---

> Forecast number of Real entries in factors: 56843947

> Forecast number of Integer entries in factors: 1346263

> Forecast maximum frontal size: 4491

> Number of nodes in Assembly tree: 40355

25c24

< without numerical pivot: 39886389

---

> without numerical pivot: 83158500

30c29

< data compression (without numerical pivoting) 38556283

---

> data compression (without numerical pivoting) 81832219

34c33

< Number of data compresses 6

---

> Number of data compresses 0

36c35

< for the assembly processes 1.367403e+08

---

> for the assembly processes 1.555713e+008

39c38

< counting multiply-add pair as two operations 2.821416e+10

---

> counting multiply-add pair as two operations 1.176029e+011

42,44c41,43

< Number of entries in factors: 30868603

< Storage for real data in factors: 30868603

< Storage for integer data in factors: 1281306

---

> Number of entries in factors: 56843947

> Storage for real data in factors: 56843947

> Storage for integer data in factors: 1346263

46c45

< completion of the numerical phase: 38556283

---

> completion of the numerical phase: 81832219

51c50

< data compression: 39886389

---

> data compression: 83158500

55c54

< Order of the largest frontal matrix: 2595

---

> Order of the largest frontal matrix: 4491

58c57

< the father node because of numerical pivot: 4

---

> the father node because of numerical pivot: 5

67,69c66,68

< Number of block pivots in factors 40066

< Number of zeros in triangle of factors: 201318

< Number of zeros in rectangle of factors: 1792324

---

> Number of block pivots in factors 40355

> Number of zeros in triangle of factors: 123876

> Number of zeros in rectangle of factors: 1851432

73c72

< for the assembly processes: 1.367403e+08

---

> for the assembly processes: 1.555713e+008

76,79c75,78

< counting multiply-add pair as two operations: 2.821416e+10

< Minimum value of the scaling factor (MC64): 1.000000e+00

< Maximum value of the scaling factor (MC64): 1.000000e+00

< Maximum modulus of matrix entry: 8.144467e+04

---

> counting multiply-add pair as two operations: 5.747331e+010

> Minimum value of the scaling factor (MC64): 3.504037e-003

> Maximum value of the scaling factor (MC64): 1.000000e+000

> Maximum modulus of matrix entry: 8.144467e+004

[Christine Tobler]

I'm seeing two differences:

1) The reordering method is METIS instead of AMD. We recently introduced the METIS reordering method, which is more efficient than AMD in many cases. Try calling

spparms('reorder', 4)

in MATLAB R2019a, this will use the previous default ordering.

2) There is no scaling in R2019a, but there was in R2011b. For now, I would assume that this is a result of the different reordering method, but that could be interesting to check.

It would make sense to me that the scaling avoids the NaN return values - I'm not sure why a different ordering would mean that no scaling is used, though.

If you would be willing to attach the matrix here, I would be interested to investigate what's going on here directly.

[Mathieu Walsh]

1) spparms('reorder', 4);

newc=Amatlul\rhs;

Warning: Unknown keyword parameter reorder.

2) I can't attach the matrix here, because it is 25 mb. Is there another way for me to send it to you?

[Walter Roberson]

The reorder keyword needs R2019a.

[Mathieu Walsh]

My bad I forgot that I had matlab 2018b on this computer. I tried the command on my other computer and I still get a vector full of NaN. I will give you the result from the spy function

[Christine Tobler]

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]

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