15. Dynamic Programming, 15.2 連鎖行列積

# 15.2 連鎖行列積

n個の行列の列として<A_1, A_2, A_3, ..., A_n>が与えられたとき、積 $$ A_1 A_2 ... A_n $$ を計算することを考える。

行列の積は結合的なので計算順序によらず解は同じになる。
計算順序の曖昧さは括弧で解消する。

```math
積 A_1 A_2 A_3 A_4 を完全に括弧付けする5通りの方法 \\
(A_1 (A_2 (A_3 A_4))) \\ 
(A_1 ((A_2 A_3) A_4)) \\ 
((A_1 A_2) (A_3 A_4)) \\ 
((A_1 (A_2 A_3)) A_4) \\ 
(((A_1 A_2) A_3) A_4) \\
がある
```

連鎖行列に対する括弧付けの仕方が積の計算コストに劇的な影響を与えることがある。まず、2つの行列の積の計算コストを考えてみよう。

```
MATRIX-MULTIPLY(A, B)
if A.columns != B.rows
  error "error"
else Cを新しいA.rows * B.columns 型行列とする
  for i = 1 to A.rows
    for j = 1 to B.columns
      c_ij = 0
      for k = 1 to A.columns
        c_ij = c_ij + a_ik * b_kj
  return C
```

2つの行列AとBを掛けることができるのは、それらが両立可能な行列の場合、すなわち A の列数が B の行数に等しいときだけである。 A が p * q 型の行列で B が q * r 型の行列ならば、結果 C は p * r 型の行列である。
計算時間を支配するのは`c_ij = c_ij + a_ik * b_kj`の回数pqrである。

以下では実行時間をスカラー乗算回数で評価する。

ref: [Infographics: Operation Costs in CPU Clock Cycles](http://ithare.com/infographics-operation-costs-in-cpu-clock-cycles/)

## 括弧付けの仕方によって連鎖行列の計算コストが変化する例

サイズが3の連鎖行列 <A_1, A_2, A_3>, A_1:10x100, A_2:100x5, A_3:5x50 とする。

((A_1 A_2) A_3) => A_1 A_2 の計算に 10x100x5 = 5000回、結果は B(サイズ10x5) とする。B A_3 の計算に 10x5x50 = 2500回。総計算回数7500。

(A_1 (A_2 A_3)) => A_2 A_3 の計算に 100x5x50 = 25000回、結果は C(サイズ100x50) とする。A_1 C の計算に10x100x50 = 50000回。総計算回数75000。

10倍計算回数に差が出る。

## 連鎖行列積問題

n個の行列の連鎖 <A_1, A_2, ... , A_n> が与えられたとき、スカラー乗算回数を最小化するに括弧付けする問題である。ただし、i = 1, 2, ..., n に対して行列 A_i を p_{i-1} \times p_i 型とする。
連鎖行列積問題は行列積の計算結果を求めるわけではないことに注意せよ。最小コストを達成する行列積の計算順序を決定することが目的である。

## 異なる括弧付けの個数

連鎖行列積問題を動的計画法で解く前に総当りでは効率のよいアルゴリズムを得ることができないことを確認しておこう。n個の行列の列に対する異なる括弧付の個数を P(n) によって表す。n = 1の場合には行列が1つしかないので、この行列積を完全に括弧付ける方法は1つだけである。n >= 2の場合には、完全な括弧付は2つの完全な括弧付の積であり、任意のk = 1, 2, ..., n-1 に対して、k番目とk+1番目の行列の間で2つの部分積へ分割できる。したがって、漸化式

```math
P(n) = \left\{
\begin{array}{ll}
1 & (n = 1)\\
\sum_{k=1}^{n-1}{P(k)P(n-k)} & (n \geq 2)
\end{array}
\right.
```

を得る。章末問題12-4で類似する漸化式の解がカラタン数列になることを証明するが、この数列は$$ \Omega(4^n/n^{3/2}) $$の速さで増大する。漸化式の解が$$ \Omega(2^n) $$であることの証明は練習問題15.2-3で。したがって総当たりは適さない。

## 動的計画法の適用

1. 最適解の構造を特徴づける
2. 最適解の値を再帰的に定義する
3. 最適解の値を計算する
4. 計算された情報から1つの最適解を構成する

## 第1段階: 最適括弧付けの構造

i <= j の時、積 A_i A_{i+1} ... A_j を計算した結果である A_{i..j} と記す。

i < j ならば、積 A_i A_{i+1} ... A_j のどの完全な括弧付けも、ある整数 k(i <= k < j) によって積を A_k と A_{k+1} の間で分割する。すなわち、ある k の値に対して、A_{i..k} と A_{k+1..j}をまず計算した後、それらの積を取って最終的な積 A_{i..j} を得る。したがって、

```math
(A_i A_{i+1} ... A_k)(A_{k+1} ... A_j) を計算するコスト = \\
(A_i A_{i+1} ... A_k) を計算するコスト \\
+ (A_{k+1} ... A_j) を計算するコスト \\
+ A_{i..k} A_{k+1..j} の積を計算するコスト
```

この問題が持つ部分構造最適性を説明する。積 A_i A_{i+1} ... A_j の最適括弧付けは、この積を A_k と A_{k+1} の間で分割すると仮定する。このとき、 A_i A_{i+1} ... A_j に対する最適括弧付けの前半の部分連鎖行列 A_i A_{i+1} ... A_k に対する部分括弧付は A_i A_{i+1} ... A_k に対する最適括弧付けでなければならない。同様に A_k A_{k+1} ... A_j についても。

問題に対する最適解を、その部分問題に対する最適解から構成できることを、上で述べた部分構造最適性を用いて証明しよう。連鎖行列積問題の自明でないインスタンスに対する任意の解は積を2つに分割し、任意の最適解はその中に部分問題インスタンスの最適解を含む。

## 第2段階: 再帰的な解

1 <= i <= j <= n に対して、 A_i A_{i+1} ... A_j の最小コスト括弧付を決定する問題を部分問題として取り上げる。

m[i, j] を行列 A_{i..j} の計算に最小限必要なスカラー乗算回数とする。元の問題である A_{1..n} を最も効率よく計算するときのコストは定義より m[1, n] である。

m[i, j] を再帰的に定義する。 i = j ならば 0。第一段階の考察した部分構造最適性を用いる。 i < j のとき、i <= k < j となるkで、最適な括弧付けは A_k と A_{k+1} の間で分割すると仮定する。各行列 A_i は p_{i-1} \times p_i 型だから、行列積 A_{i..k} A_{k+1..j} は p_{i-1} p_k p_j 回のスカラー乗算で計算できる。したがって

```math
m[i, j] = m[i, k] + m[k+1, j] + p_{i-1} p_k p_j 
```

が成立する。この再帰方程式では未知の値 k を既知である仮定している。しかし、k の取りうる値は j - 1 通りあり、最適括弧付けはこれら全てを調べて最良のものを求めればよい。したがって、積 A_i A_{i+1} ... A_j の計算に最小限必要なコストの再帰的定義は

```math
m[i, j] = \left\{
\begin{array}{ll}
0 & (i = j)\\
\min_{i \leq k < j} \{m[i, k] + m[k + 1, j] + p_{i-1} p_k p_j\} & (i < j)
\end{array}
\right.
```

となる。連鎖行列積問題は最適な計算順序をを決定する問題なので m[i, j] だけでは不十分で、分割地点 k の値(s[i, j]として定義する)が必要である。

## 第3段階: 最適コストの計算

再帰で素直に組むとロッド切り出し問題と同様に指数時間かかってしまう。

部分問題が比較的少ないことに気づくことが重要。 1 <= i <= j <= n を満たす i と j の組の数に1つの部分問題が対応し、全体で 

```math
2Cn + 1Cn = \Theta (n^2)
```

個しかない。再帰アルゴリズムではその再帰🎄の様々な分岐において同じ部分問題に何度も出会うことがある。この部分問題重複制がこの問題に動的計画法を適用できる可能性の高さを示す第2の特徴である(第1の特徴は部分構造最適性である)。

```
MATRIX-CHAIN-ORDER(p)
n = p.length - 1
m[1..n, 1..n] と s[1..n-1, 2..n] を新しい表とする
for i = 1 to n
  m[i, i] = 0
for l = 2 to n // l は連鎖の長さ
  for i = 1 to n - l + 1
    j = i + l - 1
    m[i, j] = \infty
    for k = i to j - 1
      q = m[i, k] + m[k+1, j] + p_{i-1} p_k p_j
      if q < m[i, j]
        m[i, j] = q
        s[i, j] = k
return m and s
```

3重ループなのでアルゴリズムの実行時間はO(n^3)。実行時間が \Omega (n^3) でもあることを練習問題15.2-5で証明する。表mと表sを格納するために \Theta (n^2) 領域が必要である。

実際の例は教科書参照

1..6 => k = 3
A_1 A_2 A_3 A_4 A_5 A_6 = (A_1 A_2 A_3) (A_4 A_5 A_6)
1..3 => k = 1, 4..6 => k = 5
(A_1 (A_2 A_3))((A_4 A_5) A_6)

## 第4段階: 最適解の構成

```
PRINT-OPTIMAL-PARENS(s, i, j)
if i == j
  print "A"i
else print "("
  PRINT-OPTIMAL-PARENS(s, i, s[i, j])
  PRINT-OPTIMAL_PARENS(s, s[i, j] + 1, j)
  print ")"
```

## 練習問題

### 15.2-1

与えられた次元列の連鎖行列積に対する最適括弧付けを求める

```scala
object Main extends App {
  def printMatrix(matrix: Array[Array[Int]]): Unit = println(matrix.map(array => array.mkString(",")).mkString("\n"))

  val p = Array(5, 10, 3, 12, 5, 50, 6)
  val n = p.length - 1
  val m = Array.ofDim[Int](n, n)
  val s = Array.ofDim[Int](n, n)
  for (i <- Range(0, n)) m(i)(i) = 0
  for (l <- Range(2, n + 1)) {
    for (i <- Range(0, n - l + 1)) {
      val j = i + l - 1
      m(i)(j) = Int.MaxValue
      for (k <- Range(i, j)) {
        val q = m(i)(k) + m(k + 1)(j) + p(i) * p(k + 1) * p(j + 1)
        if(q < m(i)(j)) {
          m(i)(j) = q
          s(i)(j) = k
        }
      }
    }
  }
  printMatrix(m)
  printMatrix(s)
}
```

<5, 10, 3, 12, 5, 50, 6>

  j\i |    0 |    1 |    2 |    3 |    4 |    5 |
 :--: | :--: | :--: | :--: | :--: | :--: | :--: |
    0 |    0 |    - |    - |    - |    - | - | 
    1 |  150 |    0 |    - |    - |    - | - | 
    2 |  330 |  360 |    0 |    - |    - | - | 
    3 |  405 |  330 |  180 |    0 |    - | - | 
    4 | 1655 | 2430 |  930 | 3000 |    0 | - | 
    5 | 2010 | 1950 | 1770 | 1860 | 1500 | 0 | 

  j\i |    0 |    1 |    2 |    3 |    4 |    5 |
 :--: | :--: | :--: | :--: | :--: | :--: | :--: |
    0 |    - |    - |    - |    - |    - | - | 
    1 |    0 |    - |    - |    - |    - | - | 
    2 |    1 |    1 |    - |    - |    - | - | 
    3 |    1 |    1 |    2 |    - |    - | - | 
    4 |    3 |    1 |    2 |    3 |    - | - | 
    5 |    1 |    1 |    3 |    3 |    4 | - | 

実際の分割方法
1. A_0 A_1 A_2 A_3 A_4 A_5
2. i:0, j:5 => k = 1 => (A_0 A_1)(A_2 A_3 A_4 A_5)
3. i:2, j:5 => k = 3 => (A_0 A_1)((A_2 A_3)(A_4 A_5))

実際の計算回数
1. A_0 A_1 => 150
2. A_2 A_3 => 180
3. A_4 A_5 => 1500
4. A_{2..3} A_{4..5} => 3 * 5 * 6 = 90
5. A_{0..1} A_{2..5} => 5 * 3 * 6 = 90
6. 総コスト = 150 + 180 + 1500 + 90 + 90 = 2010

### 15.2-2

行列の列とMatrixChainOrderが計算した表sと添字i,jが与えられたとき、実際に最適な順序で連鎖行列積の計算を行う再帰アルゴリズムを設計せよ。

```scala
object Main extends App {
  def matrixChainMultiply(matrixArray: Array[Array[Array[Int]]], s: Array[Array[Int]], i: Int, j: Int): Array[Array[Int]] = {
    if(i == j) matrixArray(i)
    else {
      val k = s(i)(j)
      matrixChainMultiply(matrixArray, s, i, k) * maxtrixChainMultiplay(matrixArray, s, k + 1, j)
    }
  }

  val matrixArray = Array.ofDim[Int](5, 8, 8)
  val s = Array.ofDim[Int](8, 8)
  matrixChainMultiplay(matrixArray, s, 0, 5 - 1)
}
```

### 15.2-3

漸化式 15.6

```math
P(n) = \left\{
\begin{array}{ll}
1 & (n = 1)\\
\sum_{k=1}^{n-1}{P(k)P(n-k)} & (n \geq 2)
\end{array}
\right.
```

の解が

```math
\Omega (2^n)
```

であることを置換え法を用いて示せ。

1. P(1) = 1
2. P(2) = 1
3. P(3) = P(1) P(2) + P(2) P(1) = 2
4. P(4) = P(1) P(3) + P(2) P(2) + P(3) P(1) = 5
5. P(5) = P(1) P(4) + P(2) P(3) + P(3) P(2) + P(4) P(1) = 14
6. P(6) = P(1) P(5) + P(2) P(4) + P(3) P(3) + P(4) P(2) + P(5) P(1) = 42
7. P(7) = ... = 132 > 2^7
8. ...
9. P(n) = P(1) P(n-1) + ... + P(n-1) P(1) >= 2 * P(n-1)

P(n) は n >= 7 ならば 2^n より大きい。
n >= 3 ならば P(n) は 2 * P(n-1) より大きい。
以上より
n >= 7 ならば P(n) は 2^n より大きい。

### 15.2-4

連鎖行列積問題の部分問題グラフを描け。

ホワイトボードでがんばろう。

```math
頂点数 = \sum_{i=1}^n (n - i) = \frac{n(n + 1)}{2} \\
```

### 15.2-5

MATRIX-CHAIN-ORDER の呼び出しの中で、表の要素 m[i, j] が他の要素のために参照される回数を R(i, j) で表す。表全体での参照回数の総和は

```math
\sum_{i=1}^n \sum_{j=i}^n R(i,j) = \frac{n^3-n}{3}
```

であることを示せ。

MATRIX-CHAIN-ORDER の定義より

```math
\sum_{i=1}^n \sum_{j=i}^n R(i,j) = \sum_{l=2}^n \sum_{i=1}^{n-l+1} \sum_{k=i}^{j-1} 2 \\
= \sum_{l=2}^n \sum_{i=1}^{n-l+1} \sum_{k=i}^{i+l-2} 2 \\
= \sum_{l=2}^n \sum_{i=1}^{n-l+1} 2(l-1) \\
= \sum_{l=2}^n 2(l-1)(n-l+1) \\
= \sum_{l=1}^{n-1} 2 l (n-l) = \sum_{l=1}^{n-1} (2 n l - 2 l^2) \\
= 2 n \sum_{l=1}^{n-1} l - 2 \sum_{l=1}^{n-1} l^2 \\
= 2 n * \frac{(n-1)n}{2} - 2 * \frac{(n-1)n(2n-1)}{6} \\
= n^3 - n^2 - \frac{2n^3-3n^2+n}{3} \\
= \frac{n^3-n}{3}
```

### 15.2-6

要素数が n の式に対する完全な括弧付けはちょうど n - 1 個の括弧対を持つことを示せ

完全に括弧付けされているの定義は『それが単一の行列である、または、2つの完全に括弧付けされた行列の積が括弧で囲まれている』である。よって、明らか