Dealing with (near) rank deficiency efficiently?

[NicolasL-S]

Maybe this issue should be directed towards LinearAlgebra.jl.

Iterative solvers require lots of linear solves Ax=b, but these are often quite useless when A is ill-conditioned. Looking at some packages, I have seen some people handling the problem by checking the condition number doing something like:

if cond(A) < threshold
    solve Ax = b
else
    do something else
end

The problem is that computing cond(A) is often much longer than doing the actual solve. Is there a faster way of doing this? Also, LAPACK will still spend time checking whether A is full-rank. Could this be avoided?

An alternative could be

try
    solve Ax = b
catch
    do something else
end

Here, the problem is that a near rank-deficient A will not error (and we have no control over the tolerance used to assess whether solve Ax = b has failed or not). (also, try catch takes time and should generally be avoided as I've heard).

[NicolasL-S]

To those interested, I have solved my issue in a very easy way.

1. Compute the svd of the matrix.
2. Compute the condition number = largest / smallest singular value.
3. Solve the system using the same svd.

[oscardssmith]

Just for the record, this isn't a great solution since the svd is a relatively expensive way of solving a linear problem. We should be able to get something better for this.

[vpuri3]

Computing the SVD indeed very expensive, and I wont recommend that solution as you can't tell apriori how long this method would take. Unfortunately, there is no magic pill to numerical linear algebra, and there are no general purpose methods for solving large systems.

Krylov methods aren't that great either unless you have a preconditioner. As the name suggests a preconditioner "conditions" the linear problem, ie solves a problem with a reduced condition number allowing krylov solvers to converge orders of magnitude faster.

[vpuri3]

In general, if you know the answers to questions like "what is the nature of your linear system? is it symmetric positive definite? does it have a nullspce? does it emerge from a PDE? What aspect of the system is causing the large condition number? is it diagonal dominant? what is the ratio of the largest to smallest entries in the matrix?", it is possible to recommend a reasonable preconditoner, or write a custom solving algorithm.

And if your system is truly without structure, then you can always try GMRES with some sort of sparse factorization as preconditioner.

[ChrisRackauckas]

If this is a more robust solver, RobustSVD() might be a nice algorithm to have around. Not every algorithm is about being fast: some are about being stable. Speed vs stability is a trade-off and a good library should have options that let the user choose how to best use the trade-off there.

[ChrisRackauckas]

That pseudocode seems a bit off though @NicolasL-S. You compute the condition number but then don't use it?

[NicolasL-S]

Thanks for the help! I was writing my own Anderson method for fixed-point iterations (I know some already exist in Julia). Similar to the LBFGS, the method accumulates past iterates. Sometimes the matrix A (dense, but not that large) may become quite ill conditioned. A way to address this is keeping only recent iterates when its condition number becomes too high. A description of the algorithm would go something like

while the condition number of A is too high:
    compute the SVD of A
    use the SVD to compute the condition number
    if the condition number is too high, reduce the size of A by dropping the oldest iterate
end
Use the svd to solve Ax=b

That said, I'm not sold on this solution yet (maybe using QR and not reducing A would also work well, I haven't tried yet). Since the goal is accelerating convergence, this linear solve (and computation of the condition number) shouldn't start dominating in terms of compute time. I see your point about the SVD being slow. I did it this way because this is what is implemented in dense.jl for the 2-norm condition number, so I figured the same time would be needed using cond(A). Faster solutions would definitely be better. Note also that precision may be important since numerical errors accumulate iteratively and slow convergence.

[ChrisRackauckas]

I don't think there's something actionable here? We now have a system for specifying conditioning, and we have SVDFactorization, and I think that's all of the tools here. If there is a RobustSVD algorithm (which "reduce the size of A by dropping the oldest iterate" doesn't sound too generalizable) then it would be nice to see, but the description sounds like it's a solution for a specific problem that is essentially changing the problem, so I don't see a robustsvd here.