Considering design matrix $X \in \mathbb{R}^{n\times p}$ $(n>p)$ and response $y\in \mathbb{R}^{n}$. The sandwich estimator can be calculated directly using $$(X^TX)^{-1}X^T diag(r^2) X (X^TX)^{-1}$$ where $r$ is the residuals. The sandwich estimator then is the the diagonal of the result. To avoid calculating the inverse, we can also get the sandwich estimator by first projecting $diag(r)$ onto the column space of $X$ and then take the diagonal of the product of the transpose of the coefficients and the coefficients (Note that the coefficients is a ${p\times n}$ matrix).
Since we have already had the QR decomposition when solving the linear system, i.e. $$Q^T X \beta = \begin{pmatrix} R \\ 0 \end{pmatrix}$$ $$Q^T y = \begin{pmatrix} \tilde{y_1} \\ \tilde{y_2} \end{pmatrix}$$ where $R$ is a $p\times p$ upper triangular matrix, then we can reuse $R$ and get $Q^T diag(r)$ and back substitute to get the coefficients.
Observe that the residual is actually $$r = Q \begin{pmatrix} 0 \\ \tilde{y_2} \end{pmatrix}.$$ Can we write $Q^T diag\left(Q \begin{pmatrix} 0 \\ \tilde{y_2} \end{pmatrix}\right)$ as some form so that we only need minimal number of operations on $\tilde{y_2}$ to get $Q^T diag(r)$? In general, is there any other procedure to utilize the QR decomposition efficiently to get the sandwich estimator?
Edit: $\tilde{y_1}$ is the first $d$ elements and $\tilde{y_2}$ is the last $(n-d)$ elements of $Q^T y$.
