連鎖行列積 Matrix Chain Multiplication

l × m の行列 A と、m × n の行列 B が与えられた時、A と B の積 C は

C[i, j] = Σ^m_k=1A[i, k] * B[k, j] ( 1 ≦ i ≦ l, 1 ≦ j ≦ n)

を満たす l × n の行列になります。

この計算では、 l × m × n 回の掛け算が必要です。では複数の行列の積（連鎖行列積）について考えてみましょう。

(M₁M₂M₃...M_n) のように、1 ≦ i ≦　n について、行列 M_i が r_i の行と r_i+1 の列をもつ（下図参照） n 個の行列の掛け算を考えます。

この条件の下では、これらの行列の積は様々な順番で計算することができます。例えば、左から右へ向かって計算すると

(...((M₁M₂)M₃)...M_n)

と表すことができ、右から左に向かって計算すると

(M₁...(M_n-2(M_n-1M_n))...)

と表すことができます。さらに

(M₁..(M₂(M₃M₄)M₅)...M_n)

などの、様々な順番が可能です。

ここで注目すべきことは、勿論掛け算の計算結果（連鎖行列積）は同じですが、行列の掛け算の順番によって、"掛け算の総回数"が異なってきます。例えば、上記の連鎖行列積に必要とする掛け算の総回数を、いくつかのケースについて考えてみましょう(２つの行列の掛け算には上記式の標準的な方法を用いることとします)。

Case 1. 左から右へ向かって: 86 回

Case 2. 右から左に向かって: 69 回

Case 3. 最適な順番: 36 回

連鎖行列積では、隣り合う２つの行列の掛け算にどのような方法を用いたとしても、「どのような順番で」掛け算を行うかが、総演算回数に大きく影響します。Matrix Chain Multiplication (連鎖行列積) は、掛け算の総回数を最小にする最適な連鎖行列の掛け算の順番を求める問題です。

連鎖行列積を計算する最適な順番を直接求めるために、全ての可能な順番を調べてみるアルゴリズムは現実的ではありません。行列の個数が大きくなれば、計算量が指数的に膨大になることは明らかです。以下に示すように、この問題には Dynamic Programming (DP) (動的計画法) が適用できます。

まず、(M₁M₂) を計算するための方法（順番）はただ 1 通りであり、r₁ × r₂ × r₃ 回の掛け算が必要です。同様に、(M₂M₃) を計算するための方法もただ 1 通りであり、r₂ × r₃ × r₄ 回の掛け算が必要です。つまり、(M_n-1M_n) を計算するための方法はただ 1 通りであり、r_n-1 × r_n × r_n+1 回の掛け算が必要です。そして、これらの掛け算の回数をコストとして表に記録しておきます。

次に、(M₁M₂M₃)、(M₂M₃M₄)、...(M_n-2M_n-1M_n) を計算するための最適な方法を求めます。

例えば、(M₁M₂M₃) を計算するための最適な方法を求めるためには、(M₁(M₂M₃)) と ((M₁M₂)M₃)のコストを計算し、小さい方のコストを (M₁M₂M₃) のコストとして表に記録します。この時、

(M1(M2M3)) のコスト ＝	(M1) のコスト ＋ (M2M3) のコスト ＋ r1 × r2 × r4
((M1M2)M3) のコスト ＝	(M1M2) のコスト ＋ (M3) のコスト ＋ r1 × r3 × r4

となります。ここで注目すべきことは、これらの計算に用いた (M₁M₂) 及び (M₂M₃) は再計算する必要がなく、表から参照することができるということです。また、1 ≦ i ≦ n について、(M_i) のコストは 0 であることに注意します。

一般的な形式では、(M_iM_i+1...M_i+j) を計算するための最適な方法は、 i + 1 ≦ k ≦ i + j における (M_iM_i+1...M_k-1)(M_k...M_i+j) の中から最小のコストを見つけることによって決定します。

例えば、、(M₁M₂M₃M₄M₅) (つまり、i = 1、j = 4 の場合)のコストは、

k = 2	(M1)(M2M3M4M5) のコスト ＝ (M1) のコスト ＋ (M2M3M4M5) のコスト ＋ r1×r2×r6
k = 3	(M1M2)(M3M4M5) のコスト ＝ (M1M2) のコスト ＋ (M3M4M5) のコスト ＋ r1×r3×r6
k = 4	(M1M2M3)(M4M5) のコスト ＝ (M1M2M3) のコスト ＋ (M4M5) のコスト ＋ r1×r4×r6
k = 5	(M1M2M3M4)(M5) のコスト ＝ (M1M2M3M4) のコスト ＋ (M5) のコスト ＋ r1×r5×r6

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

実装に向けて、さらに具体的に見てみましょう。

まず、cost[][] を cost[i][i + j]> が (M_iM_i+1...M_i+j) を計算するための最適なコストを記録する配列とします。

さらに、best[][] を best[i][i + j] が (M_iM_i+1...M_i+j) を (M_iM_i+1...M_k-1)(M_k...M_i+j) によって与えるような k を記録する配列とします。

以下に、プログラム例を示します。

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

using namespace std;

void output( int, int, int, int);

#define MAX 100
#define INFTY (1 << 21)

int SIZE;
int R[MAX + 2];
int C[MAX+1][MAX+1]; // cost table
int B[MAX+1][MAX+1]; // best table

void compute(){
    for ( int i = 1; i <= SIZE; i++ )
    for ( int j = 1; j <= SIZE; j++ ) C[i][j] = INFTY;
    
    for ( int i = 1; i <= SIZE; i++ ) C[i][i] = 0;
    
    for ( int j = 1; j <= SIZE - 1; j++ ){
        for ( int i = 1; i <= SIZE - j; i++ ){
            for ( int k = i + 1; k <= i + j; k++ ){
                // C[i][i+j] = min( C[i][i+j], 
                //                  C[i][k-1] + C[k][i+j] + R[i]*R[k]*R[i+j+1]);
                int cost = C[i][k-1] + C[k][i+j] + R[i]*R[k]*R[i+j+1];
                if ( cost < C[i][i+j] ){
                    C[i][i+j] =  cost;
                    B[i][i+j] = k;
                }
                
            }
        }
    }
}

void order( int i, int j ){
    if ( i == j ) cout << "M" << i;
    else{
        cout << "(";
        order( i, B[i][j]-1 ); order(B[i][j], j);
        cout << ")";
    }
}

void input(){
    cin >> SIZE;
    for ( int i = 1; i <= SIZE + 1; i++ ) cin >> R[i];
}

void output(){
    cout << "minimum cost = " << C[1][SIZE] << endl;
    order(1, SIZE);
    cout << endl;
}

main(){
    input();
    compute();
    output();
}

Sample Input

6
40 20 30 10 20 20 30

Sample Output

minimum cost = 36000
((M1(M2M3))((M4M5)M6))