連鎖行列積問題の考え方

連鎖行列積問題の解法の考え方についてまとめます

例題として、AIZU ONLINE JUDGEの問題を引用します。

問題

n個の行列の連鎖 M1,M2,M3,...,Mn が与えられたとき、
スカラー乗算の回数が最小になるように積 M1M2M3...Mnの計算順序を決定する問題
を連鎖行列積問題(Matrix-Chain Multiplication problem)と言います。

n個の行列について、行列 Mi の次元が与えられたとき、積 M1M2...Mn
の計算に必要なスカラー乗算の最小の回数を求めるプログラムを作成してください。

入力

入力の最初の行に、行列の数 nが与えられます。
続く n 行で行列 Mi(i=1...n) の次元が空白区切りの２つの整数 r、c で与えられます。
r は行列の行の数、cは行列の列の数を表します。

出力

最小の回数を１行に出力してください。

行列の積の考え方

下図のようなl × mの行列Aとm × nの行列Bの積はl × nとなり、
各要素の$C_{ij}$は次の式で得られる

$C_{ij} = \displaystyle \sum_{ i = 1 }^{ m } a_{ik}b_{kj}$

なので、この計算では、l × m × n回のかけ算が必要になる。

複数行列積について考える

下図のような$(M_1,M_2,M_3..M_6)$のような、$M_i$が$p_{i-1}行$$p_i$行の行列であるような、
n個の行列のかけ算を考えます。

これらの行列は様々な順番で計算できるので、行列同士のかけ算の順番によって、
かけ算の計算回数が変わってきます。

連鎖行列問題の解き方考察

連鎖行列において、全ての可能な順番をを調べる方法は$0(n!)$になります。
ただ、小さな部分問題として分割できるので、動的計画法を使うことができます。

動的計画法で解くための考え方

$(M_1,M_2)$を計算するための方法は1通りであり、$p_0 × p_1 × p_2$回のかけ算が必要です。
同様に$(M_2,M_3)$を計算する方法も1通りで、$p_1 × p_2 × p_3$回のかけ算が必要です。
これらのかけ算をコストとして、記録しておきます。

次に、$(M_1,M_2,M_3)$、$(M_2,M_3,M_4)$、...$(M_{n - 2},M_{n - 1},M_n)$ を計算するための最小計算回数つまり最適解を求めます。

例として、$(M_1,M_2,M_3)$の最適解を求めるには、$(M_1(M_2×M_3))$と$((M_1×M_2)M_3))$のコストを計算し、
そのうちの小さい方のコストを$(M_1,M_2,M_3)$のコストとして記録します。

計算量を考えると、

$(M_1(M_2×M_3))$のコスト = $(M_1)$のコスト + $(M_2M_3)$のコスト + $p_0 × p_1 × p_3$
$((M_1×M_2)M_3))$のコスト = $(M_1M_2)$のコスト + $(M_3)$のコスト + $p_0 × p_2 × p_3$

となります。

上記の計算に関して、$(M_1M_2)$のコストと$(M_2M_3)$のコストは記録してあるので、再計算する必要がありません。

また、$1 \leqq i \leqq n$について、$(M_i)$のコストは0になります。

連鎖行列の最適解

一般に連鎖行列の最適解$(M_iM_{i + 1}...M_j)$の最適解は、
$i \leqq k < j$における$(M_iM_{i + 1}...M_k)(M_{k + 1}...M_j)$の最小コストになります

例をあげると
$(M_1M_2M_3M_4M_5)(i = 1,j = 5)$の最適解は

$(M_1)(M_2M_3M_4M_5)$のコスト = $(M_1)$のコスト + $(M_2M_3M_4M_5)$のコスト + $p_0 × p_1 × p_5$(k = 1)のとき

$(M_1M_2)(M_3M_4M_5)$のコスト = $(M_1M_2)$のコスト + $(M_3M_4M_5)$のコスト + $p_0 × p_2 × p_5$(k = 2)のとき

$(M_1M_2M_3)(M_4M_5)$のコスト = $(M_1M_2M_3)$のコスト + $(M_4M_5)$のコスト + $p_0 × p_3 × p_5$(k = 3)のとき

$(M_1M_2M_3M_4)(M_5)$のコスト = $(M_1M_2M_3M_4)$のコスト + $(M_5)$のコスト + $p_0 × p_4 × p_5$(k = 4)のとき

のうちの最小値になります。

上記のアルゴリズムの実装例

以下のような役割の変数を定義します。

m[n + 1][n + 1]	m[i][j]を$(M_iM_{i + 1}...M_j)$を計算するための 最小のかけ算の回数とする2次元配列
p[n + 1]	$M_i$がp[i - 1] × p[i]の行列となるような1次元配列

これらの変数を用いて、m[i][j]の値を以下のように設定します。

if(i == j){
  m[i][j] = 0 
}

if(i < j){
 m[i][j] = $min_{i \leqq k < j}\{m[i][k] + m[k + 1][j] + p[i - 1] × p[k] × p[j]\}$
}

これを疑似アルゴリズムで書くと以下のようになります。

matrixMultiplication()
  for i = 1 to n
    m[i][j] = 0

  for l = 2 to n
    for i = 1 to n - l + 1
      j = i + l - 1
      m[i][j] = INFTY
      for k = i to j - 1
      m[i][j] = min(m[i][j],m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j])

計算量の考察

上記の疑似コードでは、対象とする行列の数lを2からnまで増やしながら、その範囲を指定するiとjを決めています。
さらにその中で、iからjの範囲でkの位置を決めています。
forの3重ループになっているので、$O(n^3)$になります。

例題の解答

上記をふまえて、例題の解答をc++で書きます。

#include <stdio.h>
#include <iostream>
#include <algorithm>

using namespace std;
#define N 100

int main(void){
    int n;
    int p[N + 1],m[N + 1][N + 1] = {};
    
    cin >> n;
    for(int i = 1;i <= n;++i){
        cin >> p[i - 1] >> p[i];
    }
    for(int l = 2;l <= n;++l){
        for(int i = 1;i <= n - l + 1;++i){
            int j = i + l - 1;
            m[i][j] = __INT_MAX__;
            for(int k = i;k <= j - 1;++k){
                m[i][j] = min(m[i][j],m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j]);
            }
        }
    }
    cout << m[1][n] << endl;
}

参考

プログラミングコンテスト攻略のためのアルゴリズムとデータ構造

最小コストソートのアルゴリズムと例題とc++の解答

最小コストソートの考え方を理解できるようにまとめます。

最小コストソートを利用する問題

AIZU OINLINE JUDGEの最小コストソートに関する問題を例題として利用します。

問題

重さ wi(i=0,1,...,n−1) の n 個の荷物が１列に並んでいます。
これらの荷物をロボットアームを用いて並べ替えます。
１度の操作でロボットアームは荷物 i と荷物 j を持ち上げ、それらの位置を交換することができますが、
wi+wj のコストがかかります。
ロボットアームは何度でも操作することができます。

与えられた荷物の列を重さの昇順に整列するコストの総和の最小値を求めてください。

入力

１行目に整数 n が与えられます。２行目に n 個の整数 wi(i=0,1,2,...n−1) が空白区切りで与えられます。

出力

最小値を１行に出力してください。

入力例 1

5
1 5 3 4 2

出力例 1

7

最小コストソートの考え方

例として、{4,3,2,7,1,6,5}という数列の最小コストを考えます。

この数列において、各々の数値が最終的にどの位置に置かれればいいかということに絞って、 次のような入れ替えを行う。

まず先頭の4からを3番目のindexの7と交換して、

7,3,2,4,1,6,5

続いて、7を最大index6の5と交換して

5,3,2,4,1,6,7

続いて、5をindex4の1と交換して、

1,3,2,4,5,6,7

続いて、対象となる1はすでにソートの位置に置かれているので、ここで一連の流れが一旦終了することになります。

先頭から順にソートする作業を再開する、次のindex1の3に着目してindex2の2と交換する

1,2,3,4,5,6,7

2は適切な位置に置かれているので、ここで一連の流れが終了します。

つづいて、まだソート完了フラグが立っていない、index5の6に着目しますが、
6は正しい位置に置かれているので、これでソート終了ということになります。

ここで、行った処理の流れをまとめると、

①4 → 7 → 5 → 1 → 4
②3 → 2 → 3
③6 → 6

この①②③というサイクルにおいての最小コストについて考えます。

サイクルの長さが1のケース

③のようなケースでは入れ替える必要がないので、コストは0になります。

サイクルの長さが2のケース

②のようなサイクルの長さが2のケースでは、1回の交換でそれぞれの要素が最終位置に到達するので、
対象となる2つの要素の和がコストになります。
②のケースだと3 + 2 = 5がコストになる。

サイクルの長さが3以上のケース

①のケース{4 7 5 1}を考えます。
①のケースでは、以下のようにこのサイクルにおける最小コスト1を使って、それ以外の要素を移動することで、最小コストの19になります。

cost = 1 + 5 = 6
4 7 1 5

cost = 1 + 7 = 8
4 7 5 1 

cost = 1 + 4 = 5
4 1 5 7 

総コストはcost = 6 + 8 + 5 = 19となる

つまり、このケースでは、サイクルの中で最小値を使って、それ以外の要素を移動することが、
最適な方法となる。

最小コスト移動の総コストを考える

サイクルの中の要素を$w_i$サイクル内の要素の個数をnとすると、
コストの総和 + (最小値 * (個数 - 2))と一致するので、

$\displaystyle \sum w_i + (n - 2) * min(w_i)$

となります。

最小コスト移動を行なった場合の各サイクルの総移動コスト

サイクル①のコストは0サイクル②のコストは5、
サイクル③のコストは19なので、この数列の最小総移動コストは24ということになる。

が、この考え方には反例があります。

反例を考える

数列{2,1,8,10,7,9}を例に先ほどのアルゴリズムを適応してみます。

まず、index0の先頭2とindex1の1を交換します。

1 2 8 10 7 9

交換した1は交換する必要がないので、コスト3となりここでこのサイクルは終了

続いて、8 10 7 9についてソートします。

ここではサイクルの長さが4になるので、最小コストソートを利用します。

まず、7と9を交換(コスト16)

8 10 9 7

続いて、7と10を交換(コスト17)

8 7 9 10

最後に、7と8を交換(コスト15)

7 8 9 10

総コストは、3 + (16 + 17 + 15) = 51となります。

最小コストを他のサイクルから借りてくる方法

8 → 10 → 9 → 7 → 8のサイクルについて、7と1を交換して、
8 → 10 → 9 → 1 → 8のサイクルとして、コストの総和を考えると

28 + 2 * 1に1と7の2回分の交換を加えて、このサイクルのコストは46となり、
合計しても49となるので、先ほどのアルゴリズムよりもコストが低くなります。

つまり、サイクルの外から借りてきた要素を使って、移動した方がコストが小さくなる場合があります。

最小コストをサイクルの外から借りてきた場合の計算量を考える

サイクルの外の要素を最小値をxとすると、
$2 * (min(w_i) + x$だけコストが増加するが、
$(n - 1) * (min(w_i) - x)$だけコストが減ります。
この時のコストは

$\displaystyle \sum w_i + (n - 2) * min(w_i) + 2 * (min(w_i) + x - (n - 1) * (min(w_i) - x)$ $ = \displaystyle \sum w_i + min(w_i) + (n + 1) * x$

となる。

以上を考えると、各サイクルのコストは、要素全体の最小値を借りた場合と借りない場合を計算して、
コストが小さい方を採用する必要があります。

最小コスト問題の実装のポイント

上記の仕組みを利用して、プログラムを組む上でのポイントをまとめます。

ソートされた数値がどの位置に置かれるか保存する配列を持つ

ソートされた数列の各々の数値がどの位置に入るかを記録することで、
数値を交換する動作を計算なしで行うことができます。

例えば、以下のようなコードを考えます。

#define MAX 1000
// 最大入力値
#define VMAX 10001

int n;
//数列入力用
int A[MAX];
// 入力数列の中の最小値
int s;
// 入力したAをソートした配列
int B[MAX];
// 数値をINDEXとしてソートした数列がどの位置に置かれるのかを記録する
int T[VMAX];

以下のように、ソートされた数列の各々の数値が何番目に置かれるのかを記録します。
計数ソートの考え方と同じです。

for(int i = 0;i < n;++i){
    B[i] = A[i];
}
// Bをソートする
sort(B,B + n);

for(int i = 0;i < n;++i){
    // ソートされた並び順ごとに順番を入れておく
    // ソートされた数列が1 3 5 7 9とするとT[1] = 0,T[3] = 1のように順番が入る
    T[B[i]] = i;
}

配列T[]は交換するサイクルにおいて、以下のように使えます。

for(int i = 0;i < n;++i){
    if(V[i]) continue;
    // 初回は0が対象になる
    int cursor = i;
    // 1サイクルでかかるコスト
    int S = 0;
    // 1サイクルにおける最小コスト、取りうる最大数値にしておく
    int m = VMAX;
    // サイクルの長さ
    int an = 0;
    // 1サイクル処理
    while(1){
        // ソート済みフラグ
        V[cursor] = true;
        an++;
        // 対象のコスト
        int v = A[cursor];
        // 最小コスト
        m = min(m,v);
        // 最小の値においての交換を加味しないサイクルの合計コスト
        S += v;
        // 対象が置かれるべき位置を割り出すT[v]に置かれるべき位置を入れる
        cursor = T[v];
        // サイクルが終了していたら、次に移る
        if(V[cursor]) break;
    }
    // サイクル内の最小値を利用したコスト計算(Σwi + (n - 2) * min(wi)
    int cost1 = S + (an - 2) * m;
    // 数列の最小コストを借りた場合のコスト(Σwi + min(wi) + (n + 1) * min(全数列)
    int cost2 = S + m + (an + 1) * s;
    // 最小値のコストを外から借りた場合とコストを比較する
    ans += min(cost1,cost2);
}

全体のソースコード

c++の全体のコードを載せます。

#include <iostream>
#include <algorithm>
using namespace std;

#define MAX 1000
// 最大入力値
#define VMAX 10001

int n;
//数列入力用
int A[MAX];
// 入力数列の中の最小値
int s;
// 入力したAをソートした配列
int B[MAX];
// 数値をINDEXとしてソートした数列がどの位置に置かれるのかを記録する
int T[VMAX];

int solve(){
    // コスト
    int ans = 0;
    // 交換終了フラグ
    bool V[MAX] = {};
    for(int i = 0;i < n;++i){
        B[i] = A[i];
    }
    // Bをソートする
    sort(B,B + n);
    
    for(int i = 0;i < n;++i){
        // ソートされた並び順ごとに順番を入れておく
        // ソートされた数列が1 3 5 7 9とするとT[1] = 0,T[3] = 1のように順番が入る
        T[B[i]] = i;
    }
    
    for(int i = 0;i < n;++i){
        if(V[i]) continue;
        // 初回は0が対象になる
        int cursor = i;
        // 1サイクルでかかるコスト
        int S = 0;
        // 1サイクルにおける最小コスト、取りうる最大数値にしておく
        int m = VMAX;
        // サイクルの長さ
        int an = 0;
        // 1サイクル処理
        while(1){
            // ソート済みフラグ
            V[cursor] = true;
            an++;
            // 対象のコスト
            int v = A[cursor];
            // 最小コスト
            m = min(m,v);
            // 最小の値においての交換を加味しないサイクルの合計コスト
            S += v;
            // 対象が置かれるべき位置を割り出すT[v]に置かれるべき位置を入れる
            cursor = T[v];
            // サイクルが終了していたら、次に移る
            if(V[cursor]) break;
        }
        // サイクル内の最小値を利用したコスト計算(Σwi + (n - 2) * min(wi)
        int cost1 = S + (an - 2) * m;
        // 数列の最小コストを借りた場合のコスト(Σwi + min(wi) + (n + 1) * min(全数列)
        int cost2 = S + m + (an + 1) * s;
        // 最小値のコストを外から借りた場合とコストを比較する
        ans += min(cost1,cost2);
    }
    return ans;
}

int main(){
    cin >> n;
    // 最小値
    s = VMAX;
    for(int i = 0;i < n;++i){
        cin >> A[i];
        s = min(s,A[i]);
    }
    int ans = solve();
    cout << ans << endl;
}

参考

プログラミングコンテスト攻略のためのアルゴリズムとデータ構造

c++で8パズルを単純な幅優先探索で解く方法

aizu online judgeから問題を引用します。

8パズルの問題

8 パズルは１つの空白を含む 3×3のマス上に 8 枚のパネルが配置され、
空白を使ってパネルを上下左右にスライドさせ、絵柄を揃えるパズルです。

この問題では、次のように空白を 0、各パネルを 1 から 8 の番号でパズルを表します。

1 3 0
4 2 5
7 8 6

1 回の操作で空白の方向に１つのパネルを移動することができ、ゴールは次のようなパネルの配置とします。

1 2 3
4 5 6
7 8 0

８パズルの初期状態が与えられるので、ゴールまでの最短手数を求めるプログラムを作成してください。

8パズルの解法

状態の遷移状態がそれほど多くないことから、単純な深さ優先探索や幅優先探索で解くことができる。
これらの解法の練習になります。

パズルの状態を管理する構造体・クラスを作る

諸々必要な変数を定義します

// 要素数
#define N 3
#define N2 N * N

// 上,右,下,左の順でチェックする
static const int dx[4] = {0,1,0,-1};
static const int dy[4] = {-1,0,1,0};
static const char dir[4] = {'u','r','d','l'};

struct Puzzle{
    // 各々のタイルの状態(spaceは0だけど、9に換算する)
    int tiles[N2];
    // 空白の位置を管理indexで0 - 8まで
    int space;
    // 解法までのルートを記録
    string path;
    
    // 同じ状態遷移か判定する
    bool operator < (const Puzzle &p) const {
        for(int i = 0;i < N2;++i){
            if(tiles[i] == p.tiles[i]) continue;
            // 最短を求めるためにタイル数値が大きいと一致していると判定する(上から1,2,3と並ぶので)
            return tiles[i] > p.tiles[i];
        }
        return false;
    }
};

ポイント

タイルの表現方法

タイルは配列で管理します。
indexを割り振ると

1[0] 2[1] 3[2]
4[3] 5[4] 6[5]
7[6] 8[7] s[8]

となります。

枝刈りの方法

mapを使用して、その状態を保存するかどうか判定します。
最短経路を求める必要があるので、保存したマップの中でindexが小さいタイルの数値が大きい場合は、
非効率とみなして、一致していると判定するようにします。

パズルが一致しているかチェックする関数

indexで0から保存しているので、+1してチェックしています

/**
    パズルの一致チェック
*/
bool isCorrect(Puzzle p){
    // パズルが順番通り一致しているかチェック
    for(int i = 0;i < N2;++i){
        // 1,2,3,4,5,6,7,8の順で並ぶ
        if(p.tiles[i] != (i + 1)){
          return false;
        }
    }
    return true;
}

幅優先探索を行う関数

// 幅優先探索を行う
string bfs(Puzzle p){
    // 幅優先なのでqueue
    queue Q;
    // 同じ遷移状態がないか判定するためにmapを利用
    map V;
    Puzzle u,v;
    p.path = "";
    Q.push(p);
    // 最初の状態遷移を保存
    V[p] = true;
    
    // 幅優先探索
    while(!Q.empty()){
        u = Q.front();
        Q.pop();
        
        // 正解チェック
        if(isCorrect(u)){
            return u.path;
        }
        // 空白の位置を算出する
        int sx = u.space % N;
        int sy = u.space / N;
        
        // スペースの周りの四方のタイルを動かす
        for(int r = 0;r < 4;++r){
            // 動かすタイルのindexを算出する(上,右,下,左の順)
            int tx = sx + dx[r];
            int ty = sy + dy[r];
            // 範囲外チェック
            if(tx < 0 || ty < 0 || tx >= N || ty >= N){
                continue;
            }
            v = u;
            // 入れ替え
            swap(v.tiles[u.space],v.tiles[ty * N + tx]);
            // 入れ替えた場所がspaceになる
            v.space = (ty * N) + tx;
            
            // 1.同じ状態がでなかったら
            if(!V[v]){
                V[v] = true;
                v.path += dir[r];
                Q.push(v);
            }
        }
    }
    return "unsolvable";
}

ポイント

ポイントは、1のmapを使って、最短の経路を探しているところです。
探索するか探索しないかの状態チェックを行うのに先ほどのPuzzle構造体で定義したoperator関数を使っています。

8マップのように探索量がそれほどでもないなら、これでいけるようです。

全体のソースコード

最後に全体のソースコードを載せます

#include <iostream>
#include <cmath>
#include <string>
#include <map>
#include <queue>

using namespace std;

// 要素数
#define N 3
#define N2 N * N

// 上,右,下,左の順でチェックする
static const int dx[4] = {0,1,0,-1};
static const int dy[4] = {-1,0,1,0};
static const char dir[4] = {'u','r','d','l'};

struct Puzzle{
    // 各々のタイルの状態(spaceは0だけど、9に換算する)
    int tiles[N2];
    // 空白の位置を管理indexで0 - 8まで
    int space;
    // 解法までのルートを記録
    string path;
    
    // 同じ状態遷移か判定する
    bool operator < (const Puzzle &p) const {
        for(int i = 0;i < N2;++i){
            if(tiles[i] == p.tiles[i]) continue;
            // 最短を求めるためにタイル数値が大きいと一致していると判定する(上から1,2,3と並ぶので)
            return tiles[i] > p.tiles[i];
        }
        return false;
    }
};

/**
    パズルの一致チェック
*/
bool isCorrect(Puzzle p){
    // パズルが順番通り一致しているかチェック
    for(int i = 0;i < N2;++i){
        // 1,2,3,4,5,6,7,8の順で並ぶ
        if(p.tiles[i] != (i + 1)){
          return false;
        }
    }
    return true;
}

// 幅優先探索を行う
string bfs(Puzzle p){
    // queueを使って、パズルの状態遷移を管理する
    queue Q;
    // 同じ遷移状態がないか判定するためにmapを定義する
    map V;
    Puzzle u,v;
    p.path = "";
    Q.push(p);
    // 最初の状態遷移を保存
    V[p] = true;
    
    // 幅優先探索
    while(!Q.empty()){
        u = Q.front();
        Q.pop();
        
        // 正解チェック
        if(isCorrect(u)){
            return u.path;
        }
        // 空白の位置を算出する
        int sx = u.space % N;
        int sy = u.space / N;
        
        // スペースの周りの四方のタイルを動かす
        for(int r = 0;r < 4;++r){
            // 動かすタイルのindexを算出する(上,右,下,左の順)
            int tx = sx + dx[r];
            int ty = sy + dy[r];
            // 範囲外チェック
            if(tx < 0 || ty < 0 || tx >= N || ty >= N){
                continue;
            }
            v = u;
            // 入れ替え
            swap(v.tiles[u.space],v.tiles[ty * N + tx]);
            // 入れ替えた場所がspaceになる
            v.space = (ty * N) + tx;
            
            // 同じ状態がでなかったら
            if(!V[v]){
                V[v] = true;
                v.path += dir[r];
                Q.push(v);
            }
        }
    }
    return "unsolvable";
}

int main(){
    Puzzle in;
    
    for(int i = 0;i < N2;++i){
        cin >> in.tiles[i];
        if(in.tiles[i] == 0){
            // 空白は9として管理
            in.tiles[i] = N2;
            in.space = i;
        }
    }
    
    string answer = bfs(in);
    cout << answer.size() << endl;
    
    return 0;
}