Heya Prof. Johnson! 
In most circumstances, you should try to avoid computing the Gram matrix A^TA at all, as it squares the condition number of A.
Absolutely, this is what I’m trying to avoid, but at the moment there doesn’t seem to be another way of approaching this. I am solving a proximal operator for a quadratic (and other similar cases) for ADMM, so it requires solving a ‘tilted’ least squares problem (i.e., a sum-of-squares + a linear function), which would require a factorization of the Gram matrix. Essentially, what’s an easy way of computing the solution to
for several b and c? (Where x here is the optimization variable, of course.)
As mentioned by @jichan, one possibility is to backsolve with R and R^T (where QR = A), but this is also not great for computational performance as back-solving is the time consuming step in each iteration.
(As a side note: I’m also not worried about losing much precision due to a large condition number, especially for sparse-direct solving, since I’ll be using ADMM to solve the problem to a rather moderate precision.)