第四章、 數值方法
4.1 Lanczos Method
首先，考慮標準特徵値問題
μ
=
Ax x (4.1.1)
A 為一個 N x N 的矩陣，μ為特徵值，X 為特徵向量，如果我們從任意向
量v1開始，在這 Lanczos 演算法中，我們會利用連續的向量，當做新的基底
1 1 1
2
1 1 1 , 1 1 1
p+ = +p+ p+ ap+ = +p+ p+
b w w u Hu (4.1.15)
我們使用此方法決定其他三條對角線(tridiagonal)矩陣 T 的其他元素。
所以
此為一個相似變換(similar transformation)，其中一種特性為 A 與 T 有相同 特徵値，便可以利用 Lanczos 方法去找出ui使 A 能夠轉成 T，讓對角化簡單
1 1
但事實上，正因為三對角線(Tridiagonal)矩陣 T 的關係，使正交性快速的減 低，T 的特徵值是相近於 A 的特徵值，並且在某些情況之下，有可能產生 的新基底是與舊基底線性相關(linear dependent)的，所以 T 的特徵值也有可 能並不接近 A 的特徵值，因此，Lanczos 演算法仍然不是很穩定。
Lanczos 所需改善的地方：
1.避免正交性的損失。
2.使現有基底的正交。
3.移除假的特徵值。
因此，基本的 Lanczos 的演算法中，做為大型矩陣計算方面，需要再做許多 修改，接下來介紹的 ARPACK，克服了某些 Lanczos 的問題。
4.2 ARPACK 簡介
ARPACK 是一個用來求解大規模特徵值問題的軟體。ARPACK 是由
Fortran77 所寫的副程式的集合，主要可以解決來自各應用領域的大規模對 稱、非對稱與廣義特徵值問題。可以求特定數量的特徵值，其中包含最大 或最小實部特徵值、最大或最小虛部特徵值、以及最大或最小振幅的特徵 值。ARPACK 軟體是基於隱式重開始 Arnoldi/Lanczos 方法(Implicitly Restarted Arnoldi/Lanczos Method)。Arnoldi/Lanczos 是用來求解大規模 n*n 矩陣的幾個特徵值的重要方法，特別適合用於大規模稀疏矩陣的特徵值問 多的矩陣，例如解決固態物理中所使用到的 Tight-Binding Method、解決微 分方程的有限差分法(Finite Difference Method)，都是矩陣中具有相當多為零 的元素，使用 ARPACK 才會有較佳的結果。
主要功能包括：
1. 求解單雙精度實對稱矩陣特徵值問題
2. 求解雙精度複數矩陣特徵值問題(Hermitian & non-Hermitian) 3. 廣義特徵值問題
4. 奇異值分解(SVD)問題
最新軟體版本可上 http://www.caam.rice.edu/software/ARPACK/ 自行下載使 用，網站上也有線上的 User Guide 可以供大家使用。
Implicit Restart Arnoldi Method (IRAM)
圖 4.2.1、IRAM 流程圖
Arnoldi Factorization
在介紹 IRAM 之前，首先要瞭解 k 階 Arnoldi Factorization。先考慮一個 一般特徵值問題，定義：如果待求矩陣A∈Cn n× ，且滿足
T k = k k+ k k
AU U H f e (4.2.1)
其中Uk∈Cn k× 為正交歸一的行向量，稱為 Arnoldi 向量，且U fHk k =0，
∈ k k×
k C
H 是上海森堡矩陣(upper Heisenberg matrix)，為下對角線以上皆有值 的矩陣，eTk =
(
L 0 1 0 L)
在第 k 個位置才有 1 其它為零，此關係為 A 的 k 階 Arnoldi Factorization，為 Krylov 子空間的一種應用。若 A 是 Hermitian，則Hk是實數對稱三對角線矩陣(real symmetric tridiagonal matrix)，此關係稱
為 A 的 k 階 Lanczos 分解，所對應的行向量Uk稱為 Lanczos 向量。理論上
Filter polynomials
假設在n n× 矩陣 A 中有特徵值μi與特徵向量xi，而有興趣的是前 k 個特
與 A 有關，則 f 為 A 的特徵多項式 積，直接處理將此 Filter 多項式乘上向量，稱為顯式重開始(Explicit
Restarting)。若矩陣很大時，顯式重開始勢必會花費許多的時間，因此，隱 式重開始便因應而生。
Implicit Restart Arnoldi Method
針對顯式重開始的不方便，因此發展了隱式重開始，先將 Arnoldi 分解
由此可以看出，藉由隱式位移 QR 分解，可以得到與顯式重開始相同效 果的基底，詳細推導內容請參見附錄。由於e QTm 的前 k-1 項為零，所以此 m 階的前 k 項可以等價改寫成 k 階 Arnoldi 分解，或等同只取前 k 項，本質上 是一樣的：
T
k k k k k
+ = + ++ +
AU U H f e (4.2.7)
其中，殘餘向量fk+ =U ek+ k+1βˆk+fmσ ，此時的 k 階 Arnoldi 分解，可以包含 此 m 階 Arnoldi 分解中，有興趣的特徵值所展延的特徵空間，並且起始向量 也將不想要的特徵值過濾掉，進而再增加 p 階使 Arnoldi process 產生下一個 新的 m 階 Arnoldi 分解，再進行隱式位移 QR 法，直到收斂為止。利用所得 到的該 Krylov 子空間的基底向量與顯式重開始v1所產生的空間相同的特 性，便可以大大減少產生 Krylov 子空間所需要的計算量。此為 ARPACK 重 要的精神。
根據上述方法，可編成下列演算法：
( )
an m step Arnoldi Factorization For l until convergence
Computer and select set of p shifts based upon or perhaps other information
For j p
Beginning with the k step Arnoldi factorization
apply p additional step of the Arnoldi process to obtain a new m step
特徵空間、Krylov 子空間與稀疏矩陣紀錄方法
眾所皆知的，一個向量可以展延(span)成一維空間，二個獨立向量可以 展延(span)成二維空間，以此類推，當有 N 個獨立向量時，便可以展延成無 限維空間，在物理上稱為 Hilbert 空間(Hilbert space)。
特徵空間是用來描述一個由 A 的特徵值相對應特徵向量所形成的空
4.2.2 Test Result
我們是使用有限差分法(finite difference method)[15]將微分方程轉成矩陣，
再利用 ARPACK 去執行對角化所產生的結果，其正解應為 0.5、1.5、2.5...
等，以下是測試的結果
圖 4.2.1、矩陣 10 萬*10 萬的特徵值
圖 4.2.2、ARPACK 的矩陣大小與時間分析圖
圖 4.2.3、LAPACK 的矩陣大小與記憶體分析圖
圖 4.2.4、ARPACK 的矩陣大小與記憶體分析圖
圖 4.2.5、LAPACK 的矩陣大小與記憶體分析圖
可以看出以 ARPACK 與 LAPACK 相比之下，記憶體大小大幅度的降低，計 算大小大幅度的增加。
ARPACK 的所有參數皆列在附錄 A 中，其中最主要需要更改的參數如 下所示：
參數 功能描述
n 矩陣大小
nev 求解特徵值的個數
ncv Arnoldi 分解長度，即基底數量；預設值=4*nev which 指定求解特徵值的特性(參考附錄 A) iparam(3) 允許疊代的最大次數；預設值=3*n
其中，主要控制運算時間長短的兩個參數為 ncv 與 iparam(3)，增加基底
數 ncv 能降低運算時間，其主要是將待解矩陣分割成小矩陣的大小，ncv 越 大，所產生小矩陣的大小越大，疊代次數會較少，但是每次疊代後解小矩 陣的時間越多；反之，ncv 越小，所產生的小矩陣的大小越小，解小矩陣的 時間越少，但疊代次數會較多。舉例來說：若待解矩陣為 1000*1000 的大 小，ncv 假設是 200，則每次的小矩陣為 200*200，可能的疊代次數假設是 5 次，即解 5 次 200*200 大小的矩陣；若 ncv 假設是 10，則每次的小矩陣 為 10*10，可能的疊代次數假設是 100 次，即解 100 次 10*10 的矩陣。
疊代次數(iparam(3))只需要修改到滿足程式未判定疊代次數不夠即可。
圖 4.2.6、ARPACK 參數修改對時間的影響，前項參數為 ncv(basis)=4*nev，
後項參數為 iparam(3) =3*N
在經過一連串的測試之後，針對求一維簡諧振子最低的 15 個能量，測 試出基底數量取 12*nev (nev=15)是可以節省最多時間的，在此提供給大家 參考。但是一般來說，不同的矩陣對於參數的優化並沒有一個固定的公式，
端看不同矩陣的稀疏程度、特徵值的結構與個數的而決定，因此，針對任 何一個新的問題，便必須經由反覆的測試才能得到最佳化的參數。
第五章、總結與未來工作
之後，可以慢慢增加計算複數矩陣，或是加入能平行處理矩陣的 PARPACK，都是未來的計算課題。
參考文獻
[1]Novoselov, K.S. et al. "Electric Field Effect in Atomically Thin Carbon Films", Science, Vol 306 (5696), p. 666-669, 2004
[2] Novoselov, K.S. et al. "Two-dimensional atomic crystals", PNAS, Vol 102 (30), p.
10451-10453, January 26, 2005
[3] Mark Wilson, Phys. Today 59 January 21 (2006)
[4] A.K. Geim & K.S. Novoselov, Nature Materials 6, 183-191 (2007)
[5] J. R. Williams, L. C. DiCarlo, and C. M. Marcus, Science 317, 638 (2007) [6] D. A. Abanin and L. S. Levitov, Science 317, 641 (2007)。
[7]Barbaros Özyilmaz et al.，Phys. Rev. Lett. 99, 166804 (2007)
[8]Marek Korkusinski,”Application of the effective band-orbital model to semiconductor nanostructure”
[9]Gene Howard Golub, Charles F. Van Loan,”Matrix computations”(1996)
[10]R. B. Lehoucq, D. C. Sorensen, C. Yang,” ARPACK Users' Guide:Solution of Large Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods”(1997)
[11]A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov and A. K. Geim, Rev.
Mod. Phys. 81, 109 (2009)
[12] K. Nakada, M. Fujita, G. Dresselhaus, M. S. Dresselhaus, Phys. Rev. B 54, 17954 - 17961 (1996)
[13] Peter Y. Yu,Manuel Cardona,”Fundamentals of Semiconductors:Physics and Materials Properties”(springer-Verlag 1996)
[14] D. J. Chadi and M. L. Cohen Calculations of the Valence Bands phys. stat. sol. (b) 68, 405 (1975)
[15] L. V. Fausett, “Applied Numerical Analysis Using Matlab”, Prentice-Hall(1999) [16] R. Saito, M. Fujita, G. Dresselhaus, M.S. Dresselhaus , Phys. Rev. B 46, 1804 - 1811
(1992)
[17] R. Saito, G. Dresselhaus, M.S. Dresselhaus ,”Physical Properties of Carbon Nanotubes”,(Imperial College Press, Convent Garden, London 1998)
[18]A. Varykhalov, J. Sa´nchez-Barriga, A. M. Shikin, C. Biswas, E. Vescovo, A. Rybkin, D.
Marchenko, and O. Rader, Phys. Rev. Lett 101 157601(2008)
附錄 A：隱式位移 QR 法
one step of single shifted QR algorithm
proof easily
singular
non singular singular μ
( ) ( )
( ) ( )( )
附錄 B：ARPACK 程式參數
若設定了以上變數，則需相對修改 maxn、maxnev、maxncv 滿足下式：
maxn ≧ n maxncv ≧ ncv maxnev ≧ ncv ncv ≧nev+1 (在 C code 中並無 max 這些參數)
根據需要，下列變數可以依照其功能修改：
lwork1：陣列 work1 的大小，設定為工作的陣列大小，至少應被設置 為 ncv*(ncv+8)。
tol：收斂標準，此值用來控制收斂精確度。tol 越小，達到停止標準所需要 的計算量越大
ido：逆通訊標誌，進入**aupd 之前要設為 0
bmat：表示解決的是標準特徵值問題(bmat=’I’)或是廣益特徵值問題 (bmat=’G’)
iparam(1)：是否在隱式 QR 法中使用位移。若不清楚可設定為(iparam(1)=1) iparam(3)：允許的 IRAM 疊代的最大次數
iparam(7)：ARPACK 的運算模式