Juliaでナップサック問題を解いてみた

前回↓記事でJulia入門して、もうちょっとアルゴリズムを書いてみたい！となりましたのでそういえば使ったことない遺伝的アルゴリズムをやってみました。

とりあえず簡単なもので試したいので、数理最適化の初歩であるナップサック問題を扱いました。とりあえず動きそうなものが出来たので良かったです。

Githubはこちら

ナップサック問題

ナップサック問題とは、以下のような組み合わせ最適化問題のことです。

n 種類の品物（各々、価値 vi、重量 wi）が与えられたとき、重量の合計が W を超えない範囲で品物のいくつかをナップサックに入れて、その入れた品物の価値の合計を最大化するには入れる品物の組み合わせをどのように選べばよいか

ウィキペディアより引用

この手の問題はNP困難です。ようするにパラメータが増えることで計算量が膨大に膨れ上がり実質的に解くことができない問題なので、近似解を効率良く得る工夫が求められます。

このナップサック問題に関しては全探索法以外での解法としては効率の良い貪欲法が見つかっています。が、今回はこれを遺伝的アルゴリズムで解いてみます。

遺伝的アルゴリズムとは

遺伝的アルゴリズムは生物の進化の仕組みを模したアルゴリズムです。有性生物の配合とその選択によって最も優秀な遺伝子をもつ個体を探索する方法です。アルゴリズムは以下の順番で動いていきます。

1. 遺伝子を生成する

2. 最も良い遺伝子を評価する

3. 良い成績を残す遺伝子を選別する

4. 残った遺伝子を配合する

5. 配合した遺伝子を突然変異させ、微小な変化を加える

6. 2に戻り、遺伝子をどんどん先鋭化させていく（世代交代）

遺伝子の作り方や評価方法を工夫することで様々な問題を解くことが可能です。組み合わせ最適化に限らず他の数理最適化問題で良く使われます。

一般的なパラメータ推定でよく出てくる局所的最適化手法とは異なる収束条件を持ち、問題によって得意不得意があります。個人的な考えとしては以下のような違いがあると思います。

• 局所的最適化

  • 近傍付近では最も良いパラメータを選択しやすい。ただし、正解とは全く誤ったパラメータに収束する可能性がある

• 遺伝的アルゴリズム

  • 必ず最適なパラメータを選択することはないが、正解付近には近づきやすい

ようするに、遺伝的アルゴリズムは安定的に70~80点を取るのが得意なアルゴリズムです。

Juliaによる実装

初期条件

まず、シミュレーションで使うパラメータを入力します。

• generation_iter: 世代交代の回数

• cost_limit: ナップサックのコストの上限

• parent_size: 一世代の親集団の個数

• children_size: 一世代で選択される遺伝子の個数

• item_size: ナップサックに入れる品物の候補数

• mutation_prob: 突然変異を起こす確率

• knapsack: ナップサックの設定

Base.@kwdef struct KnapsackGAParameters
    seed_number::Union{Int, Nothing} = nothing

    generation_iter::Int = 50
    cost_limit::Int = 50
    parent_size::Int = 100
    children_size::Int = 10
    item_size::Int = 20
    mutation_prob::Float16 = 0.5

    knapsack::Dict{String, NTuple{20, Int}} = Dict(
        "cost" => (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20),
        "value" => (2, 3, 5, 5, 6, 9, 10, 12, 15, 15, 16, 16, 17, 19, 20, 22, 24, 28, 29, 32),
    )
end

遺伝子

ナップサック問題での遺伝子は各品物を入れる or 入れないをバイナリで表現したバイナリコーディングを使います。即ち、全て0であれば品物を一つも入れずに、逆に全て1であれば品物を全てナップサックに入れた状態を表します。

mutable struct Gene
    keep_idx::Int

    parents::Array{Int, 2}
    childerens::Array{Int, 2}
    scores::Array{Int, 1}
end

function build_param(params::KnapsackGAParameters)
    return Gene(
        params.parent_size - params.children_size,
        Int.(Random.bitrand((params.parent_size, params.item_size))),
        zeros(Int, params.children_size, params.item_size),
        zeros(Int, params.parent_size),
    )
end

評価

遺伝子の評価は以下の2条件で行います。

• コストがオーバーしない

• 品物の総価値が高い

function evaluation!(params::KnapsackGAParameters, gene::Gene)
    costs::Matrix{Int} = sum(gene.parents' .* params.knapsack["cost"], dims=1)
    values::Matrix{Int} = sum(gene.parents' .* params.knapsack["value"], dims=1)

    gene.scores = [
        if costs[idx] <= params.cost_limit values[idx]
        else params.cost_limit - costs[idx]
        end
        for idx in 1:params.parent_size
    ]

    best_gene_id::Int = argmax(gene.scores)
    return best_gene_id, costs[best_gene_id], values[best_gene_id]
end

選択

現在の遺伝子から成績が良いものを数個選別します。ここでは、トーナメント選択という近年メジャーなものを使いました。

function tournament_selection!(params::KnapsackGAParameters, gene::Gene)
    block::Int = div(params.parent_size, params.children_size)
    for idx in 1:block
        idx_max::Int = argmax(gene.scores[(idx - 1) * block + 1: idx * block])
        gene.childerens[idx, :] = gene.parents[(idx - 1) * block + idx_max, :]
    end
end

交配

選択された遺伝子を配合します。同じ番号の遺伝子の内どちらか一方を受け継ぎ、新しい遺伝子を生成します。

function uniform_crossbreeding!(params::KnapsackGAParameters, gene::Gene)
    gene_code::Int = -1
    parent_idx::Matrix{Int} = sample(1:params.children_size, (gene.keep_idx, 2))

    for p_idx in 1:gene.keep_idx
        for i_idx in 1:params.item_size
            if rand() <= 0.5
                gene_code = gene.childerens[parent_idx[p_idx, 1], i_idx]
            else
                gene_code = gene.childerens[parent_idx[p_idx, 2], i_idx]
            end
            gene.parents[p_idx, i_idx] = gene_code
        end
    end
end

突然変異

一定の確率でランダムに一つの遺伝子を選び、符号を入れ替えます。

function mutation!(params::KnapsackGAParameters, gene::Gene)
    for p_idx in 1:gene.keep_idx
        if rand() <= params.mutation_prob
            mutation_num = rand(1:params.item_size)
            gene.parents[p_idx, mutation_num] = div(gene.parents[p_idx, mutation_num] + 1, 2)
        end
    end
end

実行結果

これらを一つのパイプラインにし、実行します。

using Random
using StatsBase

include("parameters.jl")
include("gene.jl")

function knapsack_ga_process()
    # initiation
    params::KnapsackGAParameters = KnapsackGAParameters()
    # MersenneTwister(params.seed_number)
    gene::Gene = build_param(params)
    best_gene_id::Int = -1

    for idx in 1:params.generation_iter
        # evalution
        best_gene_id, best_cost, best_value = evaluation!(params, gene)

        # selection
        tournament_selection!(params, gene)

        # crossbreeding
        uniform_crossbreeding!(params, gene)

        # mutation
        mutation!(params, gene)

        # insertion eugenic gene
        insert_children!(params, gene)

        println(
            "Generation $(idx): cost: $(best_cost) value: $(best_value)"
        )
    end
    println("Best gene: ", gene.parents[best_gene_id, :])
end

すると、以下のような結果と、その組み合わせが出てきます。（めっちゃ収束早いw、何か間違えていたら教えてください。）

Generation 1: cost: 41 value: 58
Generation 2: cost: 49 value: 75
Generation 3: cost: 49 value: 75
Generation 4: cost: 50 value: 78
Generation 5: cost: 50 value: 78
Generation 6: cost: 50 value: 78
Generation 7: cost: 50 value: 78
Generation 8: cost: 50 value: 78
Generation 9: cost: 50 value: 78
Generation 10: cost: 50 value: 78
Generation 11: cost: 50 value: 78
Generation 12: cost: 50 value: 78
Generation 13: cost: 50 value: 78
Generation 14: cost: 50 value: 78
Generation 15: cost: 50 value: 78
Generation 16: cost: 50 value: 78
Generation 17: cost: 50 value: 78
Generation 18: cost: 50 value: 78
Generation 19: cost: 50 value: 78
Generation 20: cost: 50 value: 78
Generation 21: cost: 50 value: 78
Generation 22: cost: 50 value: 78
Generation 23: cost: 50 value: 78
Generation 24: cost: 50 value: 78
Generation 25: cost: 50 value: 78
Generation 26: cost: 50 value: 78
Generation 27: cost: 50 value: 78
Generation 28: cost: 50 value: 78
Generation 29: cost: 50 value: 78
Generation 30: cost: 50 value: 78
Generation 31: cost: 50 value: 78
Generation 32: cost: 50 value: 78
Generation 33: cost: 50 value: 78
Generation 34: cost: 50 value: 78
Generation 35: cost: 50 value: 78
Generation 36: cost: 50 value: 78
Generation 37: cost: 50 value: 78
Generation 38: cost: 50 value: 78
Generation 39: cost: 50 value: 78
Generation 40: cost: 50 value: 78
Generation 41: cost: 50 value: 78
Generation 42: cost: 50 value: 78
Generation 43: cost: 50 value: 78
Generation 44: cost: 50 value: 78
Generation 45: cost: 50 value: 78
Generation 46: cost: 50 value: 78
Generation 47: cost: 50 value: 78
Generation 48: cost: 50 value: 78
Generation 49: cost: 50 value: 78
Generation 50: cost: 50 value: 78
Best gene: [1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1]

全探索していないので、どれだけ優秀か分かりませんwが、とりあえずかなりの速度で解を求めることができました。

所感としては、for文が早く動くJuliaはこの手の逐次計算が得意であり、力を発揮しやすいと思います。（今回はかなり簡単なものなので恩恵はあまり感じれていませんw）

ということで、juliaで簡単なアルゴリズムを書いてみましたの記事でした。

最後までお読みいただきありがとうございました。