cond (and normestinv) return not consistent values

[rfourquet]

Sorry for the bad phrasing of the title.
The following reproduces a test case at line 1102 test/sparsedir/sparse.jl, which fails with PR JuliaLang/julia#14772:

julia> srand(5)
       Ac = sprandn(20,20,.5) + im* sprandn(20, 20,.5)
       cond(Ac,1) - cond(full(Ac),1)
30.262691632830524

julia> length(Set([Base.SparseArrays.normestinv(Ac) for x in 1:50]))
3

[andreasnoack]

@mfasi What is known about the magnitude of the error? This one seems a bit large, but it matched the result I see in MATLAB with the same matrix.

[rfourquet]

The test was commented out in JuliaLang/julia#14772, it should be re-enabled when this issue is fixed.

[KristofferC]

FWIW I had a similar (same) problem here half a year ago: JuliaLang/julia#12734 (comment). It failed when I added a test above it so I guess it is sensitive to the random seed.

[jiahao]

Generally condition number computations can be very imprecise, being loose upper bounds on the true condition number that are easy to compute rather than a precise estimate.

However, in this case, it looks like the discrepancy is in the computation returned by LAPACK.gecon!.

julia> cond(full(Ac), 1) #Dense case, uses the LU factorization cond(lufact(full(Ac)), 1), calling LAPACK.gecon!
100.3697356478622

julia> norm(full(Ac), 1)
19.86009075965159

julia> norm(inv(full(Ac)), 1)
6.577634959558674

julia> norm(full(Ac), 1) * norm(inv(full(Ac)), 1) #!= cond(full(Ac), 1) !
130.6324272806925

julia> cond(Ac, 1) #Essentially the same as the preceding
130.63242728069272

julia> SparseMatrix.normestinv(Ac) #Estimate is essentially identical to the exact answer
6.5776349595586865

julia> norm(Ac, 1) #Essentially identical to dense norm
19.86009075965159

julia> SparseMatrix.normestinv(Ac) * norm(Ac, 1) # = cond(Ac, 1)
130.63242728069272

The norm of A passed into LAPACK.gecon! is correct, being norm((A[:L] * A[:U])[A[:p],:],p) = 19.86009075965159.

@andreasnoack thoughts?

[andreasnoack]

Interesting. According to the LAPACK manual, xGECON is based on

N. J. HIGHAM, FORTRAN codes for estimating the one-norm of a real or complex matrix, with applications to condition estimation, ACM Trans. Math. Softw., 14 (1988), pp. 381-396.

which reports that the procedure is doing well for random matrices, but badly for some non-random matrices. In Applied Numerical Linear Algebra, Demmel writes that a good estimate of norm(inv(A),1) "...should be within a factor of ten", since "This is all one needs for an error bound which tells you about how many decimal digits of accuracy that you have..." so we should probably just adjust the test.

[jiahao]

It looks like xGECON is known to predict arbitrarily bad underestimates of the condition number for some pathological matrices (pp. 385-6 of (pdf)). So AFAICT the underestimation of the condition number is not unexpected and there's no reason to suspect that any of the algorithms are broken.

[andreasnoack]

True, but this matrix doesn't seem to be similar to the problematic matrices in the paper. The matrix in this issue is random and has a small condition number. I agree that the algorithms are probably not wrong so adjusting the tests and updating the documentation is probably the right solution.

Apparently, it is conjectured that no bounds exist for cheap estimates of the condition number, but that in "most cases" it is within a factor of ten.

[jiahao]

All I meant to say is that I was wrong in thinking that condition number estimates always overestimate.

Perhaps we could solicit feedback from @higham to be sure that our understanding so far is correct.

[higham]

I don't see anything of concern above.
The LAPACK estimator always produces a lower bound, modulo roundoff.
It's the ratio cond(Ac,1)/cond(full(Ac),1) that is relevant (but not reported above), and it is typically less than about 3.

[jiahao]

Thanks @higham for the comment, much appreciated.