6

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.

Community
  • 1
Lepidopterist
  • 712
  • 5
  • 23
  • What you describe sounds rather counter-intuitive indeed. Can you please give us an actual numerical example? Please also see the implementation note I did in your previous thread. Finally just to clarify by inverse you really mean the `inv()` not the pseudo-inverse `pinv()`. – usεr11852 Apr 22 '15 at 02:10
  • I can do so soon. I note that rather than inv() I'm using (MATRIX)^-1 – Lepidopterist Apr 22 '15 at 02:49
  • That is most probably doing an SVD, so you are not even using a QR decomposition.... (I am 100% sure how MATLAB computes the power of a matrix but through an SVD is the obvious way) In any case try the operators I commented in the previous post; MATLAB have a few amazing stuff and the `mldivide` function, i.e. `A\b`, is definitely one. – usεr11852 Apr 22 '15 at 03:45
  • Is chol inefficient in MATLAB? Is mldivide doing something smarter? – Lepidopterist Apr 22 '15 at 03:55
  • No, `chol` is not inefficient in MATLAB; yes, `mldivide` might do something smarter, check the documentation of it. Solving a linear system using Cholesky is only one of the things it might do. – usεr11852 Apr 22 '15 at 04:22
  • What happened to that example? :) – usεr11852 May 07 '15 at 01:57
  • I got it to work. Maybe I should close the question. Thank you for the mldivide suggestion, it works well. – Lepidopterist May 07 '15 at 13:22
  • Cool! Later today I will write up the comments as an answer so the whole thread is easily read from other users too. – usεr11852 May 07 '15 at 15:55
  • @usεr11852saysReinstateMonic I provided a numerical example (in R though) in the [this](https://stackoverflow.com/questions/59459602/why-is-cholesky-decomposition-not-giving-me-the-same-result-as-simply-inverting) question. Although I've written it in R it should be very straightforward to read – Euler_Salter Dec 23 '19 at 18:32

1 Answers1

2

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.

usεr11852
  • 33,608
  • 2
  • 75
  • 117