I am coding a simple univariate GP and was planning to update the posterior sequentially every time a new sample is taken. It is well known that for a partitioned Gaussian as follows
\begin{equation} \left(\begin{array}{l} y_{a} \\ y_{b} \end{array}\right) \sim N\left(\left(\begin{array}{l} \mu_{a} \\ \mu_{b} \end{array}\right), \left(\begin{array}{ll} \sum_{a a} & \sum_{a b} \\ \sum_{b a} & \sum_{b b} \end{array}\right)\right) \end{equation}
the conditional distribution is a Gaussian with mean and var given by
\begin{equation} \mu_{a|b}=\mu_a + \Sigma_{ab}\Sigma_{bb}^{-1}(y_b-\mu_b)\\ \Sigma_{a|b}=\Sigma_{aa}-\Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba} \end{equation}
where $\Sigma_{ab}=\Sigma_{ba}^T$ and $\Sigma$ is given by the kernel function of x values.
So in the Bayesian Optimization / GP context, $a$ would be the set of GP function values and $b$ the observed values. So if the values of $b$ are not given all at once, I would like to update the GP sequentially and keep an updated version of the GP at each step. My thinking was that the most computationally efficient approach would be to update it with each variable individually at each time step (as opposed to updating a growing vector of observed values and using this to update the original GP). This would remove the need to calculate $\Sigma_{bb}^{-1}$ for a growing vector.
Instead $\Sigma_{aa}$ would be updated at each iteration, essentially becoming $\Sigma_{(a|b_1)|b_2}$ for the second iteration. However since $\Sigma_{ab}$ and its transpose are calculated via the kernel function I'm confused if this makes sense and how I would even do it since I'm no longer referencing the x values used to create the initial $\Sigma_{aa}$. So do GPs have to be updated with all observed values only once? Then recalculated when a new observation is taken?