動態(tài)規(guī)劃求解矩陣連乘的最優(yōu)時間復雜度

給定一系列矩陣A=(Ai)N?1i=0求他們的連續(xù)相乘結果 A=ΠN?1i=0Ai尋找最優(yōu)的相乘組合，使得計算 A 所需的時間復雜度最小。

由于矩陣乘法要求兩個相乘矩陣的維度滿足：第一個矩陣的列數(shù)要與第二個矩陣的行數(shù)相同。所以我們只要用N+1個數(shù)字就能表示所有矩陣的維度了，這里我們用 d 來表示這 N+1 個數(shù)字， 其中di 和 di+1 分別表示第 i 個矩陣的行數(shù)和列數(shù)。

給定一個矩陣序列 A， 我們并不需要真正計算矩陣乘法，而是給出最優(yōu)時間復雜度和矩陣相乘順序。因此，我們真正的輸入是d。這里我們暫且不考慮輸出矩陣相乘的順序，先以求最優(yōu)時間復雜度為目標解決這個優(yōu)化問題。

如果你不想看分析過程可以直接看后面的算法實現(xiàn)部分。

如果我們用C(?) 表示某個矩陣連乘序列的最優(yōu)時間復雜度，那么它一定滿足下面的公式：
C(ΠN?1i=0Ai)=min{C(Πkm=0Am)+C(ΠN?1n=k+1An)+d0?dk+1?dN}N?2k=0.(1)
類似地，對于任意的正整數(shù) u 和 v， 其中 0≤u≤v≤N?1， 我們有：
C(Πvi=uAi)=min{C(Πkm=uAm)+C(Πvn=k+1An)+du?dk+1?dv+1}v?1k=u.(2)
而且我們還知道：
C(Ai)=0,i=0,…,N?1.(3)
那么，如果我們用矩陣 B 中元素 b(u,v) 表示 C(Πvi=uAi)， 我們從公式 (2) 可以看出，b(u,v) 只依賴于 b(u,u:1:v?1) 和 b(u+1:1:v,v) ，其中u:1:v?1 表示從 u 以步長 1 增長到 v?1。另外，對角線上元素都為0。故此，可以借助一個二維數(shù)組來尋找u=0 ， v=N?1時的最優(yōu)時間復雜度，每一次找到的時間復雜度記錄下最優(yōu)的 k 值就可以知道如何劃分矩陣來相乘了。

利用動態(tài)規(guī)劃的思想解決矩陣序列連乘問題的算法本身的時間復雜度跟B 矩陣的計算有關，B 矩陣需要計算其整個上三角部分，我們逐步推導：

第 1 列: 計算B(0,1): 需要 1?0=1 次計算。

第 2 列: 計算B(1,2): 需要 2?1=1 次計算。
第 2 列: 計算B(0,2): 需要 2?0=2 次計算。

第 3 列: 計算B(2,3): 需要 3?2=1 次計算。
第 3 列: 計算B(1,3): 需要 3?1=2 次計算。
第 3 列: 計算B(0,3): 需要 3?0=3 次計算。

?

第 N?1 列: 計算B(N?2,N?1): 需要 (N?1)?(N?2)=1 次計算。
第 N?1 列: 計算B(N?3,N?1): 需要 (N?1)?(N?3)=2 次計算。
第 N?1 列: 計算B(N?4,N?1): 需要 (N?1)?(N?4)=3 次計算。
?
第 N?1 列: 計算B(0,N?1): 需要 (N?1)?(0)=N?1 次計算。

所以計算時間復雜度為：

O({1}+{1+2}+…+{1+2+…+(N?1)})=O(∑n=1N?1∑r=1nr)=O(∑n=1N?1(n+1)n2)=O(∑n=1N?1n2+∑n=1N?1n)=O(N3)
時間復雜度的推到請參考這個鏈接 https://en.wikipedia.org/wiki/Faulhaber%27s_formula

雖然算法時間復雜度為O(N3)， 我們只需要存儲一個矩陣就可以了，所以空間復雜度是 O(N2)。

完整的C++實現(xiàn)如下：

#include 
#include 
using namespace std;

// 尋找最優(yōu)時間復雜度 B，以及最優(yōu)劃分 K
void find_best_complexity(vector &B, vector &K, const int *d, int N){  

  B.resize(N*N);
  K.resize(N*N);

  for (int i = 0; i < N; i++){
    B[i*N + i] = 0;
  }

  for (int v = 1; v < N; v++){
    for (int u = v - 1; u > -1; u--){

      int best_cmp = INT_MAX;
      int best_k;
      for (int k = u; k < v; k++){
        int current_cmp = d[u] * d[k + 1] * d[v + 1] + B[u*N + k] + B[(k+1)*N + (v)];
        if (current_cmp < best_cmp){
          best_cmp = current_cmp;
          best_k = k;
        }
      }

      K[u*N + v] = best_k;
      B[u*N + v] = best_cmp;
    }
  }
}

// 輸出最優(yōu)時間復雜度下矩陣的相乘順序
void print_uv(int u, int v, vector &K, int &N){
  if (u==v){
    return;
  }
  int k = K[u*N+v];
  print_uv(u, k, K, N);
  print_uv(k + 1, v, K, N);

  printf("%4d  ", u);
  printf("%4d  ", k);
  printf("%4d  n", v);
}

// 舉個例子
int main(int argc, char** argv){

  int d[] = {1, 2, 3, 1, 5};
  int N = (sizeof(d) / sizeof(d[0])) - 1;

  vector B;
  vector K;
  find_best_complexity(B, K, d, N);


  printf("###############################################################n");
  printf("# Bn");
  for (int u = 0; u < N; u++){
    for (int v = 0; v < N; v++){
      if (v < u){
        printf("%4d  ", -1);
      }
      else{
        printf("%4d  ", B[u*N + v]);
      }
    }
    printf("n");
  }

  printf("###############################################################n");
  printf("# Kn");
  for (int u = 0; u < N; u++){
    for (int v = 0; v < N; v++){
      if (v < u){
        printf("%4d  ", -1);
      }
      else{
        printf("%4d  ", K[u*N + v]);
      }
    }
    printf("n");
  }

  printf("###############################################################n");
  printf("# ordern");
  print_uv(0, N - 1, K, N);

    return EXIT_SUCCESS;
}

上述代碼用printf是為了更好的格式化輸出。