rcond warning differs from computed value of rcond

Suppose I have a matrix M.

M = [1, 2, 3; 3, 4, 5; 5, 6, 7]

M =
     1     2     3
     3     4     5
     5     6     7

This is obviously a singular matrix, so the following orange warning is no surprise if I try to invert the matrix.

>> inv(M)
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  7.031412e-18. 

I can compute the value of RCOND displayed in the warning with the MATLAB command rcond:

>> rcond(M)
ans =
   7.0314e-18

So now I want to solve a large system Ax=b where A is sparse (which is 5185x5185 and I unfortunately can't easily share) and I get the singular matrix warning.

>> A\b
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  6.163871e-22.

My question is, why does the following computation not agree with what MATLAB outputs in the warning?

>> rcond(full(A))
ans =
   6.7160e-10

Any ideas on how I can isolate this issue? I have a hard time believing this large discrepancy is due to using sparse matrices, but have no other leads on what could be happening.

Unfortunately this forum is flooded with people asking why their specific problem is ill-conditioned so searching did me no good.

[Bjorn Gustavsson]

A quick scamble gives me nothing much to go after except this little pebble of a lead from the "online documentation" of rcond:

"rcond is a more efficient but less reliable method of estimating the condition of a matrix compared to the condition number, cond."

So my question(hint?) would be to check what you get if you use cond?

[Andrew Shi]

Thanks for your comment Bjorn. Unfortunately I get roughly the same thing, that rcond and 1/cond are about equal.

>> rcond(full(A))
ans =
   6.7160e-10
>> cond(full(A), 1)
ans =
   1.7786e+09
>> 1/ans
ans =
   5.6225e-10   

It might be worth nothing that cond and rcond both take way longer to compute than it does for MATLAB to spit out that singular matrix warning.

>> tic; cond(full(A), 1); toc
Elapsed time is 4.616524 seconds.
>> tic; rcond(full(A)); toc
Elapsed time is 1.177728 seconds.

[Bjorn Gustavsson]

If I recall correctly there is svds that (I used looong ago) that can calculate some of the singular values of sparse matrices. If you use that to calculate the largest and the smalles for your sparse, does that give the same as cond?

[Andrew Shi]

Interesting. I might have to brush up on some of my ill-conditioning from numerical linear algebra, but that extremely dominant largest singular value really scares me.

Still though, MATLAB had to come up with that number I can't seem to match somewhere, I'm just completely perplexed as to how to dig into it since it's just a warning when I do backslash and not a function I can step into.

K>> svds(A, 4, 'largest')
ans =
   1.0e+04 *
   2.683270833904093
   0.000102082444469
   0.000069453011130
   0.000069378920438
K>> svds(A, 4, 'smallest')
ans =
   1.0e-04 *
   0.399225026087114
   0.396332797639595
   0.389051148324654
   0.387694825665328

[Bjorn Gustavsson]

Ok, that gives you another estimate of the condition-number (~6.7e8 - just for laughts and giggles?)

When it comes to ill-conditioned linear inverse problems I'm an enthuastic proponent of explicit calculations of the svd - just to get the entire singular-value spectrum and both the data and model-space eigenvectors. Since I work in physics the problems I solve always have noise so some filtering/damping of solutions are typically required and a standard Tikhonov solution (0th or 2nd order) is my first stab at these type of problems. A useful matlab-toolbox for this is regtools - it has a good documentation too. (In short: It is not the big singular values that is the problem - it is all the very small ones that lead to noise-amplification)

[Andrew Shi]

Thanks Bjorn. That makes sense, I certainly do need to think about my problem some more.

This was my first time asking a question on MATLAB answers, and I really appreciate the responsiveness both of you!

Answer 1

[Bruno Luong]

When your matrix is ill-conditionned, everything computing that is related to the subspaces of the smallest eigen values are affected by numerical noise and truncation, including the estimation of the smallest eigen value.

(However the largest engen space is still stable numerically).

So when the condition number is large, the estimation of the condition number itself become unstable. Since the condition number is estimated as the ratio of the largest eigen values (stable) and the smallest ones, which is instable.

Rule of thumb if true your condition number is >= 1e10, the estimate of condition number by various methods become meaningless.

Bellow that it's OK.

[Andrew Shi]

Thanks Bruno. I can certainly buy what you're saying as an explanation for the discrepancy. However, I still wonder how MATLAB computes the estimate for RCOND that they print in the warning.

[Bruno Luong]

MATLAB function RCOND uses LAPACK

It's based on 1-norm of the matrix and its inverse. More or less it computes this, without explicitly computing inv(A)

rcond(A) = 1 /(max(sum(abs(inv(A)),1))*max(sum(abs(A),1)))

It seems RCOND estimated by INV use CONDEST. You can see in the documentation the reference paper.

[Andrew Shi]

I'm a bit confused by this answer. I'm looking for something I can type into the command window that's a function of my matrix A that can get the same value of RCOND that's printed in the warning.

But as you state, in your answer, I can't actually compute that since there's an inv(A) in there.

Could you clarify the meaning of your last sentence?

[Bruno Luong]

MATLAB selects appropriate algorithm according to the form of your matrix, whereas it dense or sparse. When you do '\' there is a huge decision tree that kicks in to decide which algo to use, and possibly estimate the RCOND differently.

The bottom line message I want pass to you is DO not try to reproduce MATLAB since it's complicated. You should use one of the function RCOND, CONDEST and threshold ot at 1e-10 (or 1e10) BEFORE MATLAB warning kicks in.

But if you want at all cost to reproduce exactly MATLAB behavious, then look at various info I pointed to you and do your own reverse engineering, I doubt anyone can have exact details, unless they works for TMW.

[Andrew Shi]

Fair enough. Thanks for your answers Bruno, I appreciate it.

[Bjorn Gustavsson]

To add my recommendation to the argument (fwiw): When you have an ill-conditioned/ill-posed/mixed-determined problem (and it is not so hard-connected to a limiting value of the condition-number of the matrix) you should preferably do your own regularization, since you know what regularizations are appropriate for your problem. In my problems it is typically appropriate with smooth solutions, which I can regularize towards using 2nd order Tikhonov regularization, and adapt the regularization to the estimated noise-level of the data. These considerations matlabs \ or pinv cannot know about, therefore I have to do that work.

[Bruno Luong]

Unshamed self-promote: I have create a simple PSSEUDOINVERSE class based on Tikhonov regularization https://www.mathworks.com/matlabcentral/fileexchange/25453-pseudo-inverse

Users still have to decide how strong the regularization he/she wants. Selecting the right regularization is quite hard to make it automatic and robust for all kind of problem.