@property<EOL><INDENT>def woodbury_chol(self):<DEDENT>
if self._woodbury_chol is None:<EOL><INDENT>if self._woodbury_inv is not None:<EOL><INDENT>winv = np.atleast_3d(self._woodbury_inv)<EOL>self._woodbury_chol = np.zeros(winv.shape)<EOL>for p in range(winv.shape[-<NUM_LIT:1>]):<EOL><INDENT>self._woodbury_chol[:, :, p] = pdinv(winv[:, :, p])[<NUM_LIT:2>]<EOL><DEDENT><DEDENT>elif self._covariance is not None:<EOL><INDENT>raise NotImplementedError("<STR_LIT>")<EOL>B = self._K - self._covariance<EOL>tmp, _ = dpotrs(self.K_chol, B)<EOL>self._woodbury_inv, _ = dpotrs(self.K_chol, tmp.T)<EOL>_, _, self._woodbury_chol, _ = pdinv(self._woodbury_inv)<EOL><DEDENT>else:<EOL><INDENT>raise ValueError("<STR_LIT>")<EOL><DEDENT><DEDENT>return self._woodbury_chol<EOL>
return $L_{W}$ where L is the lower triangular Cholesky decomposition of the Woodbury matrix
$$
L_{W}L_{W}^{\top} = W^{-1}
W^{-1} := \texttt{Woodbury inv}
$$