Backslash error ''Warning: Matrix is singular to working precision."
Ältere Kommentare anzeigen
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?
Akzeptierte Antwort
Walter Roberson
am 14 Jun. 2019
0 Stimmen
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.
2 Kommentare
Mathieu Walsh
am 14 Jun. 2019
Walter Roberson
am 14 Jun. 2019
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.
Weitere Antworten (1)
Christine Tobler
am 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).
8 Kommentare
Mathieu Walsh
am 14 Jun. 2019
Walter Roberson
am 14 Jun. 2019
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
am 15 Jun. 2019
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
am 15 Jun. 2019
Walter Roberson
am 15 Jun. 2019
The reorder keyword needs R2019a.
Mathieu Walsh
am 15 Jun. 2019
Bearbeitet: Mathieu Walsh
am 18 Jun. 2019
Christine Tobler
am 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
am 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
Kategorien
Mehr zu Data Preprocessing finden Sie in Hilfe-Center und File Exchange
Produkte
am 14 Jun. 2019
am 18 Jun. 2019
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
