# Quantifying the failure of the Cholesky factorization test for indefinite matrices

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?

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!

• The exact product $A^*A$ is semidefinite, not indefinite. The matrix ATA that you pass to Cholesky -- who knows? It depends on the computational errors inside np.dot. So we cannot infer anything from this example. – Federico Poloni Mar 6 '17 at 6:35
• You are correct -- I mixed up indefinite and semidefinite. I also did not read the question carefully enough -- the interesting case is specifically the instability of the algorithm for indefinite matrices, which I did not address. – supergra Mar 7 '17 at 17:56