MATLAB – Why Cholesky Factorization and Forward Substitution Might Be Less Accurate than Inversion

Question

I recently asked this question asking for an efficient way to compute the Mahalanobis distance (without calculating the inverse). The accepted solution was to use the Cholesky factorization and forward selection.

However, it seems this is less accurate in an important case. Say that $p>n$ and that one has the following iterative algorithm (Abramovich estimator)

$$V_{(n+1)}=\sum_{i=1}^n \frac{x_ix_i^T}{x^T V_{(n)}^{-1}x}+\epsilon \cdot I$$

Note that the $\epsilon$ allows each iterative to be invertible in the next step.

However, when in my experiments in MATLAB I have seen that while Cholesky factorization is indeed faster than computing the inverse, the solution involving the inverse is more accurate. Is this anomalous or is this a well known phenomena? I had understood that the Cholesky factorization should be much more accurate in addition to being more efficient than the inverse.

Answer 1

For most cases Cholesky factorization should be a faster and more numerically stable method for solving a linear system of equation such as $Ax=b$ given that $A$ is describing a positive definite matrix. The standard workhorse behind the solution of linear systems is the QR decomposition; it does need the system $A$ to be positive definite or even square.

Some difference in performance might be found in relative small matrices as direct analytical solutions might be employed but this usually is not the case for systems larger than $3\times3$ or $4\times4$. In general in terms of performance Cholesky decomposition is approximatelly twice as fast as LU decomposition, LU decomposition is approximatelly twice as fast as QR decomposition and QR decomposition is approximatelly twice as fast as SVD decomposition. All them have different conditions that need to be so they are applicable but all them can be used to solve a linear system.

Given you are using MATLAB, MATLAB's mldivide operator will automatically make certain checks and compute the solution given the optimal decomposition from the ones described above (and some additional too, eg. check if you have a Hessenberg matrix and use a banded solver) . Powering a matrix (^-1) is a rather inefficient way to solve a linear system of equations. It will most probably resort to do the singular value decomposition of the matrix. This is stable numerically but very slow.