Cholesky Decomposition and simple (!) gaussian process regression

cholesky.jl

function chol_(A::Matrix, L::Matrix)
    for i = 1:size(A,1)
        s = 0
        @simd for j = 1:i-1
            @inbounds L[i,j] = (A[i,j] - sum([L[i,k]*L[j,k] for k = 1:j-1])) / L[j,j]
            @inbounds s += L[i,j]^2
        end
        L[i,i] = sqrt(A[i,i] - s)
    end
    L
end

cholesky(A) = chol_(A, zeros(A))
function cholesky!(A)
    chol_(A, A)
    for i = 1:size(A,1)
        @simd for j = 1:i-1 @inbounds A[j,i] = 0 end
    end
    A
end

gp.jl

function learn_GP(X::Matrix, y::Vector, K::Matrix, noise::Real)
    L = cholesky!(K + noise * eye(size(K,1)))
    alpha = L' \ (L \ y)
    alpha, (x::Vector -> [kernel(X[i,:], x) for i in size(X,1)])
end

function learn_GP(X::Matrix, y::Vector, kernel::Function, noise::Real)
    K = zeros(size(X,1), size(X,1))
    for i = 1:size(X,1), j = 1:i
        @inbounds K[i,j] = kernel(X[i,:], X[j,:])
        @inbounds K[j,i] = K[i,j]
    end
    learn_GP(X, y, K, noise)
end

function predict_GP(alpha::Vector, kernel::Function, x::Vector)
    s = 0
    @simd for i = 1:size(alpha,1)
        @inbounds s += alpha[i] * kernel(x)
    end
    s
end