ナップサック問題の解き方（Wの数値が大きいとき）【分枝限定法】

ナップサック問題よ、私は帰ってきた！

はじめに

• 以前、ナップサック問題の解き方を書いたけど、W(ナップサックの最大許容重量)が大きいとき（$10^9$程度）は対応できなかったから、Wが大きいときも対応できる解法を書く。
• 今回は新しいアプローチを考える。分枝限定法と言うらしいけど、正直、名前なんてどうでも良い。大事なのはアルゴリズムだよね。ちなみに、wiki日本語版へのリンク張ったけど、日本語版の解説は人に理解して貰うつもりゼロの、教科書での定義を単に書き写しただけっぽい感じ。wiki英語版(c++での実装例解説付き)を見た方がマシです。ブラウザの日本語翻訳機能が向上しているから問題なく読めるでしょ。

新しいアプローチ

• 前提では、ナップサックに入れられるN個の荷物は、それぞれ、入れるか入れないか（0個or1個）の判断しか出来ない。
• この条件を緩和して、0個〜1個の間の実数を可能とする前提（これを実数条件と呼び、元の条件を整数条件と呼ぶ）での答えを考える。i番目の荷物の個数[実数]をxi個、i番目の荷物の価値をvi、i番目の荷物の重さをwiとする。
• 当然、この条件で価値を最大化しようとすれば、vi(価値)/wi(重さ){i=0,1,2,...,N-1}の大きい順に並べて、重さの合計がWになるまでナップサックに詰めるだけ。説明の簡略化のため、既にiはvi/wiが降順になるように並べ替え済みとする。
• 仮にΣwi<=Wなら、全ての荷物をナップサックに入れればよいから、答えはΣviとなり、Σwi>Wであれば、Σwi(i in 0..<(n-1)) + wnxn = W となる整数nと実数xnが簡単に算出できるから、その値を用いて、答えはΣvi(i in 0..<(n-1)) + vnxn　となる。
• まず、以前解説したの動的計画法で解けなかった、Wが大きい(>$10^8$)ケースの入力値で試してみる。なお、下記入力値は単位重さあたりの価値の大きい順に並び替え済み。

10 936447862 // 荷物の個数10個、ナップサックの最大許容重量9.3e8
832 42367415 // 1つ目の荷物 価値832，重さ42367415
853 243593312
854 810169801
642 727293784
691 957981784
294 687140254
369 977358410
333 932608409
139 870916042
101 685539955

• 個数を実数(0.3個とかを許容する)とする場合に解を求めるコードは、以下の通り。

let (N,W) = [readLine()!.split(separator:" ").map{Int($0)!}].map{($0[0],$0[1])}[0]
var vwfs:[(v:Int,w:Int,f:Float)] = [] // v:価値、w：重さ、f:価値÷重さ

for _ in 0..<N {
    let (v,w) = [readLine()!.split(separator:" ").map{Int($0)!}].map{($0[0],$0[1])}[0]
    let f = Float(v) / Float(w) // Float値だからf
    vwfs.append((v,w,f))
}

vwfs.sort(){$0.f > $1.f} // 重量あたり価値の大きい順に並び替え。ただし、重量あたり価値が同じ場合は、重さの小さい順。

var xs:[Float] = [] // xs:個数（実数）の配列
var vsum:Float = 0 // 実数条件を満たす価値合計
var vsumN:Int = 0 // 整数条件を満たす価値合計

var res:Int = W // ナップサックの許容残量

for i in 0..<N {
    if res == 0 { // xが0となるケース
        xs.append(0)        
    } else if res >= vwfs[i].w { // xが1となるケース
        res -= vwfs[i].w
        vsum += Float(vwfs[i].v)
        vsumN += vwfs[i].v
        xs.append(1)
    } else { // xが実数となるケース(forループ中、最大1回)
        let x = Float(res) / Float(vwfs[i].w)
        res = 0
        vsum += Float(vwfs[i].v) * x
        xs.append(x)
    }
}

print(xs) // [1.0, 1.0, 0.8029022, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
print(vsum) // 2370.6785
print(vsumN) // 1685

• 上記の通り、実数条件では、i=0,1の荷物およびi=2の荷物の0.8029022個までリュックに詰めたときに解2370.6785となる。
• このとき、整数条件では、i=0,1の荷物をリュックに詰めたときはナップサックの許容残量を満たし、価値1685となる。この価値1685を整数条件での暫定解とする。なお、当然、「実数条件での解」≧「整数条件の暫定解」となる。
• そして、最終的に求めるべきナップサック問題の解に対しては、「実数条件での解」≧「ナップサック問題の解」≧「整数条件の暫定解」の関係がある。
• 次の段階では、i=2の荷物のみ整数条件とする2パターン(x2=0or1)について、それぞれのパターンで、i=0,1,3,...,9の荷物の個数のみ実数条件として、xsの組み合わせを考える。
• 一般化のため、求める個数配列xiについて、特定のiについて0or1を固定したときの「xs(整数条件)」「実数条件の解」「整数条件の解」「xs(実数条件)のときにxiの値が実数となるインデックスi」を求める関数solve()を作成した。
• solveへ与える引数は整数配列で、配列の値は、「1:個数を1で固定」「0:個数を0で固定」「-1:個数が未定」とする。例えば、最初の状態では、全て未定のため、[-1,-1,-1,-1,-1,-1,-1,-1,-1,-1]となり、x2=0のときは、[-1,-1,0,-1,-1,-1,-1,-1,-1,-1]、x2=1のときは、[-1,-1,1,-1,-1,-1,-1,-1,-1,-1]
• ついでに、元のプログラムでvar vwfs:[(v:Int,w:Int,f:Float)] = [] // v:価値、w：重さ、f:価値÷重さとしてた配列vwfsのfは削った。並び替えのsortを調整すれば済む話なので、余計なモノを削ったということ。

let (N,W) = [readLine()!.split(separator:" ").map{Int($0)!}].map{($0[0],$0[1])}[0]
var vws:[(v:Int,w:Int)] = [] //v:価値、w：重さ

for _ in 0..<N {
    let (v,w) = [readLine()!.split(separator:" ").map{Int($0)!}].map{($0[0],$0[1])}[0]
    vws.append((v,w))
}

vws.sort{Float($0.0)/Float($0.1) > Float($1.0)/Float($1.1)}

typealias zeroOneNonfixs = [Int] // リュックに入れないことが確定している場合は0，入れることが確定している場合は1，未定の場合は-1とする。
var zons_root = zeroOneNonfixs(repeating:-1,count:N) 

func solve(_ zons:zeroOneNonfixs) -> ([Int],Float,Int,Int,Float)? { // zonは、zeroOneNonfixed型。なお、W,vwfsはグローバル定数として取得するため、引数に含めない

    var xs = [Int](repeating:0,count:N) // xs:個数（整数条件）の配列-解答①
    var vsum:Float = 0 // 実数条件を満たす価値合計-解答②
    var vsumN:Int = 0 // 整数条件を満たす価値合計 -解答③
    var k:Int = -1 // 解答②の実数条件で実数が格納されるインデックス -解答④  
    var f:Float = 0 // xs:個数（実数条件）の配列のインデックスkの要素の値-解答⑤
    var ans:([Int],Float,Int,Int,Float)? // 解答(①②③④)を格納する変数
    
    var ivws_r:[(i:Int,v:Int,w:Int)] = [] // zonの未定部分だけで作成する配列
    var res:Int = W // ナップサックの許容残量
    
    for (i,zon) in zons.enumerated() {
        let (v,w) = vws[i]
        switch zon {
            case -1:
                ivws_r.append((i,v,w))   
            case 1:
                res -= w
                vsum += Float(v)         
                vsumN += v
                xs[i] = 1
            default:break
        }
    } 
    
    if res < 0 {return ans} //　解なしで解答
    
    for (i,v,w) in ivws_r {
        if res >= w {
            res -= w
            vsum += Float(v)
            vsumN += v
            xs[i] = 1
        } else {
            let x = Float(res) / Float(w)
            res = 0
            vsum += Float(v) * x
            k = i
            f = x
            break
        }
    }

    ans = (xs,vsum,vsumN,k,f)

    return ans
}

print(zons_root) // [-1, -1, -1, -1, -1, -1, -1, -1, -1, -1]
print(solve(zons_root) ?? "解なし") // ([1, 1, 0, 0, 0, 0, 0, 0, 0, 0], 2370.6785, 1685, 2, 0.8029022)
                                    // 「解なし」となるのは、zonsの1となっている荷物の重さだけでWを超過するとき

• よし！元のプログラムと結果が一致することを確認！
• じゃあ、次は、i=2の荷物のみ整数条件とする2パターン(x2=0or1)をとく。

var zons_k2_0 = zons_root
var zons_k2_1 = zons_root
zons_k2_0[2] = 0
zons_k2_1[2] = 1

print(solve(zons_k2_0) ?? "解なし") // ([1, 1, 0, 0, 0, 0, 0, 0, 0, 0], 2259.201, 1685, 3, 0.8943939)
print(solve(zons_k2_1) ?? "解なし") // ([1, 0, 1, 0, 0, 0, 0, 0, 0, 0], 1979.8331, 1686, 1, 0.34447023)

• 上記のように最初の条件rootから、k2_0,k2_1のパターンを作ったことを分枝すると呼ぶ。枝を増やすとも言う。
• このときの出力をみると、整数条件を満たす解(k2_0:1685,k2_1:1686)のうち、暫定解(root:1685)を超える数値(k2_1:1686)があるので、最も大きい数値(1868)を暫定解としてアップデートする。実数条件をみたす解(k2_0:2259.201,k2_1:1979.8331)はまだ暫定解(1868)を超えているので、k2_0,k2_1のパターンはまだ残す。
• もし、分枝したパターンの実数条件解が暫定解以下となった場合は、そのパターンの先は分枝しても整数解が暫定猪飼を超えることは無いので、分子不要となる。そのように「実数条件<暫定解」となったパターンを捨てることを「枝を刈る」と表現する。この枝を刈ることを以て、「分枝限定法」と呼ばれる。
• それでは、k2_0,k2_1はまだ枝を刈らないので、更に分枝する。分枝は、それぞれのkの位置(k2_0:インデックス3 , k2_1:インデックス1)で行う。

// k2_0から分枝
var zons_k2_0__k3_0 = zons_k2_0
var zons_k2_0__k3_1 = zons_k2_0
zons_k2_0__k3_0[3] = 0
zons_k2_0__k3_1[3] = 1

print(solve(zons_k2_0__k3_0) ?? "解なし") // ([1, 1, 0, 0, 0, 0, 0, 0, 0, 0], 2154.2017, 1685, 4, 0.67901826) -- まだ、枝刈りしない
print(solve(zons_k2_0__k3_1) ?? "解なし") // ([1, 0, 0, 1, 0, 0, 0, 0, 0, 0], 2058.0432, 1474, 1, 0.6846931)  -- まだ、枝刈りしない


// k2_1から分枝
var zons_k2_1__k1_0 = zons_k2_1
var zons_k2_1__k1_1 = zons_k2_1
zons_k2_1__k1_0[1] = 0
zons_k2_1__k1_1[1] = 1

print(solve(zons_k2_1__k1_0) ?? "解なし") // ([1, 0, 1, 0, 0, 0, 0, 0, 0, 0], 1760.07, 1686, 3, 0.115373805)  -- まだ、枝刈りしない
print(solve(zons_k2_1__k1_1) ?? "解なし") // 解なし

• これ以上を手作業で行うのは面倒なので、キューを用いて行う。とはいえ、配列を使った「なんちゃてキュー」を使用する。

// print(zons_root) // [-1, -1, -1, -1, -1, -1, -1, -1, -1, -1]

var queue:[zeroOneNonfixs] = []

var max_vsumN = 0 // 解答用の価値総額
var ans_xs = zons_root // 解答用のzons配列

queue.append(zons_root)
while !queue.isEmpty {
    let zons = queue.removeFirst()
    guard let (xs,vsum,vsumN,k,f) = solve(zons) else {continue} // 解が無ければスキップ -- 枝刈り①
    // print(xs,vsum,vsumN,k)
    if vsumN > max_vsumN { max_vsumN = vsumN ; ans_xs = xs } //整数解が暫定解を超えたら、暫定解を塗り替え
    if vsum < Float(max_vsumN) {continue} // 実数解が暫定解未満の場合、その枝に将来はないのでスキップ -- 枝刈り②
    
    //分枝作成
    if k < 0 {continue} // 分枝せずにスキップ
    var xs_k_0 = zons
    var xs_k_1 = zons
    xs_k_0[k] = 0
    xs_k_1[k] = 1
    queue.append(xs_k_0)
    queue.append(xs_k_1)
}

print(max_vsumN)
// print(ans_xs) // [1, 0, 1, 0, 0, 0, 0, 0, 0, 0]

• 完成だ！

実際の問題に適用してみる

• 以前、ＡＣ出来なかったABC032のd問題で、前記のコードを提出したら、40ms位で通った！よかったよかった。
• でも、写経がわりに、「Pythonによる 数理最適化入門」を読んで書いた、下記コードの方が速くて、1/3以下の時間で終わった。

struct KnapProb { // 暫定解
    let capacity:Int //許容重量
    let items:[(v:Int,w:Int)] //荷物の価値(v)と重さ(w)：ソート済みで格納
    let ratio:[Double] // 各荷物のv/wの値
    var zeroOnes:[Int] // 確定した条件を格納(確定前は-1)
    var lb:Int = -1 // 下界（貪欲法での解）
    var ub:Double = -1.0 // 上界（線形緩和での解）
    var xlb:[Int] // 下界を得た条件
    var xub:[Double] // 上界を得た条件(実数となる可能性があるためFloat)
    var bi:Int = -1 // 線形緩和の解を得たときに実数となった荷物の番号

    init(_ capacity:Int,_ items:[(v:Int,w:Int)],_ zeroOnes:[Int] = []){
        self.capacity = capacity
        self.items = items
        var ratio:[Double] = []
        for (v,w) in items {
            ratio.append(Double(v) / Double(w))
        }
        self.ratio = ratio
        if zeroOnes.isEmpty {
            self.zeroOnes = [Int](repeating:-1,count:items.count)
        } else {
            self.zeroOnes = zeroOnes
        }
        self.xlb = [Int](repeating:0,count:items.count)
        self.xub = [Double](repeating:0,count:items.count)
    }
    
    mutating func get_bounds(){ // 　zeroOnesを前提として、lbとubを計算する。その過程で、xlb,xub,biも計算する
        var wsum = 0
        var ritems:[(index:Int,v:Int,w:Int)] = []
        for (i,e) in zeroOnes.enumerated() { // zeroOnesの条件をxlb,xluに反映
            switch e {
                case 1:
                    xlb[i] = 1
                    xub[i] = 1
                    wsum += items[i].w
                case 0:
                    xlb[i] = 0
                    xub[i] = 0
                case -1:
                    ritems.append((i,items[i].v,items[i].w))
                default:break
            }
        }
        var rcap = capacity - wsum // 残り許容重量        
        if rcap < 0 {lb = -1;ub = -1;return} // zeroOnes条件の解なし
        for (i,v,w) in ritems {
            if w <= rcap {
                rcap -= w
                xlb[i] = 1
                xub[i] = 1
            } else {
                xub[i] = Double(rcap) / Double(w)
                bi = i
                break
            }
        }
        lb = 0
        for (i,x) in xlb.enumerated() { lb += Int(x) * items[i].v }
        ub = 0
        for (i,x) in xub.enumerated() { ub += x * Double(items[i].v) }
    }
}

func kpSolve(_ capacity:Int,_ items:[(v:Int,w:Int)])->(Int,[Int]){
    var queue:[KnapProb] = []
    var root = KnapProb(capacity,items)
    root.get_bounds() // bestは既知の暫定解のうち、最大の下限を持つ暫定解
    var best = root
    queue.append(best)
    while !queue.isEmpty {
        // print(queue.count)
        var p = queue.removeFirst()
        p.get_bounds()
        // print("best:",best.bi,best.zeroOnes)
        // print("試行p:",p.bi,p.zeroOnes)
        if p.ub <= Double(best.lb) {continue} //試行pの上限がbestの下限以下の時は無視（枝切り）
        // 以下は、p.ub > best.lb の時
        
        // best更新
        if p.lb > best.lb {best = p} // 試行ｐの下限がbestの下限を超えたので、best更新
        
        // 子問題作成(=分枝作業)
        if p.ub > Double(p.lb) {
            let k = p.bi
            var temp01_k0 = p.zeroOnes
            temp01_k0[k] = 0
            var temp01_k1 = p.zeroOnes
            temp01_k1[k] = 1            
            let p0 = KnapProb(capacity,items,temp01_k0)
            let p1 = KnapProb(capacity,items,temp01_k1)
            queue.append(p0)
            queue.append(p1)

        }
    }
    
    return (best.lb,best.xlb)
    
}
 
let (N,W) = [readLine()!.split(separator:" ").map{Int($0)!}].map{($0[0],$0[1])}[0] 

var items:[(Int,Int)] = []
for _ in 0..<N {
    let (a,b) = [readLine()!.split(separator:" ").map{Int($0)!}].map{($0[0],$0[1])}[0]
    items.append((a,b))
}

items.sort{Float($0.0)/Float($0.1) > Float($1.0)/Float($1.1)}
// print(items)
// for (v,w) in items {
//     print(Float(v)/Float(w))
// }

let ans = kpSolve(W,items)
print(ans.0)
print(ans.1)//

• ロジック自体は何も変わりないから、インスタンス内にsolve関数をおいて、インスタンス内のプロパティに対して、計算したから？もしくは、Floatを使わずにDoubleを使ったから？
• 検証が面倒なので、何もしてないけど、不思議だなぁ。

さいごに

• 以前､ギブアップしてから2か月くらいかかっちゃったよ。
• まぁ、理解しちゃえば、大したことしてないけどね。端的に言えば、一番最初のbit全探索について、分岐を枝刈りしただけだからね。
• まあ、終わって良かった。
• スッキリしたから、AtCoderのd問題以下の過去文を解く地道な作業に戻るよ。