Computing the (simplest) sandwich estimator incrementally

The biglm package in R does {incremental, online, streaming} linear regression for data potentially larger than memory. This isn’t rocket science: accumulating $X^{T}X$ and $X^{T}Y$ is trivial; the package just goes one step better than this by using Alan Miller’s incremental $QR$ decomposition code to reduce rounding error in ill-conditioned problems. 

The code also computes the Huber/White heteroscedasticity-consistent variance estimator (sandwich estimator). Someone wants a reference for this. There isn’t one, because it’s too minor to publish, and I didn’t have a blog ten years ago.  But I do now. So:

The Huber/White variance estimator $A^{-1}B{A}^{-1}$, where $A^{-1}=(X^{T}X{)}^{-1}$ and $B={\left(X^T\right(Y-\hat{\mu }\left)\right)}^{\otimes 2}$

The $(i,j)$ element of $B$ is
$$\sum _{k=1}^N{x}_{ki}(y_k-x_k\hat{\beta })x_{kj}(y_k-x_k\hat{\beta })$$

Multiplying this out, we get
$$\sum _{k=1}^N{x}_{ki}{x}_{kj}{y}_k^2$$
and about $2p$ terms that look like
$$\sum _{k=1}^N{x}_{ki}{x}_{kj}{x}_{k\ell }{y}_k{\hat{\beta }}_{\ell }$$
and about $p^2$ terms that look like
$$\sum _{k=1}^N{x}_{ki}{x}_{kj}{x}_{k\ell }{x}_{km}{\hat{\beta }}_{\ell }{\hat{\beta }}_m$$

We can move the $\beta$s outside the sums, so the second sort of terms look like $${\hat{\beta }}_{\ell }\left(\sum _{k=1}^N{x}_{ki}{x}_{kj}{x}_{k\ell }{y}_k\right)$$ and the third sort look like $${\hat{\beta }}_{\ell }\left(\sum _{k=1}^N{x}_{ki}{x}_{kj}{x}_{k\ell }{x}_{km}\right){\hat{\beta }}_m$$

Now if we define $Z$ to have columns $x_i{x}_j$ and $x_{i}y$ (for all $i,j$), the matrix $Z^{T}Z$ contains all the $x$ and $y$ pieces needed for $B$.  The obvious thing to do is just to accumulate $Z^{T}Z$ in R code, one chunk at a time.

If you were too convinced of your own cleverness you might realise that $(X,Z)$ could be fed into the $QR$ decomposition as if it were $X$, and that you’d get $Z^{T}Z$ For Free! Where ‘for free’ means at $O\left(\right(p^2{)}^3)$ extra computing time plus the mental anguish of reconstructing $Z^{T}Z$ from the $QR$ decomposition.  It’s not a big deal, since the computation is dominated by the $O\left(np\right)$ cost of reading the data, but it does look kinda stupid in retrospect.

I suppose that means I’ve learned something in ten years.