DPOCON, DGECON with result

[ilya-fursov]

Hello,

I've been testing dpocon on a 10x10 symmetric positive definite matrix, and it resulted in a condition number 14.5376, which is different from Matlab result 16.4048.
When I applied dpocon to the inverse of the original matrix, the result was correct: 16.4048.
Is it the intended behaviour or a bug?

I attach c++ code which defines the test matrix and LAPACK procedure calls to demonstrate what I described above. The code #includes "lapacke.h" and "cblas.h", and should be linked with appropriate libraries. I tested it on Netlib's LAPACK and reference BLAS, as well as ATLAS BLAS on Ubuntu, and on Windows with Cygwin.

P.S. I also tested dgecon on the same matrix, with the same wrong result 14.5376.

[Jenniferfhughes]

SPOCON is a simplified interface to the JLAPACK routine spocon. This interface converts Java-style 2D row-major arrays into the 1D column-major linearized . DGECON estimates the reciprocal of the condition number of a general * real matrix A, in either the 1-norm or the infinity-norm.

[ilya-fursov]

Ok, I don't know why you mention JLAPACK here, but the original question was about a case where
condition number from dpocon has an error 13% compared to the correct condition number.

What is the reason of this error?

[Julien Langou]

Hi Ilya,

Your test code looks good to me.

Please remember that DPOCON and DGECON are condition number estimators. What matters (in general) for the condition number is the order of magnitude. Is the condition number 1, 10, 100, 1000, 10000, etc. ? Whether the condition number is 14 or 16 is not that important (in general). In any case, DPOCON returns an estimate of the condition number. It does not return the condition number. I quote Higham and Tisseur [2]: "The [...] LAPACK estimators both produce estimates that in practice are almost always within a factor [...] 3 [...], of the quantities they are estimating [3], [4], [5]." So in your case an answer between 5 and 48 would be deemed to be acceptable. The original algorithm is from Hager [1].

Note 1: If you want to compute the exact condition number, you have to compute inv(A) explicitly, and then computing the 1-norm. This is (much) more expensive than POCON.

Note 2: I have no idea why using POCON on inv(A) returns a very accurate answer, while POCON on A returns a not so accurate answer.

Julien.

[1] William W. Hager, "Condition Estimates," SIAM J. Sci. Stat. Comput. 5, 1984, 311-316, 1984.

[2] Nicholas J. 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., Vol. 21, 1185-1201, 2000.

[3] N. J. Higham, A survey of condition number estimation for triangular matrices, SIAM Rev., 29 (1987), pp. 575–596.

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

[5] N. J. Higham, Experience with a matrix norm estimator, SIAM J. Sci. Stat. Comput., 11 (1990), pp. 804–809.

[ilya-fursov]

Thanks for the references,
This makes sense now.

I was originally confused since
- LAPACK online reference and dpocon.f don't say anything about precision of the estimate (and strictly speaking almost everything we calculate on computers is just an estimate).
- Matlab gave an exact estimate for the same example.

Ilya

[Julien Langou]

Matlab gave an exact estimate for the same example.

If you use `cond`, you get the `exact` condition number.

If you use `condest`, you get same as LAPACK xPOCON (or xGECON).

As a matter of fact, if we type condest(A) in matlab, we get 14.5376; while we obtain 16.4048 for cond(A,1), cond(inv(A),1) and condest(inv(A)). Same numbers as produced by your C++ code.

Cheers,
Julien.