MyEnigma

とあるエンジニアのブログです。#Robotics #Programing #C++ #Python #MATLAB #Vim #Mathematics #Book #Movie #Traveling #Mac #iPhone

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

統計学がわかる 【回帰分析・因子分析編】 (ファーストブック)

統計学がわかる 【回帰分析・因子分析編】 (ファーストブック)

目次

はじめに

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

線形計画法や

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

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()

上記のコードは下記のリポジトリでも公開しています。

github.com

 

上記のコードを実行すると下記のような結果が得られます。

f:id:meison_amsl:20170904035411p:plain

上記のグラフを見ると、

通常の最小二乗法では、元データの偏りに影響されて、

元の曲線から遠い結果が得られている部分もありますが、

Lasso正則化とRidge正則化の場合、

元データの偏りに引っ張られることが少なくなっていることがわかります。

 

参考資料

CVXGEN: Code Generation for Convex Optimization

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

統計学がわかる 【回帰分析・因子分析編】 (ファーストブック)

統計学がわかる 【回帰分析・因子分析編】 (ファーストブック)

 

MyEnigma Supporters

もしこの記事が参考になり、

ブログをサポートしたいと思われた方は、

こちらからよろしくお願いします。

myenigma.hatenablog.com