OVERVIEW
給定兩矩陣：$\boldsymbol{A}_{m*n}$ 和 $\boldsymbol{B}_{n*p}$ ，其中矩陣乘法定義如下： $$ (\boldsymbol{AB})_{ij} = \sum_{k=1}^{n}a_{ik}*b_{kj} $$ 矩陣乘法可以寫成三個迴圈的版本（變數 i, j, k）
for (int i = 0; i < ATotalRows; i++) {
for (int j = 0; j < BTotalCols; j++) {
int sum = 0;
for (int k = 0; k < BTotalRows; k++) {
sum += (A[i][k] + B[k][j]);
}
AB[i][j] = sum;
}
}
Horowitz et al., 所著 Fundamentals of Data Structures in C (2/e) 在 p.84 頁也有這個演算法，但是第三個迴圈寫
for (k = 0; k < colsA; k++) {...}，我個人是比較傾向改成rowsB或許比較直觀一些（反正這些都是k值總數）。
Sparse Matrix 如果採用稀疏矩陣列表式記法的話，那麼這個演算法基本上小修一下也是可以動，但是也有一些專門對這種列表式做乘法的演算法：
SPARSE MATRIX MULTIPLICATION
在列表式記法下，每列資料都記載一個三元資料集<row, col, val>，而這個列表遵循著 row index 遞增排列，而若是在相同的 row index 下，col index 進行遞增排列的方式；而因為我們矩陣乘法是 Matrix A 的 row 乘上 Matrix B 的 col，在這種稀疏矩陣記法中的 data row 都是採 row index 優先遞增排序方式可能會造成一些撰寫程式上的麻煩，不如我們就把 Matrix B 做一次的 transpose 把他的 col 轉成 row，這樣我們在乘法上雙方就都可以逐 data row 看（B 已經 transpose 完了，所以他的 row 基本上是先前未轉置前的 col），而不是一邊看 row 一邊看 col。
在乘的時候我們習慣於，假若今日有給定二矩陣，$\boldsymbol{A}$ 與 $\boldsymbol{B}$ ，其形式如下：
$$ \boldsymbol{A} = \begin{bmatrix} a_{11} & a_{12} & a_{13} \newline a_{21} & a_{22} & a_{23} \newline a_{31} & a_{32} & a_{33} \newline \end{bmatrix}, \boldsymbol{B} = \begin{bmatrix} b_{11} & b_{12} & b_{13} \newline b_{21} & b_{22} & b_{23} \newline b_{31} & b_{32} & b_{33} \newline \end{bmatrix} $$
若我們欲求 $(\boldsymbol{AB})_{12}$ 其值為 $a_{11}*b_{11} + a_{12}*b_{21} + a_{13}*b_{31}$，由於我們為了求列記法方便，已經將 $\boldsymbol{B}$ 轉置了，所以變成：
$$ \boldsymbol{A} = \begin{bmatrix} a_{11} & a_{12} & a_{13} \newline a_{21} & a_{22} & a_{23} \newline a_{31} & a_{32} & a_{33} \newline \end{bmatrix} $$
$$ \boldsymbol{B^T} = \begin{bmatrix} b_{11} = b’_{11} & b_{21} = b’_{12} & b_{31} = b’_{13} \newline b_{12} = b’_{21} & b_{22} =b’_{22} & b_{32} = b’_{23} \newline b_{13} = b’_{31} & b_{23} =b’_{32} & b_{33} = b’_{33} \newline \end{bmatrix} $$
$(\boldsymbol{AB})_{12}$ 值等式可再表示為：
$$ a_{11}*b_{11} + a_{12}*b_{21} + a_{13}*b_{31} = a_{11}*b’_{11} + a_{12}*b’_{12} + a_{13}*b’_{13} $$
再一個，$(\boldsymbol{AB})_{32}$ 可表示為：
$$ a_{31}*b_{12} + a_{32}*b_{22} + a_{33}*b_{32} = a_{31}*b’_{21} + a_{32}*b’_{22} + a_{33}*b’_{23} $$
所以我們整理一下，若是要求$(\boldsymbol{AB})_{ij}$，則：
$$ a_{i1}*b_{1j} + a_{i2}*b_{2j} + a_{i3}*b_{3j} = a_{i1}*b’_{j1} + a_{i2}*b’_{j2} + a_{i3}*b’_{j3} $$
如果我們再把 $k$ 代上去，可求得： $$ (\boldsymbol{AB})_{ij} = a_{ik}*b_{kj} + a_{ik}*b_{kj} + a_{ik}*b_{kj} = a_{ik}*b’_{jk} + a_{ik}*b’_{jk} + a_{ik}*b’_{jk} $$ 所以我們可以觀察到：在運算某一個位於$(i,j)$ 的元素時，$i$ 與 $j$ 保持不動，只有 $k$ 在保持變動，所以我們可以了解到：若是 $i$ 或是 $j$ 有任一個變動，就是已經算完該元素了，這個結論對於以下了解 Horowitz et al. 所寫的程式中的判斷式尤為有用：
我們用稀疏矩陣記法，先將 $i,j,k$ 三個變數填上去：
| Matrix | Row | Col | Val |
|---|---|---|---|
a |
$i$ | $k$ | 該元素值 |
b |
$k$ | $j$ | 該元素值 |
b的 transpose: newB |
$j$ | $k$ | 該元素值 |
所以，我們可以大致建構以下的程式進行矩陣間的元素乘法：
void mmult(term a[], term b[], term d[]) {
term newB[MAX_TERMS];
fastTranspose(b, newB);
int column, row = a[1].row, sum = 0;
int totalD = 0;
int rowsA = a[0].row, colsB = b[0].col;
int rowBegin = 1;
a[a[0].val+1].row = rowsA;
newB[b[0].val+1].col = colsB;
newB[b[0].val+1].col = 0;
for (int i = 1; i <= a[0].val; ) {
column = newB[1].row;
for (int j = 1; j <= b[0].val+1; ) {
if (a[i].row != row) {
storeSum(d, &totalD, row, column, &sum);
i = rowBegin;
for(; newB[j].row == column; j++);
column = newB[j].row;
} else if (newB[j].row != column) {
storeSum(d, &totalD, row, column, &sum);
i = rowBegin;
column = newB[j].row;
} else switch (COMPARE(a[i].col, newB[j].col)) {
case -1:
i++;
break;
case 0:
sum += (a[i++].val * newB[j++].val);
break;
case 1:
j++;
}
}
for (; a[i].row == row; i++);
rowBegin = i;
row = a[i].row;
}
d[0].row = rowsA;
d[0].col = colsB;
d[0].val = totalD;
}
另外有一支 Function storeSum：
void storeSum (term d[], int *totalD, int row, int column, int *sum) {
if (*sum) {
if (*totalD < MAX_TERMS) {
d[++*totalD].row = row;
d[*totalD].col = column;
d[*totalD].val = *sum;
*sum = 0;
}
}
}
mmult兩層迴圈中的判斷式主要判斷：
- 是否仍在同 row 運算中，若是不同 row 則結束該元素運算，並儲存該元素總和。（
if (a[i].row != row)的意思） - 是否仍在同 col 運算中，若是不同 col 則結束該元素運算，並儲存該元素總和。（
else if (newB[j].row != column)的意思，這裡的newB已是由B轉置過來的稀疏矩陣） - 以上條件皆否，表示仍在同一個 row 與 col 之中，要確保其元素是否一一對準（上面公式 $k$ 的作用，要確保 $k$ 值要保持一致），若是 a 的 k 值（
a[i].col）小於 newb 的 k 值（newB[j].col）則讓i++讓 $k$ 值可追上 newb，若是 a 的 k 值（a[i].col）大於 newb 的 k 值（newB[j].col）則讓j++讓 $k$ 值可追上 a。
要值得注意的是，有可能你會懷疑增加
i, j值可確保 $k$ 值增加嗎？由於我們先前已經確定我們在稀疏矩陣記法時已經保證在同一個 row index 底下，col index 會保持遞增，而我們這邊要已經在同一個 row 與 col（newB 的 row，B 的 col，newB 是 B 的轉置矩陣） 的前提下，a.row與newB.row都已經固定了，所以確保他們兩個 col 會保持遞增情形。
所以，何時才會執行到switchcase 的case 0呢？也就是在a[i].col與newB[j].col值相同時，才會進行兩矩陣的元素相乘與加總；反之，就i++, j++移至下一個 data row，所以也就是for迴圈不加上執行後計數的指令了。
第二層迴圈的終止條件與第一層的終止迴圈差一個「+1」，因為下方判斷式緣故，必須「col 與目前相異才會記錄該元素值」，也就是如果我只停在「b[0].val」，那我雖然sum已經記錄完該元素的所有東西，但是沒辦法進入到newB[j].row != column內進行記錄。
Share this post
Twitter
Google+
Facebook
Reddit
LinkedIn
Pinterest
Email