Swiftでガウス過程

ガウス過程回帰をSwiftで実装します。

ガウス過程は線形回帰の重みパラメータを積分消去し、カーネル関数（基底関数の内積）のみで関数を表現するモデルです。
確率モデルなので、予測に対する分散も表現できます。

Pythonで実装したものがあるので、これを参考にしてSwiftでも実行できるようにしたいと思います。

Calc()は前回の記事で作成した行列計算をまとめたクラスです。
今回は負の対数尤度、対数尤度をガウスカーネルの各ハイパーパラメータで微分したものまで書きました。

import Foundation

class GaussianProcess {
    let calc = Calc()
    let X: [[Double]]
    var theta: [Double]
    var y: [Double]
    
    
    init(X: [[Double]]){
        self.X = X
        self.theta = [Double.random(in: 0 ..< 1), Double.random(in: 0 ..< 1), Double.random(in: 0 ..< 1)]
        self.y = [Double](repeating: 0, count: X.count)
    }
    
    func kernel(X_i: [Double], X_j: [Double]) -> Double{
        let r: Double = calc.vsum(vector: calc.vpow(vector: calc.vmv(vector: X_i, vector_: X_j), power: 2))
        let k: Double = exp(self.theta[0]) * exp(-r / exp(self.theta[1]))
        return k
    }
    
    func Kernel_matrix() -> [[Double]] {
        let N: Int = self.X.count
        var K = [[Double]](repeating: [Double](repeating: 0, count: N), count: N)
        for i in 0..<N {
            for j in 0..<N {
                K[i][j] = self.kernel(X_i: self.X[i], X_j: self.X[j])
            }
        }
        K = calc.mpm(matrix: K, matrix_: calc.sxm(scalar: exp(self.theta[2]), matrix: calc.eye(n: N)))
        return K
    }
    
    func partial_kernel(X_i: [Double], X_j: [Double]) -> [Double] {
        let r: Double = calc.vsum(vector: calc.vpow(vector: calc.vmv(vector: X_i, vector_: X_j), power: 2))
        var partial_k = [Double](repeating: 0, count: self.theta.count)
        partial_k[0] = exp(theta[0]) * exp(-r / exp(self.theta[1]))
        partial_k[1] = exp(theta[0]) * exp(-r / exp(self.theta[1])) * exp(self.theta[1]) * r
        partial_k[2] = exp(theta[2])
        return partial_k
    }
    
    func partial_Kernel_matrix() -> [[[Double]]] {
        let N: Int = self.X.count
        var partial_K = [[[Double]]](repeating: [[Double]](repeating: [Double](repeating: 0, count: N), count: N), count: N)
        for i in 0..<N {
            for j in 0..<N {
                let partial_k = self.partial_kernel(X_i: self.X[i], X_j: self.X[j])
                for k in 0..<self.theta.count {
                    partial_K[k][i][j] = partial_k[k]
                }
            }
        }
        return partial_K
    }

    func log_likelihood(K: [[Double]], K_inv: [[Double]], det_K: Double) -> Double {
        let N: Double = Double(K.count)
        let L: Double = -(N/2)*log(2*Double.pi) - (1/2)*log(det_K) - (1/2)*calc.mdot(matrix: calc.mdot(matrix: calc.transpose(matrix: [self.y]), _matrix: K_inv), _matrix: [self.y])[0][0]
        return -L //負の対数尤度
    }
    
    func partial_L(K_inv: [[Double]]) -> [Double] {
        let dK: [[[Double]]] = self.partial_Kernel_matrix()
        var partial_L = [Double](repeating: 0, count: self.theta.count)
        let K_inv_dot_y: [Double] = calc.mdot(matrix: K_inv, _matrix: [self.y])[0]
        for i in 0..<self.theta.count {
            partial_L[i] = calc.trace(matrix: calc.mdot(matrix: K_inv , _matrix: dK[i])) - calc.mdot(matrix: calc.mdot(matrix: calc.transpose(matrix: [K_inv_dot_y]), _matrix: dK[i]), _matrix: [K_inv_dot_y])[0][0]
        }
        return partial_L
    }
}