Yes, sure, it is possible to do it with two factorizations; one which initially solves the underdetermined system (e.g., via a Moore-Penrose Pseudoinverse or other least-norm type solution) and one which solves the resulting least squares problem, but just computing the two factorizations would override any performance gains from simply factorizing the Gram matrix and reusing this factorization (even if its numerical accuracy is lower, which doesn’t affect ADMM much due to its generally low-to-modest precision). For context, factorizing the matrix takes on the other of the same time as the entire algorithm takes to converge with the known factorization.
Which brings me back to the original question: as per, e.g., this paper, section 3.2, CHOLMOD should support a method for computing the Cholesky factorization of the Gram matrix of A without needing to directly compute the Gram matrix. Is there any way of accessing this method from Julia?
EDIT: Ah, I see (the post became a little clearer after the edit). Even in this case, though, the dimension of the problem is very close to square, and this method would still require two solutions of the problem, which would also place a large constant in the runtime of the problem (this is similar to reusing the R matrix of the QR factorization of A).