Is there an iterative way to evaluate least squares estimation? Suppose to have a set of data $\{y_i, u_i\}_{i=1}^m$, where $y_i \in \mathbb{R}$ and $u_i \in \mathbb{R}^n$.
The claim is that 
$$y_i = u_i^\top \theta + \varepsilon$$
where $\theta \in \mathbb{R}^n$ is a vector of unknown parameters and $\varepsilon$ is some kind of zero-average noise (e.g. Gaussian noise).
If I build the followings:
$$Y = \left[\begin{array}{c}y_1\\y_2\\\vdots\\y_m\end{array}\right],
U = \left[\begin{array}{c}u_1^\top\\u_2^\top\\\vdots\\u_m^\top\end{array}\right], $$
then I can find a least squares estimation of $\theta$:
$$\hat\theta = (U^\top U)^{-1}U^\top Y.$$
Anyway, suppose that $n$ and $m$ are very huge. Then solving the previous problem numerically is very hard due to memory consumption.
My question is: is there some iterative methods such that I can find $\hat\theta$ iteratively in order to save memory?
I hope that there exists some schema like the following:
$$\begin{align}
\hat\theta_1 & = f(y_1, u_1)\\
\hat\theta_2 & = g(\theta_1, y_2, u_2) \\
\hat\theta_3 & = g(\theta_2, y_3, u_3) \\
\vdots\\
\hat\theta_m & = g(\theta_{m-1}, y_m, u_m) = \hat\theta =  (U^\top U)^{-1}U^\top Y.\\
\end{align}$$

I recommend checking out the paper "LSMR: An iterative algorithm for sparse least squares problems" by Fong and Saunders.
This algorithm is included in scipy.
