Assuming a given matrix is of the form $A^TA$, where $A$ is a real matrix of full rank, $A^TA$ is a positive definite matrix and we can write its Cholesky decomposition.

Is there a more efficient algorithm of writing the Cholesky decomposition of $A^TA$ = $L^TL$? Is it possible to directly express elements of $L$ in terms of elements of $A$ (in a non recursive fashion)?

  • 1
    $\begingroup$ You can QR factor $A$, but probably computing Cholesky from $A^TA$ is better. $A^TA$ is symmetric, and usually smaller than $A$. $\endgroup$ – arts Dec 2 '17 at 14:19

I must caution you against explicitly forming the product $A^TA$. This process increases the condition number of the active matrix, and it is likely to lower the quality of the computed Cholesky factor.

A far safer alternative is to compute a $QR$ factorization of the matrix $A$, i.e. $A = QR$. Then $$A^T A = R^T Q^T Q R = R^T R$$ reveals that the Cholesky factorization of $A^T A$ is $A^T A = LL^T$, where $L = R^T$ is lower triangular.

In this context it is irrelevant if $A$ is a square matrix or a tall matrix.

I propose that you test both strategies using, say, a scaled copy of the one dimensional discrete Laplacian as your matrix $A$.

Up to a scaling, this matrix is tridiagonal, symmetric and Toeplitz matrix with $2$ on the main diagonal and $-1$ on the subdiagonal. You can compute the exact product $A^TA$, as well as the exact right hand side $b = A^T A x$ corresponding to a solution consisting entirely of, say, small integers. This construction ensures that all rounding errors are attributable to the factorization and triangular solves.

Compare the computed solution $\hat{x}$ using either strategy to your chosen solution $x$. The condition number of $A$ is $O(n^2)$, so the condition number of $A^T A$ should defeat single precision calculations around $n=100$.

In MATLAB you can construct dense Toeplitz matrices using the "toeplitz" command. You force MATLAB to work in single precision using the "single" qualifier.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.