Motivation

This PR resolves the following issue:

• GPSampler suddenly slows down once n_trials hits exactly 100 #6230

See also the issue in NumPy:

• BUG: Linear algebra heavy calculations experience sudden slowdown after the size of 100 numpy/numpy#29884

Basically, np.linalg.inv suddenly slows down after cov_Y_Y.shape exceeds (100, 100).
We avoid this issue by using np.linalg.cholesky.

Here are some backgrounds.
Let's denote $C \in \mathbb{R}^{N \times N}$ as the kernel matrix with the addition of a noise at the diagonal elements and $k = [k(x^\star, x_1), \dots, k(x^\star, x_N)]\in \mathbb{N}$ as the kernel vector at a new point $x^\star$.
The Gaussian process infers the posterior variance using $k^\top C^{-1} k$, necessitating the matrix inversion in our current implementation.
However, as $C$ is a positive definite matrix owing to the kernel matrix nature, we can indeed rewrite $v = C^{-1} k$ as $k = C v = L L^\top v$ where $L \in \mathbb{R}^{N \times N}$ is the lower triangular matrix yielded by the Cholesky decomposition of $C$.
We first obtain $u = L^\top v$ by solving $L u = k$ and then obtain $v$ by solving $L^\top v = u$.
Please note that each linear system can be solved with the time complexity of $O(N^2)$ because $L$ is a triangular matrix.

By doing so, we can avoid the matrix inversion, leading to a significant speedup.
The by-product of this modification is the numerical stability.
The speedup effect can be found below:

Description of the changes

• Replace np.linalg.inv with np.linalg.cholesky