Solution (source code)

= Solution

The <least-squares normal equations> are
$$
\boxed{Hx=c,\qquad H=G^TG,\quad c=G^Tb}.
$$
The <full column rank> of $G$ gives $z^THz=\|Gz\|_2^2>0$ for nonzero $z$, so $H$ is a <symmetric positive-definite matrix>. Apply the <Conjugate gradient method> to this <linear system>. To verify that its solution is the desired one, let $x_*$ solve the displayed equations. Then $G^T(Gx_*-b)=0$, and for any $x$ the <Pythagorean identity> gives
$$
\|Gx-b\|_2^2=\|Gx_*-b\|_2^2+\|G(x-x_*)\|_2^2.
$$
Thus $x_*$ is the unique minimizer of the <linear least-squares problem>.

One can apply $H$ without forming the <Gram matrix>: compute $Hp=G^T(Gp)$. More explicitly, initialize $x_0$, $r_0=c-Hx_0$, $p_0=r_0$; until the <linear-system residual> vanishes, use the <Conjugate gradient method> recurrence
$$
\alpha_j=\frac{r_j^Tr_j}{p_j^THp_j},\quad
x_{j+1}=x_j+\alpha_jp_j,\quad
r_{j+1}=r_j-\alpha_jHp_j,
$$
$$
\beta_j=\frac{r_{j+1}^Tr_{j+1}}{r_j^Tr_j},\qquad
p_{j+1}=r_{j+1}+\beta_jp_j.
$$
Each dense application of $H$ this way costs $O(mn)$ operations. The <condition number> satisfies $\kappa_2(H)=\kappa_2(G)^2$, where $\kappa_2(G)=\sigma_{\max}(G)/\sigma_{\min}(G)$; this is the conditioning cost of the <least-squares normal equations>.