Julia言語でカーネル回帰

カーネル回帰のjulia実装が研究で必要になったから１０分で実装した．

とりあえず，前に書いたカーネルを呼ぶ式．

function rbf_inners(X,Y,c)
    rbf_inner(x,y,c) = exp( -c*norm(x-y)^2 )
    d, n = size(X)
    d, m = size(Y)
    K = zeros(m,n)
    for j = 1:n
        for i = 1:m
           K[i,j] = rbf_inner(Y[:,i],X[:,j],c)
        end
    end
    return K
end

function poly_inners(X,Y,c,deg)
    poly_inner(x,y,c,deg) = (c+x'*y)^deg
    d, n = size(X)
    d, m = size(Y)
    K = zeros(m,n)
    for j = 1:n
        for i = 1:m
           K[i,j] = poly_inner(Y[:,i],X[:,j],c,deg)
        end
    end
    return K
end

function lapras_inners(X,Y,c)
    lapras_inner(x,y,c) = exp(-c*sum( abs.(x .- y)))
    d, n = size(X)
    d, m = size(Y)
    K = zeros(m,n)
    for j = 1:n
        for i = 1:m
           K[i,j] = lapras_inner(Y[:,i],X[:,j],c)
        end
    end
    return K
end

function Kernel(X, Y; c=1.0, deg=3.0, kernel="rbf")
    if kernel == "rbf"
        K = rbf_inners(X,Y,c)
    elseif kernel == "poly"
        K = poly_inners(X,Y,c,deg)
    elseif kernel == "linear"
        K = poly_inners(X,Y,0,1)
    elseif kernel == "lapras"
        K = lapras_inners(X,Y,c)
    else
        error("Kernel $kernel is not defined")
    end
    return K
end

カーネル回帰の本体はこんな感じ．正則化のパラメータがlambda．これが大きいと正則化が強くなる． cとdegがカーネルのパラメータ．cはRBFカーネルの定数で，dは多項式カーネルの次数． Lは訓練データの数．逆行列を求めるなら，\の方がいいってどこかで誰かが言っていた気がするので，invは使ってない． 機械学習民なら，カーネル回帰の式は特に何の説明もいらないと思うがわかりやすい解説はQiitaにあると思う．（これとか）

function Kreg(X_train, Y_train, X_test; lambda=1, c=1, deg=3.0, kernel="rbr")
    K = Kernel(X_train, X_train)
    L = size(K)[1]
    alpha = (K + lambda * Matrix(1.0I, L, L)) \ Y_train'
    k = Kernel(X_train, X_test)
    Y_pred = k * alpha
    return Y_pred
end

適当にデータを作ってプロットしてみる．

using LinearAlgebra
using Distributions
function make_1d_data(L)
    f(x) = x/2 + sin(x)*cos(x)+1
    dist = Uniform(-6, +6)
    x = rand(dist,(1, L))

    eps = 0.1
    noise_dist = Uniform(-eps, +eps)

    y = zeros(1, L)
    for i = 1:L
        noise = rand(noise_dist)
        y[1,i] = f(x[i]) + noise
    end
    return x, y
end

function main()
    L = 150
    N = 50
    X_train, Y_train = make_1d_data(L)
    X_test, Y_test = make_1d_data(N)
    Y_pred = Kreg(X_train, Y_train, X_test)
    L2error = norm(Y_test[:] .- Y_pred[:])^2 / size(Y_test)[2]
    println("Ferror:$L2error")
    plots(X_train, Y_train, X_test, Y_pred)
end

あっているように思います．ハイパラcとlambdaを調節したらもっと汎化すると思う． ちなみにプロットに使ったのはこちらのコード

using Plots
using Plots.PlotMeasures
function plots(x, y, x_test, y_pred)
    fnt = font(20, "Helvetica")
    plt = plot(legend=:topleft, legendfont=fnt,
                bottom_margin = 10mm,
                left_margin = 10mm,
                size=(900,500),
                xlabel ="x", ylabel="y",
                xguidefont = fnt,
                yguidefont = fnt,
                xtickfont = fnt,
                ytickfont =fnt
              )
    plot!(plt, x[:], y[:], label="train", color="blue", st=:scatter)
    plot!(plt, x_test[:], y_pred[:], label="test", color="red", st=:scatter)
    png("tmp.png")
    println("saved")
end

追記

このカーネル回帰の実装は，重回帰にしてもちゃんと動く． ためしに，下の人工データを使ってみるとちゃんとコスト関数が小さくなってる．

function make_2d_data(L)
    #f(x1,x2) = sin(x1+x2^2)
    f(x1,x2) = (x1^2+x2^2)*exp(1-x1^2-x2^2)
    dist = Uniform(-2, +2)
    x = rand(dist,(2, L))

    y = zeros(1, L)
    for i = 1:L
        y[1,i] = f(x[1,i],x[2,i])
    end
    return x, y
end

Ferror:0.008802513523177683