posted with カエレバ

伊藤真翔泳社 2018-01-24

Amazonで検索

楽天市場で検索

Yahooショッピングで検索

はじめに

以前、有名な最適化問題として

線形計画法や

myenigma.hatenablog.com

二次計画法

myenigma.hatenablog.com

ポートフォリオ最適化

myenigma.hatenablog.com

ナップサック問題

myenigma.hatenablog.com

などの説明と簡単なサンプルコードを紹介しましたが、

今回の記事では、Lasso正則化, Ridge正則化と呼ばれる

統計学や機械学習でよく使用される最適化問題の概要と、

このLasso正則化, Ridge正則化を解くサンプルコードについて紹介したいと思います。

Lasso正則化, Ridge正則化とは？

正則化とは、

データからモデルを推定する回帰分析において、

過学習を防ぐためのモデル回帰法の一つです。

www.analyticsvidhya.com

例えば、

下記のような普通の最小二乗法では

f:id:meison_amsl:20170904013935p:plain

すべてのデータにおいて、

誤差の二乗を減らすだけなので、

過学習が起きやすいですが、

下記のように、モデルパラメータをできるだけ小さくする

という項を追加することで、

f:id:meison_amsl:20170904014308p:plain

過学習が防げる場合があります。

これはモデルパラメータがモデルの複雑さを表している場合、

モデルの複雑さを減らすように回帰が行われるようになるためです。

上記の式ではモデルパラメータの絶対値の和を最小化していますが、

この正則化の方法をLasso正則化と呼び、

モデルパラメータの二乗値の和を最小化する場合は、

Ridge正則化と呼びます。

通常の最小二乗法は、

解析的に解を得ることができますが、

前述のLasso正則化とRidge正則化の最適化式は、

解析的に解を得ることができないので、

最適化手法を使って解く必要があります。

myenigma.hatenablog.com

JuliaによるLasso正則化とRidge正則化最適化サンプルコード

下記が、Juliaの最適化ライブラリであるJuMPを使った

LassoとRidge正則化サンプルコードです。

myenigma.hatenablog.com

下記のコードでは、sin関数をベースにした元データを

5次関数で最小二乗法フィッティングしています。

通常の最小二乗法と、

Lasso正則化、Ridge正則化の3つ手法でフィッティングした結果を比較しました。

また今回の最適化では非線形最適化ライブラリであるIpoptを使いました。

using JuMP
using Ipopt
using PyCall
solver = IpoptSolver()

@pyimport matplotlib.pyplot as plt


function main1()
    println("JuMP lasso sample1")

    x = [i for i in 0:0.2:5]
    y = sin.(x) + (rand(length(x)) - 0.5) * 1.5
    ty = sin.(x)

    plt.plot(x,ty, "-k", label="True")
    plt.plot(x,y, ".r", label="raw")
    plt.plot(x,normal_least_square(x, y), "-b", label="normal_least_square")
    plt.plot(x,lasso_least_square(x, y), "-g", label="lasso")
    plt.plot(x,ridge_least_square(x, y), "-c", label="ridge")
    plt.legend()
    plt.show()
end

function lasso_least_square(x, y)
    nd = 5
    lamda = 5.0
    model = Model(solver=solver)
    @variable(model, z[i=1:nd])
    @variable(model, ys[i=1:length(x)])
    @NLobjective(model, Min, sum( (ys[i] - y[i])^2 for i = 1:length(x) ) + sum( z[i]^2 * lamda for i = 1:nd ))
    for i = 1:length(x)
        @NLconstraint(model,  ys[i] == z[1] * x[i]^4 + z[2] * x[i]^3 + z[3] * x[i]^2 + z[4]* x[i] + z[5] )
    end
    status = solve(model)
    return getvalue(ys)
end

function ridge_least_square(x, y)
    nd = 5
    lamda = 2.0
    model = Model(solver=solver)
    @variable(model, z[i=1:nd])
    @variable(model, ys[i=1:length(x)])
    @NLobjective(model, Min, sum( (ys[i] - y[i])^2 for i = 1:length(x) ) + sum( abs(z[i])*lamda for i = 1:nd ))
    for i = 1:length(x)
        @NLconstraint(model,  ys[i] == z[1] * x[i]^4 + z[2] * x[i]^3 + z[3] * x[i]^2 + z[4]* x[i] + z[5] )
    end
    status = solve(model)
    return getvalue(ys)
end


function normal_least_square(x, y)
    nd = 5
    model = Model(solver=solver)
    @variable(model, z[i=1:nd])
    @variable(model, ys[i=1:length(x)])
    @NLobjective(model, Min, sum( (ys[i] - y[i])^2 for i = 1:length(x) ))
    for i = 1:length(x)
        @NLconstraint(model,  ys[i] == z[1] * x[i]^4 + z[2] * x[i]^3 + z[3] * x[i]^2 + z[4]* x[i] + z[5] )
    end
    status = solve(model)
    return getvalue(ys)
end

main1()