Lecture - Solving sparse linear system

Example 1

如果某矩陣大部分的元素為零, 只有少部分的非零元素, 那麼一個比較有效率的做法是不要將整個矩陣的元素都存下來, 只需要存非零部分就好. 這種矩陣稱為稀疏矩陣(sparse matrix).

舉例來說, 以下這矩陣只有對角線部分有值, 其他皆為零,

$$\left[\begin{array}{ccc}2& 0& 0\\ 0& 1& 0\\ 0& 0& 2\end{array}\right]$$

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$:

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

再舉個例子, 假設我們想解

$$Dx=b,$$
其中

$D$ 是個只有對角線有值的矩陣,

$D_{k,k}=1+\frac{1}{k}$,

$b$ 則是個向量,

$b_k=1$.

$D$ 是個

$n\times n$ 矩陣, 我們有三種作法來解這問題. 先分別將三種作法所需要的矩陣造出來:

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$ 向量點除矩陣對角線向量):

tic
x = b./d0;
toc

第一次 run 程式時間會比較長, 建議多跑幾次取時間的平均. 我所得的時間為

Elapsed time is 0.015570 seconds.

Method 2

我們知道解線性系統

$Dx=b$ 問題在 matlab 裡的指令為 x = b\D, 所以我們就這樣做:

tic
x = D1\b;
toc

結果是

Elapsed time is 0.646748 seconds.

Method 3

第三種一樣是反除, 只不過我們這次是反除那個稀疏矩陣

tic
x = D2\b;
toc

結果是

Elapsed time is 0.015434 seconds.

Conclusion

可以看出, 直接做跟對稀疏矩陣反除效果都很好, 不過對整個大矩陣反除就慢多了.

---

Rule:

以上這種將矩陣造出, 再 sparse 變成稀疏矩陣的作法純屬課堂效果, 絕對不可以用.

不過我們還是介紹兩個指令:
sparse: 如同前面的例子, 將一個矩陣轉成稀疏矩陣
full: 這是反過來, 將一個稀疏矩陣變成一般矩陣

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$ 則表示對角線往下移一排.

以下兩個範例:

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

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$,

$${\left[\begin{array}{cccccc}-2& 1& & & & \\ 1& -2& 1& & & \\ & 1& \ddots & \ddots & & \\ & & \ddots & \ddots & 1& \\ & & & 1& -2\end{array}\right]}_{n\times n}x=b,$$
其中

$$b_k=\frac{1}{(n+1{)}^2}\sin \left(\frac{k\pi }{n+1}\right).$$