Could you add an inline comment about this design choice?
Namely, we could use mean (posterior mean), median, min, or any kind of operators, right?
Plus, as you may already know, we often integrate out the value because we can calculate the posterior distributions for each running trial.

Thank you for the comment. I have added a few inline comments. I also left some notes on issue #6392 regarding the future implementation plan.

I think we can compute this more efficiently by reusing self._cov_Y_Y_chol, which we already have.
This would reduce the computational cost from the current
　　 $$O((\text{n\_train} + \text{n\_running})^3)$$
　　 $$= O(\text{n\_train}^3 + 3\,\text{n\_train}^2\:\text{n\_running} + 3\,\text{n\_train}\:\text{n\_running}^2 + \text{n\_running}^3),$$
to
　　 $$O(\text{n\_train}^2\:\text{n\_running} + \text{n\_train}\:\text{n\_running}^2 + \text{n\_running}^3).$$
Considering the fact that the $\text{n\_train}$ is the dominant size ($\text{n\_train} >> \text{n\_running}$), the reducing the order of $\text{n\_train}$ from 3 to 2 will improve runtime.
This will also eliminate the need for storing self._cov_Y_Y.
However, it is also acceptable to leave the current one and address this in a follow-up PR, since the current implementation itself is not incorrect.

We can calculate this as follows:

Let $A$ be the symmetric positive definite matrix (in this case, self._cov_Y_Y) and define its Cholesky decomposition as $LL^\top=A$ (in this case, $L$ corresponds to self._cov_Y_Y_chol). Let $M$ be a symmetric positive definite matrix partitioned as:
　　 $M$ = [ [ $A\:$ $X$ ] [ $X^\top\:$ $Y$ ] ]
(in this case, $X$ corresponds to self.kernel(X_running), and $Y$ corresponds to self.kernel(X_running, X_running)).
Then, the Cholesky decomposition of $M$ is provided as follows:
　　 $\mathrm{chol}(M)$ = [ [ $L\:$ $0$ ] [ $Z^\top\:$ $C$ ] ],
where $Z := L^{-1} X$ and $C := \mathrm{chol}(Y - Z^\top Z)$.
Note that we can compute $Z$ efficiently by solving $LZ = X$, since $L$ is the triangular matrix. (Use scipy.linalg.solve_triangular.)

proof.
Assume that $M$ admits a block Cholesky factorization of the form:
　　 $M$ := [ [ $L\:$ $0$ ] [ $B\:$ $C$ ] ] [ [ $L^\top\:$ $B^\top$ ] [ $0\:$ $C^\top$ ] ]
　　 = [ [ $LL^\top\:$ $LB^\top$ ] [ $BL^\top\:$ $BB^\top + CC^\top$ ] ].
Comparing this to $M$ = [ [ $A\:$ $X$ ] [ $X^\top\:$ $Y$ ] ], we obtain:
　　 $LB^\top = X$ and $BB^\top + CC^\top = Y$. 　　(1)
Thus, it follows that $B^\top = L^{-1}X$. Defining $Z := L^{-1}X$, we have $B = Z^\top$, and hence $BB^\top = Z^\top Z$.
Substituting this into the Equation (1) yields:
　　 $CC^\top = Y - Z^\top Z$
　　 $= Y - (L^{-1}X)^\top L^{-1}X$
　　 $= Y - X^\top (L^{-1})^\top L^{-1}X$
　　 $= Y - X A^{-1} X$.
Since $M$ is symmetric positive definite, the Schur complement $Y - X A^{-1} X = Y - Z^\top Z$ is also symmetric positive definite. Therefore, its Cholesky factor exists and $C$ can be chosen as $C = \mathrm{chol}(Y - Z^\top Z)$.
This completes the proof. □

Thank you so much for the precise explanation! I made the update accordingly.

I think we can also compute here little bit more efficiently by reusing self._cov_Y_Y_inv_Y, leading $O( (\text{n\_train} + \text{n\_running})^2)$ to $O(\text{n\_train}~\text{n\_running} + \text{n\_running}^2)$.
However, this is not very important and we can leave here as it is, since the overall complexity of append_running_data is bottlenecked by computing cov_Y_Y_chol.

Thank you for the helpful suggestion! I agree that the forward substitution, scipy.linalg.solve_triangular(cov_Y_Y_chol, self._y_all.cpu().numpy(), lower=True) can likely be reduced to O(n_train * n_running + n_running**2) by caching the corresponding train-only term,
scipy.linalg.solve_triangular(cov_Y_Y_chol, self._y_train.cpu().numpy(), lower=True) in _cache_matrix.

What I was less certain about was whether the backward substitution could also be reduced to the same complexity. As you mentioned, I think this would be a good optimization to consider in a future improvement.

Thank you for the reply. Indeed, yes. We need $W := AX^{-1}$ by solving $L^\top W = Z$, which takes $O(\text{n\_train}^2 + \text{n\_running})$. Also, the overall complexity of append_running_data is already bottlenecked by the computation of cov_Y_Y_chol, which we cannot improve more.
Therefore, we can achieve only a constant-factor improvement here, which is the reason why I felt this is not very important in contrast to the one I mentioned here. Also, I'm not very certain how this constant-factor improvement contribute to the overall runtime without conducting benchmarking on it.