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)?
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.