矩陣連乘

漆楚衡發表於2014-12-13

矩陣連乘

給定n個矩陣 {A1, A2, A3, ... , An} ,其中 AiAi+1 是可乘的,i = 1, 2, ... , n - 1,考察這n個矩陣的連乘積A1A2...An

由於矩陣乘法滿足結合律,故計算矩陣的連乘積可以有多種不同的計算次序。這種計算次序可以用加括號的次序來確定。若一個矩陣連乘積的計算次序完全確定,也就是說該連乘積已完全加括號,則可以此次序反覆呼叫兩個矩陣相乘的標準演算法計算出矩陣連乘積。完全加括號的矩陣連乘積可遞迴地定義如下:

  1. 單個矩陣是完全加括號的
  2. 矩陣連乘積A是完全加括號的,則A可表示為2個完全加括號的矩陣連乘積B和C的乘積並加括號,即A = (BC)

每一種完全加括號方式對應於一種矩陣連乘積的計算次序,而矩陣連乘積的計算次序與其計算量有密切關係。

因為單次矩陣乘法A*B的複雜度取決於A的高度A.h,A的寬度(B的高度)A.w,B的寬度B.w(O(A.h * A.w * B.w))。不同的加括號次序可以產生不同的計算量。

矩陣連乘積的最有計算次序問題:對於給定的相繼n個矩陣{A1, A2, ... , An}(其中,Ai的維數為Pi-1*Pi, i = 1, 2, ... , n)如何確定計算矩陣連乘積A1A2...An的計算次序(完全夾括方式),使得依次計算次序計算矩陣連乘積所需要的數乘次數最少。

問題分解(遞迴思路)

根據題目之前對於加括號的定義(一個遞迴結構)。我們可以構造出這個問題的分解方法,對於給定矩陣串As(長度為n,A0, ... An - 1),如果n > 1,我們可以選擇從第Ai(i = 1 ... n - 1)個矩陣前斷開矩陣連得到兩個子矩陣串As1 = As[0, i),As2 = As[i, n)。可以看出它們具有與原問題相同的形式。

於是,設T(As)為As連乘積的最少加括號數乘量,有

  • T(As) = 0 n = 1(單個矩陣)
  • T(As) = min(T(As[0, i)) + T(As[i, n)) + As[0].h * As[i].h * As[n - 1].w) (i = 1...n - 1) |n > 1|

最優子結構

注意到我們在上面已經利用到了這個問題的最優子結構:我們在解決As的最優解時分解了As,通過一系列T(As1),T(As2)我們就判斷出了As的最優解,而沒有考慮是否As1或As2的非最優解構造As的最優解。

我們知道對於一個矩陣串As,對於其最終結果M,始終有M.h = As[0].h, M.w = As[n - 1].w,這不受加括號次序的影響:對於As[0]參與的任意相乘過程,As[0].h都會作為新矩陣的h保留下來,並新矩陣依然處於矩陣串頭部,As[n - 1].w同理。

對於As給定斷開位置i得到的子解(As1, As2)。As1(As2)的任意加括號次序所組成的M1(M2)都有相同的h和w,考察上面的遞迴式,As[0].h * As[i].h * As[n - 1].w可以寫作M1.h * M2.h * M2.w。即這一部分獨立於As1或As2的加括號方式。

對於As給定斷開位置i,T(As)由三部分組成:

  • As1的某種數乘次數a
  • As2的某種數乘次數b
  • 合併成本M1*M2為常數

顯然,要使T(As)最小,就是要a + b最小,即使a = T(As1), b = T(As2)

分析問題(子問題重疊)

通過剛才的分析,我們已經能夠給出一個遞迴的演算法來解決這個問題。現在我們來考慮優化的可能性。

對於As,我們設x = [l,r)是一個某子問題的區間。

  • 如果我們取斷開位置i >= r或i <= l,則有x出現在i的每一個取值的子問題中並被重複計算一次
  • 如果有l = 0r = n,則x會被父問題求解一次,並有可能在另一種分割中被傳遞給子問題重複求解

不難想象,對於n越大,x越小(r - l越小),重複強度會急劇增大。

引入記憶

本題中,一個問題可以表現為一個區間([l, r)),所以我們可以用二維陣列表示記憶表。

int mem[l][r] = {Inf};    //存放T(As[l, r)),初值可以取任意大於可能解的數。

問題初始化部分:

const int n = init();

struct{
    int h;
    int w;
} As[N] = init();

不難寫出帶記憶表的遞迴搜尋演算法:

int T(int l, int r){
    if(mem[l][r] == Inf){
        for(int i = l + 1; i < r; i++){
            mem[l][r] = min(mem[l][r]
                , T(l, i) + T(i, r) + As[l].h * As[i].h * As[r - 1].w);
        }
    }
    return mem[l][r];
}

動態規劃

現在自底向上構建解

  • 問題的基本是當子問題規模為1,此時解為0
  • 當長度小於len的所有子問題都解決了,就有了足夠的資訊來構建長度為len的解

int TDP(){
    for(int i = 0; i < n; ++i){
        mem[i][i + 1] = 0;
    }
    for(int len = 2; len <= n; ++len){
        for(int l = 0; l < n; l++){
            int r = l + len;
            for(int i = l; i < r; ++i){
                mem[l][r] = min(mem[l][r]
                    , mem[l][i] + mem[i][r] + As[l].h * As[i].h * As[r - 1].w);
            }
        }
    }
    return mem[0][n];
} 

相關文章