How can I numerically solve a least square problem with L2 penalty

Question

I want to solve a least square problem that: \begin{equation*} \text{minimize } \frac{1}{2}\|Ax-b\|_2^2 + \frac{1}{2}\lambda\|x\|_2^2 \end{equation*}

From the normal equation, I have to solve the normal equation that: \begin{equation*} (A^TA + \lambda I)x = A^Tb\end{equation*}

But in my problem, $A$ is a very high dimensional matrix ($A \in R^{512^2 \times 512^2}$), so I can not use method like QR decomposition. How can I solve this problem?

If helps, $A$ is a convolution operator in image processing, but it can be represented in a matrix form.