= Solution
For the first $i$ observations, write
$$
A_i=\Sigma^{-1}+X_{1:i}^TX_{1:i},
\qquad b_i=X_{1:i}^TY_{1:i}.
$$
By <Gaussian conjugacy for a normal linear model>, the prefix posterior is
$$
\beta\mid Y_{1:i}\sim N(A_i^{-1}b_i,A_i^{-1}).
$$
Compute a <Cholesky decomposition> of $A_0=\Sigma^{-1}$ once. If $x_i^T$ is row $i$ of the <design matrix>, then
$$
A_i=A_{i-1}+x_ix_i^T,
\qquad b_i=b_{i-1}+x_iY_i.
$$
The <Rank-one Cholesky update> obtains a triangular factor $A_i=L_iL_i^T$ from $L_{i-1}$ in $O(p^2)$ operations. Two <triangular linear system>[triangular solves] give $m_i=A_i^{-1}b_i$, and, for $z_i\sim N(0,I_p)$, another solve gives
$$
\beta_i=m_i+L_i^{-T}z_i\sim N(m_i,A_i^{-1}).
$$
The initial factorization costs $O(p^3)$ and all $n$ updates and samples cost $O(np^2)$. This is within the requested $O(p^3+np^2+n^2p)$ bound.
Back to article page
