--- title: 數學軟體實作 - 稀疏矩陣 tags: 2020 Fall - 數學軟體實作 GA: G-77TT93X4N1 --- # Lecture - Solving sparse linear system ### Example 1 如果某矩陣大部分的元素為零, 只有少部分的非零元素, 那麼一個比較有效率的做法是不要將整個矩陣的元素都存下來, 只需要存非零部分就好. 這種矩陣稱為稀疏矩陣(sparse matrix). 舉例來說, 以下這矩陣只有對角線部分有值, 其他皆為零, $$ \begin{bmatrix} 2 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 2 \end{bmatrix} $$ ```matlab= A = diag([2 1 2]) % 產生一個只有對角線有值的矩陣 S = sparse(A) % 將 A 變成稀疏矩陣 ``` 執行結果為: ``` A = 2 0 0 0 1 0 0 0 2 S = (1,1) 2 (2,2) 1 (3,3) 2 ``` 可以看出兩種作法的差別, $A$ 記全部元素, $B$ 則只記非零的部分. 再來試更大的矩陣, 這次是 $10000\times 10000$: ```matlab= A = diag(rand(10^4, 1)); % 產生一個只有對角線有值的矩陣 S = sparse(A); % 將 A 變成稀疏矩陣 whos % 看看變數細節 ``` 可以看出 $A$ 跟 $S$ 都是 $1000\times 1000$ 的矩陣, 不過所需要存的記憶體則有約三千倍的差別 (Bytes 那欄). ``` Name Size Bytes Class Attributes A 10000x10000 800000000 double S 10000x10000 240008 double sparse ``` 所以, 如果一個矩陣很稀疏, 那就務必存成稀疏矩陣的格式. --- ### Example 2 再舉個例子, 假設我們想解 $$ D x = b, $$ 其中 $D$ 是個只有對角線有值的矩陣, $D_{k,k} = 1+\frac{1}{k}$, $b$ 則是個向量, $b_k=1$. $D$ 是個 $n\times n$ 矩陣, 我們有三種作法來解這問題. 先分別將三種作法所需要的矩陣造出來: ```matlab= n = 10^4; % 矩陣 size 為 nxn d0 = 1+1./(1:n)'; % 對角線自己存成一向量 D1 = diag(d0); % 造出整個矩陣 D2 = sparse(D1); % 造出稀疏矩陣 b = ones(n,1); % 等號右邊 ``` 以下我們要分別計時一下, 看三種做法速度差多少 #### Method 1 因為我們知道這問題的解是 $$ x_k = \frac{b_k}{d_{k,k}}, $$ 因此我們可以這樣做 (也就是直接把 $b$ 向量點除矩陣對角線向量): ```matlab= tic x = b./d0; toc ``` 第一次 run 程式時間會比較長, 建議多跑幾次取時間的平均. 我所得的時間為 ``` Elapsed time is 0.015570 seconds. ``` #### Method 2 我們知道解線性系統 $Dx=b$ 問題在 `matlab` 裡的指令為 `x = b\D`, 所以我們就這樣做: ```matlab= tic x = D1\b; toc ``` 結果是 ``` Elapsed time is 0.646748 seconds. ``` #### Method 3 第三種一樣是反除, 只不過我們這次是反除那個稀疏矩陣 ```matlab= tic x = D2\b; toc ``` 結果是 ``` Elapsed time is 0.015434 seconds. ``` #### Conclusion 可以看出, 直接做跟對稀疏矩陣反除效果都很好, 不過對整個大矩陣反除就慢多了. --- #### Rule: 以上這種將矩陣造出, 再 `sparse` 變成稀疏矩陣的作法純屬課堂效果, 絕對不可以用. 不過我們還是介紹兩個指令: `sparse`: 如同前面的例子, 將一個矩陣轉成稀疏矩陣 `full`: 這是反過來, 將一個稀疏矩陣變成一般矩陣 ```matlab= A = diag([2 1 2]) % 產生一個只有對角線有值的矩陣 S = sparse(A) % 將 A 變成稀疏矩陣 B = full(S) % 將 Ａ 變回來 ``` --- ## 稀疏矩陣指令 生成稀疏矩陣有幾個重要指令: #### `speye`: sparse identity matrix `speye(m, n)`: 生成一個稀疏的 $m\times n$ 矩陣, 其中主對角線都是 $1$, 其餘所有元素都是 $0$. #### `spdiags`: sparse matrix formed from diagonals `spdiags(B,d,m,n)`: 生成一個 $m\times n$ 的稀疏矩陣, $B$ 是個 $m\times 1$ 向量, $d$ 則是要將這向量放在哪個位置, $d=0$ 表示對角線, $d=1$ 表示對角線往上移一排, $d=-1$ 則表示對角線往下移一排. 以下兩個範例: ```matlab= A1 = spdiags([1;2;3], 0, 3,3) A2 = full(A1) ``` ``` A1 = (1,1) 1 (2,2) 2 (3,3) 3 A2 = 1 0 0 0 2 0 0 0 3 ``` ```matlab= B = [2 7; 4 9; 1 3]; d = [1 -1]; S = spdiags(B, d, 4, 3) A = full(S) ``` ``` S = (2,1) 7 (1,2) 4 (3,2) 9 (2,3) 1 (4,3) 3 A = 0 4 0 7 0 1 0 9 0 0 0 3 ``` ### Exercise 1 改寫 Jacobi method 的程式, 確保他執行速度夠快. ### Exercise 2 使用 Jacobi method 來解以下問題, $n=10^5$, $$ \begin{bmatrix} -2 & 1 & & & & \\ 1 & -2 & 1 & & & \\ & 1 & \ddots & \ddots & & \\ & & \ddots & \ddots & 1 & \\ & & & 1 & -2 \end{bmatrix}_{n\times n} x = b, $$ 其中 $$ b_k = \frac{1}{(n+1)^2}\sin\left(\frac{k\pi}{n+1}\right). $$