Poor accuracy of dgecon

[maxaehle]

Description
From the docs:

DGECON estimates the reciprocal of the condition number of a general
real matrix A, in either the 1-norm or the infinity-norm, using
the LU factorization computed by DGETRF.

An estimate is obtained for norm(inv(A)), and the reciprocal of the
condition number is computed as
RCOND = 1 / ( norm(A) * norm(inv(A)) ).

DGECON seems to compute fairly bad estimates of the reciprocal 1-norm condition numbers when applied to the LU factorization of the 2x2 matrix

A(p) = [ 1.0    p  ]
       [ 0.0   1.0 ]

with 0<p<1, and good estimates outside that range.

Note that the inverse of A(p) is A(-p), and both have a 1-norm of 1+abs(p), so the reciprocal of the 1-norm condition number is 1/(1+abs(p))^2.

Minimal example: dgecon_example.c.txt
Use plot_cond.py.txt for reproducing the above plot.

Note that the LU decomposition with dgetrf works well for all -3<p<3, plotted with
plot_lu.py.txt:

Also the anorm computed with dlange is fine for all -3<p<3.

Checklist

• I've included a minimal example to reproduce the issue
• I'd be willing to make a PR to solve this issue

I'd be willing to investigate the accuracy problem further if it is not yet understood and if a better accuracy appears desirable (given that dgecon only computes an estimate without specified tolerances).

[maxaehle]

Related: scipy/scipy#21620

[langou]

Hi @maxaehle,

Thanks for reporting this.

How bad can the condition number estimator be?

See Sheung Hun Cheng and Nick Higham [0]: "The price to pay for using an estimate instead of the exact condition number is that it can sometimes be a poor estimate. Experiments in [6] show that the underestimation is rarely by more than a factor of 10 (the estimate is, in fact, a lower bound), which is acceptable in practice as it is the magnitude of the condition number that is of interest. However, counterexamples for which the condition numbers can be arbritrarily poor estimates exist [6], [7]. Moreover, when the accuracy of the estimates becomes important for certain applications [9], the moethod does not provide an obvious way to improve the estimate."

[0] Sheung Hun Cheng and Nicholas Higham. Implementation for LAPACK of a Block Algorithm for Matrix 1-norm Estimation. LAPACK Working Note 152. August 13 2001.

[6] Nicholas Higham. FORTRAN codes for estimating the one-norm of a real or a complex matrix, with applications to condition estimation. (Algorithm 674.) ACM Trans. Math. Software, 14(4):381-396, December 1988.

[7] Nicholas Higham. Experience with a matrix norm estimator. SIAM J. Sci. Stat. Comp., 11:804-809, 1990.

[9] Nicholas Higham and Françoise Tisseur. A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra. SIAM J. Matrix Anal. Appl., 21(4):1185-1201, 2000.