LoginSignup
16
14

More than 5 years have passed since last update.

@_likr(Yosuke Onoue)

GPUでナップザック問題を解く話 GPGPU AdC2012編

Last updated at Posted at 2012-12-13

執筆者の忙しさが＜加速＞されてる気がするGPGPU Advent Calendarの13日目です。

「GPGPUでこんなの速くしてみた」っていうネタがまだ少ないので書くことにします。関西GPGPU勉強会 #2にて、こんな感じのタイトルで発表を行いましたが、今回はそれとは異なったアプローチの内容です。

ナップザック問題とは

Wikipedia先生のナップサック問題のページによるとこう書いています。

容量 C のナップサックが一つと、n 個の品物（各々、価値 p_i, 容積 c_i）が与えられたとき、ナップサックの容量 C を超えない範囲でいくつかの品物をナップサックに詰め、ナップサックに入れた品物の価値の和を最大化するにはどの品物を選べばよいか

競技プログラミング系ではより単純に、最大化された価値の和(最適値や最適目的関数値と呼ばれます)を求めよという問題が出題されることがあり、SPOJ 3321なんかはその例になっています。

ナップザック問題を解く代表的な方法は動的計画法といって、Web上でもたくさん解説があるほか、これについて載ってるカジュアルな書籍としてはプログラミングコンテストチャレンジブックあたりが有名でしょうか。

とりあえずGPGPUで解いてみる

とりあえず後者の最適目的関数値だけを求めるパターンのCPU版を実装してみるとザッとこんな感じになります。Item構造体は以降のソースで共通です。

CPU版
struct Item
{
    int p;
    int c;
};


int cpu_dp(const int n, const Item* items, const int C)
{
    std::vector<int> buffer(C + 1);
    for (int j = 0; j < n; ++j) {
        const Item item = items[j];
        for (int i = C; i >= item.c; --i) {
            const int v0 = buffer[i];
            const int v1 = buffer[i - item.c] + item.p;
            if (v1 > v0) {
                buffer[i] = v1;
            }
        }
    }
    return buffer.back();
}

さらにサクッとGPU版も実装します。環境はCUDA+Thrustです。

GPU版
const int THREAD_NUM = 512;


__global__ void dp_step(
        const int* buffer0,
        int* buffer1,
        const Item item,
        const int C)
{
    const int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (item.c <= i && i <= C) {
        const int v0 = buffer0[i];
        const int v1 = buffer0[i - item.c] + item.p;
        buffer1[i] = max(v0, v1);
    }
}


int cuda_dp(const int n, const Item* items, const int C)
{
    const int size = C + 1;
    thrust::device_vector<int> buffer0(size, 0);
    thrust::device_vector<int> buffer1(size, 0);

    const int thread_num = THREAD_NUM;
    const int block_num = size / thread_num + (size % thread_num != 0);
    for (int j = 0; j < n; ++j) {
        dp_step<<<block_num, thread_num>>>(
                thrust::raw_pointer_cast(buffer0.data()),
                thrust::raw_pointer_cast(buffer1.data()),
                items[j], C);
        buffer0.swap(buffer1);
    }
    return buffer0.back();
}

計算時間の比較をします。使用したマシンのCPUがXeon W3520、GPUがTesla C2050です。また、計算するナップザック問題は品物数nが10,000個、ナップザックの容量Cが約2,500,000のほんのちょっと大きな問題です。

解法 計算時間（秒）
CPU 33.751
GPU 2.179

計算時間の差が15.5倍とよくありそうな結果になりました。

解のトレース

ここまでで最適目的関数値はそこそこ高速に計算できることがわかりました。では、最初の問題定義に戻って「どの品物を選べばよいか」まで求めてみましょう。プログラムの大筋はそれぞれ上のプログラムと変わりません。

CPU版（解トレースあり）
int cpu_dp_sol(const int n, const Item* items, const int C)
{
    std::vector<int> buffer(C + 1);
    const int solution_size = n / 32 + (n % 32 != 0);
    std::vector<std::vector<int> > solution(solution_size, std::vector<int>(C + 1));
    int sol_index = 0;
    for (int j = 0; j < n; ++j) {
        const Item item = items[j];
        for (int i = C; i >= item.c; --i) {
            solution[sol_index][i] <<= 1;
            const int v0 = buffer[i];
            const int v1 = buffer[i - item.c] + item.p;
            if (v1 > v0) {
                buffer[i] = v1;
                solution[sol_index][i] &= 1;
            }
        }
        if (j % 32 == 31) {
            ++sol_index;
        }
    }
    return buffer.back();
}
GPU版（解トレースあり）
__global__ void dp_sol_step(
        const int* buffer0,
        int* buffer1,
        const Item item,
        const int C,
        int* solution)
{
    const int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (item.c <= i && i <= C) {
        const int v0 = buffer0[i];
        const int v1 = buffer0[i - item.c] + item.p;
        buffer1[i] = max(v0, v1);
        solution[i] = (solution[i] << 1) & (v1 > v0);
    }
}


int cuda_dp_sol(const int n, const Item* items, const int C)
{
    const int size = C + 1;
    thrust::device_vector<int> buffer0(size, 0);
    thrust::device_vector<int> buffer1(size, 0);
    thrust::device_vector<int> solution(size, 0);
    const int solution_size = n / 32 + (n % 32 != 0);
    std::vector<std::vector<int> > host_solution(solution_size, std::vector<int>(size));

    const int thread_num = THREAD_NUM;
    const int block_num = size / thread_num + (size % thread_num != 0);
    int sol_index = 0;
    for (int j = 0; j < n; ++j) {
        dp_sol_step<<<block_num, thread_num>>>(
                thrust::raw_pointer_cast(buffer0.data()),
                thrust::raw_pointer_cast(buffer1.data()),
                items[j], C,
                thrust::raw_pointer_cast(solution.data()));
        if (j % 32 == 31) {
            thrust::copy_n(solution.begin(), size, host_solution[sol_index].begin());
            ++sol_index;
        }
        buffer0.swap(buffer1);
    }
    thrust::copy_n(solution.begin(), size, host_solution[sol_index].begin());
    return buffer0.back();
}

solutionテーブルから解を復元する処理は実装していませんが、CPU版GPU版共通になるので気にしないことにしましょう。計算時間はそれぞれ以下のようになりました。

解法 計算時間（秒）
CPU 82.563
GPU 31.648

残念ながらこのケースでは、最適目的関数値を求める場合と比較して大量のホスト側メモリとデータ転送を必要とし、GPUの場合でも計算時間が大幅にかかってしまいました。それにしても時間がかかりすぎている気はしますが、ホスト側の環境も影響してるようで、これもこの記事の本筋ではないので放置することにしましょう。

解のトレースを高速化する

プログラムのチューニングの基本は「遅い部分を速くする」ことです。解のトレースには解情報の保持や転送に大幅な時間がかかっていて、実装面での高速化はあまり期待できなさそうです。そこで、アルゴリズムに工夫をし、根本的に保持しなければならない解情報を減らす方法を考えてみましょう。

ナップザック問題では、品物をあらかじめ、価値を容積で割った値（効用値）でソートしておくと、効用値の大きい側の大部分の品物はナップザックに詰める品物として選ばれ、逆に小さい側の大部分はナップザックに詰められないということが経験的にわかっています。この観察に基づけば、ナップザック問題で考慮しなければいけない品物の個数は全品物の個数と比べると非常に少ないことになり、この考慮しなければいけない品物の部分をコアと呼びます。

コアの考え方を使うと、最初に最適目的関数値を求めれば、コアの周りの少数の品物でだけ解のトレースを行えばよくなり、計算時間も消費するメモリ量も大幅に抑えられるでしょう。コアのサイズは品物数が10000ぐらいだとsqrt(n) ~ 2 * sqrt(n)ぐらいあれば十分です。

GPU版（解トレース改良）
struct ItemComparator
{
    __host__ __device__ bool operator()(const Item& item1, const Item& item2)
    {
        return item1.p * item2.c > item2.p * item1.c;
    }
};


int cuda_dp_core(const int n, const Item* items, const int C)
{
    const int opt = cuda_dp(n, items, C);

    thrust::device_vector<Item> dev_items(items, items + n);
    thrust::sort(dev_items.begin(), dev_items.end(), ItemComparator());
    thrust::host_vector<Item> host_items(dev_items);

    int s;
    int sum_c = 0;
    for (int j = 0; j < n; ++j) {
        sum_c += host_items[j].c;
        if (sum_c > C) {
            s = j;
            break;
        }
    }

    int delta = (int)std::sqrt(n);
    const int core_inc = 10;
    int core_opt;
    do  {
        const int core_size = delta * 2;
        const Item* core_start =
            thrust::raw_pointer_cast(host_items.data()) + s - delta;

        int p_offset = 0;
        int core_C = C;
        const int stop = s - delta;
        for (int j = 0; j < stop; ++j) {
            p_offset += host_items[j].p;
            core_C -= host_items[j].c;
        }

        core_opt = cuda_dp(core_size, core_start, core_C) + p_offset;

        delta += core_inc;
    } while (core_opt != opt);
    return opt;
}

実行時間をこれまでの2つのGPU版プログラムと比較してみましょう。

解法 計算時間（秒）
GPU（解トレースなし） 2.179
GPU（解トレースあり） 31.648
GPU（解トレース改良） 2.184

解トレースにかかる時間が大幅に短縮されて、解トレースなしの場合とほとんど変わらない時間になりました。もうちょっとひねったら論文とかにもできるんじゃないですかね。（適当

最後に

GPGPU、ひいては広義のHPCには、プロセッサの知識だけでなく問題領域の知識も重要だと感じます。そういう意味でも、プロセッサや並列計算の専門家が世のあらゆる問題領域のプログラムを高速化するといったことは不可能で、逆に問題領域の専門家がプロセッサの知識を持って自分の問題領域のプログラムを高速化しなければいけないと私は思っています。近年のGPGPUの文化は問題領域の専門家にHPCの門口を広げてくれたという点で非常に大きな役割をしてくれているように思えます。Xeon Phiとかもありますが、今後のGPGPUの動向もまだまだ楽しみですね。

16
14
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
16
14