feature request: cond() for sparse matrices

[alsam]

Julia do have condition number computation cond algorithm for dense matrices but lacks it for sparse matrices:

julia> x = sprandn(7,7,0.3);

julia> cond(x)
ERROR: no method cond(SparseMatrixCSC{Float64,Int64})

julia> y = rand(7,7)
7x7 Array{Float64,2}:
 0.105457  0.820211    0.526384    0.257885  0.0951214   0.619188  0.11465  
 0.328247  0.00401575  0.417114    0.747672  0.167338    0.132277  0.298542 
 0.624044  0.564439    0.00227799  0.883262  0.335541    0.779947  0.93982  
 0.231413  0.63107     0.523263    0.823538  0.452903    0.129218  0.0104342
 0.562113  0.536426    0.167322    0.565356  0.00539442  0.477658  0.522652 
 0.868369  0.279572    0.101924    0.435134  0.483316    0.47326   0.655955 
 0.492717  0.204044    0.230936    0.113277  0.912932    0.393085  0.269244 

julia> cond(y)
592.9028870805738

Thanks!

[tkelman]

Slow but accurate: cond(full(x))

Probably faster for large matrices but subject to iterative tolerances:

sqrt(eigs(x'*x; nev=1, which="LM", ritzvec=false)[1][1] / 
  eigs(x'*x; nev=1, which="SM", ritzvec=false)[1][1])

(denominator may have wrong sign if matrix is singular)

Edit: might be better to do what Matlab's svds does, and use [spzeros(size(x,1),size(x,1)) x; x' spzeros(size(x,2),size(x,2))] instead of the matrix product

[alsam]

Ok, thanks for recommendation!

[andreasnoack]

We should consider if the condition number estimates for the one an inf norms that LAPACK do can be generalised to the sparse case. It appears that MATLAB does something like that.

[tkelman]

Yep, would be good to have Hager-Higham condition estimates for general matrices. There's some pseudocode here http://www.cse.psu.edu/~barlow/cse451/Hager.ps

[ViralBShah]

UMFPACK does provide a rough estimate of the condition number:

Info [UMFPACK_RCOND]:  
A rough estimate of the condition number, equal
    to min (abs (diag (U))) / max (abs (diag (U))), or zero if the
    diagonal of U is all zero.

Would this be a good enough start?

[tkelman]

that's... very rough

[jiahao]

Indeed. Some of these estimates (even in LAPACK) I've found to be practically useless.

[gwhowell]

For condition estimation, another possible algorithm is being championed by Sivan Toledo
Here's a link.
http://www.eecs.berkeley.edu/~odedsc/seminars/scmc-seminar-spring2013_files/Sivan-Toledo-Slides.pdf

Allows estimate of conditioning for non square matrices ?

[gwhowell]

http://www.eecs.berkeley.edu/~odedsc/seminars/scmc-seminar-spring2013_files/Sivan-Toledo-Slides.pdf

interesting -- much faster than computing the SVD, slow sometimes .. he does survey other methods ..

[jiahao]

@gwhowell Thanks, these look like pretty interesting slides.

The basic idea has been around for awhile - Parlett's and Saad's books (amongst others) have described how you can get estimates of eigenvalues (and hence the condition number) from studying the Lanczos vectors that underly iterative solvers. The new idea here seems to be his use of the forward error to study the convergence rate of LSQR. Normally the forward error is unknown, but the idea here is to simply pick a problem with a known solution and compute the Rayleigh quotient for the smallest singular value when it is converged.

The usual caveats of iterative methods apply. The main problem (which is acknowledged in the slides but not addressed) is that the convergence rate depends strongly on how well separated the singular value you are interested in is from the other singular values. Many matrices do not have good spectral gaps for the smallest singular value, and so estimating the condition number remains a difficult problem in the general case.

This could be a feature request for IterativeSolvers.jl

[poulson]

For what it's worth, the de facto method right now for one-norm estimation is the block variant of the mentioned Hager-Higham algorithm: http://epubs.siam.org/doi/abs/10.1137/S0895479899356080?journalCode=sjmael

[gwhowell]

Hi Jack,
I think the Hager-Higham algorithm you mention is the one used in octave
and matlab -- the operator is called
condest

[edit: linke to documentation: http://octave.sourceforge.net/octave/function/condest.html]

[StefanKarpinski]

Please don't quote the copyrighted documentation of other projects here. It might qualify as fair use, but it's better not to risk any copyright infringement.

[poulson]

@gwhowell Agreed, that's what I was implying when I said "de facto". From what I recall from my last brief conversation with Nick Higham, MATLAB defaults to setting the blocksize to 2.

[jiahao]

Please don't quote the copyrighted documentation of other projects here. It might qualify as fair use, but it's better not to risk any copyright infringement.

attn: @gwhowell - we've asked you on previous occasions not to do this.

[poulson]

@StefanKarpinski @jiahao I assume that quoting documentation from permissively licensed (e.g., New BSD) projects is okay?

[tkelman]

Better to just link to things in general, if they're more than a few lines

[mfasi]

I implemented the Hager-Higham algorithm in Julia (a tentative gist) and would be happy to add it to base, but I'm not sure what the best way to do this will be. If I run cond(A) with a A = sprandn(n,n,d) I get the following

julia> cond(A,2)
ERROR: MethodError: `svdvals!` has no method matching svdvals!(::Base.SparseMatrix.SparseMatrixCSC{Float64,Int64})
 in cond at linalg/dense.jl:494
 in cond at linalg/dense.jl:493

julia> cond(A,1)
ERROR: MethodError: `cond` has no method matching cond(::Base.SparseMatrix.UMFPACK.UmfpackLU{Float64,Int64}, ::Int64)
Closest candidates are:
  cond(::Number, ::Any)
  cond{T<:Union{Complex{Float32},Float64,Float32,Complex{Float64}},S}(::Base.LinAlg.LowerTriangular{T<:Union{Complex{Float32},Float64,Float32,Complex{Float64}},S}, ::Real)
  cond{T<:Union{Complex{Float32},Float64,Float32,Complex{Float64}},S}(::Base.LinAlg.UnitLowerTriangular{T<:Union{Complex{Float32},Float64,Float32,Complex{Float64}},S}, ::Real)
  ...
 in cond at linalg/dense.jl:499

julia> cond(A,Inf)
ERROR: MethodError: `cond` has no method matching cond(::Base.SparseMatrix.UMFPACK.UmfpackLU{Float64,Int64}, ::Float64)
Closest candidates are:
  cond(::Number, ::Any)
  cond{T<:Union{Complex{Float32},Float64,Float32,Complex{Float64}},S}(::Base.LinAlg.LowerTriangular{T<:Union{Complex{Float32},Float64,Float32,Complex{Float64}},S}, ::Real)
  cond{T<:Union{Complex{Float32},Float64,Float32,Complex{Float64}},S}(::Base.LinAlg.UnitLowerTriangular{T<:Union{Complex{Float32},Float64,Float32,Complex{Float64}},S}, ::Real)
  ...
 in cond at linalg/dense.jl:499

which I think isn't better than the original ERROR: no method cond(SparseMatrixCSC{Float64,Int64}).

I think what we should have is:

• A condest method to estimate the 1-norm condition number for both dense and sparse matrices
• A cond method for AbstractMatrix (in generic.jl rather than dense.jl, maybe) which computes cond(A) if A is dense and returns a message to help the user choose between cond(full(A)) and condest(A).

Does this sound reasonable? In case it does, I'll be happy to adapt the gist to base and open a pull request. Comments on the code would be appreciated as well.

[tkelman]

I've changed my mind since April 2014, while I still think this code would be really good to have cleaned up and easily available somewhere, I'm not sure it absolutely needs to be in base unless it's important enough to be a dependency of all other code written in the language. Especially now with package precompilation, what about making a package for it, or contributing it to somewhere like https://github.com/andreasnoack/LinearAlgebra.jl (which is not yet registered, but probably will be sooner or later)?

[mfasi]

@tkelman Ok, thank you, I'll clean it up and then see where to put it. I might need it at some point to implement one of the functions in JuliaLang/julia#5840, but I'm not sure yet.
However, I think that cond(A) should return an informative message rather than an apparently unrelated error when A is a SparseMatrix. Maybe it would be enough to do something similar to what inv(A) does?

[andreasnoack]

As long as we haven't moved larger pieces of the sparse functionality out of base I think the condition number estimates for sparse matrices belong in base. It should just be called cond and without the est part. The LAPACK condition number subroutine returns estimates and you could also argue that the computed spectral condition number is an estimate (just with a small error).

[tkelman]

Let's be honest, condition estimates are a pretty niche feature. And I don't think we should call it cond if it's an estimate.

I guess we're calling the Lapack 1 and inf norm estimators in the dense case for cond(A, 1) and cond(A, Inf), but we should more clearly document that. Maybe we could put this code into base for the 1 or inf norm estimates in the sparse case, but we should really clarify the documentation that these are estimates.

[andreasnoack]

Improvements to the documentation are always welcome and I agree that it would great to clarify how the different condition numbers are computed. Preferably with links to the relevant sources/papers.

[gwhowell]

Along with a link to the relevant papers, the documentation should point
out what condition number is being estimated. For the Higham Hager
algorithm, that's L_1 condition number for example (don't know that the
current dense cond() is estimating that same number, if not, that would be
argument for using two different names for the sparse and dense cases)

On Fri, Jul 31, 2015 at 9:16 AM, Andreas Noack [email protected]
wrote:

Improvements to the documentation are always welcome and I agree that it
would great to clarify how the different condition numbers are computed.
Preferably with links to the relevant sources/papers.

—
Reply to this email directly or view it on GitHub
#104.

[andreasnoack]

In general, I don't think we should have different functions for different algorithms for the same quantity. We should pick a good default and document it. If we'd like to give choices (which I'm considering for for eigen and singular values decompositions) I think we should use a keyword argument to select the algorithm.

[mfasi]

In this case, cond(A, p) where p = 1, 2 or Inf computes the p-condition number of a matrix, the default being p = 2. The Hager-Higham algorithm estimates the 1 (and thus Inf) condition number, the case p = 2 can be implemented as the power method on A' * A.

[andreasnoack]

I wouldn't bother with the spectral condition number for now. It would be great just to have 1 and Inf. The smallest singular value would be pretty inaccurate if calculated with A'A. In principle, we should be able to calculate the spectral condition number with svds, but it appears that right now our implementation only allows us to compute the largest values.

[andreasnoack]

Closed by JuliaLang/julia#12467