Quantifying the failure of the Cholesky factorization test for indefinite matrices

Question

The Cholesky factorization is the classic test to check if a matrix is positive definite. In infinite precision it is also an exact test: A matrix has a Cholesky factorization iff it is positive definite. However, in floating point arithmetic the Cholesky test is not perfect and two types of errors can (probably) be made:

1. A false positive error: The Cholesky test passes when the matrix is indefinite.

2. A false negative error: The Cholesky test fails when the matrix is positive definite.

The conditions for false negative errors is pretty well understood. For example, if a matrix $A$ is positive definite then defining $D = \hbox{diag}(A)^{1/2}$ and $A = DHD$, a false negative error is made if $\lambda_{\min}(H) \leq -n \gamma_{n+1}/(1-\gamma_{n+1})$, $\gamma_{n} = \mathcal{O}(n u)$ with $u$ being unit round off.

This is a satisfactory result because it can also be shown that the test makes no error on positive definite matrices if $\lambda_{\min}(H) > n \gamma_{n+1}/(1-\gamma_{n+1})$. For more details see theorem 10.7, page 200 of Accuracy and Stability of Numerical Algorithms By Nicholas J. Higham.

The situation seems more difficult for quantifying false positive error. On indefinite matrices the Cholesky factorization is numerically unstable so one would expect examples where an indefinite matrix (which is far from positive definite) still passes the Cholesky test. Does anyone know of an example? Can the Cholesky factorization be used in finite precision as a test without quantifying the false positive error? Is it just a highly improbable event that an indefinite matrix causes a false positive error?

Answer 1

(This answer only addresses semidefinite matrices, not indefinite matrices, so it doesn't directly answer the question.)

False positives can easily happen. Consider building the normal equation matrix for a 3D least-squares problem with only two observations.

#!/usr/bin/python3
import numpy as np
A = np.array([[ 1.5913486472486638e-07, 3.5332365422302558e-08, 2.9810777287749696e-05],
              [-4.5829836998710627e-06, 2.8828850632468817e-05, 1.8186216299018523e-05]])
ATA = np.dot(A.T, A)

print(np.linalg.cholesky(ATA))
print(np.linalg.inv(ATA))
print(np.linalg.det(ATA))

This should be guaranteed to fail, as the exact product ATA must be semi-definite. Nevertheless, the result (on my machine) is a successful Cholesky decomposition, inverse, and non-zero determinant:

[[  4.58574569e-06   0.00000000e+00   0.00000000e+00]
 [ -2.88102610e-05   1.03573186e-06   0.00000000e+00]
 [ -1.71407671e-05   3.04239216e-05   1.39200237e-11]]
[[  1.36481647e+27   2.21729573e+26  -7.54841487e+24]
 [  2.21729573e+26   3.60224284e+25  -1.22632372e+24]
 [ -7.54841487e+24  -1.22632372e+24   4.17481532e+22]]
5.40352709685e-46

Of course the problem is not Cholesky, but the inexact representation of the input matrix.

How often does this happen?

Consider building the ATA matrix from two random vectors, with elements between 0 and 1.

n_pass = 0
N = 100000
for i in range(N):
    A = np.random.rand(2, 3)
    ATA = np.dot(A.T, A)
    try:
        np.linalg.cholesky(ATA)
        n_pass += 1
    except:
        pass

print(100*n_pass/N,"% passed")

It appears that Cholesky falsely passes about 39-40% of the time!