Gaussian process likelihood function

I'm trying to understand the likelihood function in Gaussian Process. The book by Rasmussen et al. defined Gaussian Process lml as

$$log~p(y|X) = -\frac{1}{2}y^T\alpha-\sum log L_{ii} - \frac{N}{2}log (2\pi)$$

Where $$\alpha$$ is computed from the lower triangular matrix of the cholesky decomposition (I'm omitting the noise for simplicity): $$L = cholesky(K)$$ $$\alpha = L^T\backslash L\backslash y$$

The implementation in scikit learn is:

L = cholesky(K, lower=True)
alpha = cho_solve((L, True), y_train)
log_likelihood_dims = -0.5 * np.einsum("ik,ik->k", y_train, alpha)
log_likelihood_dims -= np.log(np.diag(L)).sum()
log_likelihood_dims -= K.shape / 2 * np.log(2 * np.pi)
log_likelihood = log_likelihood_dims.sum(-1)

If we compare the $$\alpha$$ given by the equation above and sklearn implementation, they are not same. sklearn will give $$\alpha = L^{-1}y$$ while the equation says $$\alpha=L^T/(L^{-1}y)$$. My hunch is both of them are equivalent expression as the numerical value of lml is not something meaningful, i.e., we are only interested in the hyperparameters that optimizes lml. So, instead of computing the true lml, we are computing a surrogate lml?

On a similar note, I saw several different likelihood in the GPML code written in MATLAB. So, does that mean we can choose any suitable likelihood function (or cost function) to optimize hyperparameter depending on the problem description?

Actually, this code is computing $$\alpha=L^T\backslash L\backslash y$$. With cho_solve you are solving the original system taking advantage of the cholesky descomposition.