在一些物理問題的應用裡，矩陣的某些特徵值會是我們特別想知道的，在這裡我們的目標是求一個廣義特徵值問題裡任意指定的第k個特徵值（以由小到大順序排列）。

hspace{10mm} 根據線性代數裡的Sylvester''s law of inertia，我們可以藉由$LDL^T$分解來知道一個矩陣在任一個給定的值之前會有多少個特徵值。如果把給定的值之前有多少個特徵值當成一個函數f來看（這會是一個階梯函數），那麼尋找矩陣的第k個特徵值可以轉變為一個求根的問題( 解f(x)=k )。對於這個問題而言，我們的目標是以某些求根的演算法來找到一個實數值s滿足f(s)=k或k-1，這樣我們就能以某些shift-invert eigenvalue solver來找到精確的第k個特徵值。在這篇論文裡我們不著重在特定eigenvalue solver的選擇，而是將重點放在發展求根演算法以及引入一個可以減少函數值計算時間的方法(相對於$LDL^T$分解)。

hspace{10mm} 在先前的文獻裡，目前可用於階梯函數的穩定求根演算法已有二分法（Bisection method），而在這裡，首先我們將其推廣至一個適合於平行的多分法（Multiple-section method）。接著，我們也發展了一個基於單調三次分段插值的求根演算法，並將其和可平行的多分法作結合。在同樣的條件下（所採用平行節點的個數相同），相較於多分法，它能在大部份的情況下達成降低函數值總共計算次數的目的。

hspace{10mm} 另一方面，除了以$LDL^T$作為精確函數值的計算，我們還嘗試了一個用來估計特徵值個數的近似方法作為$LDL^T$分解的替代方案來達到減少單次函數值計算時間的目的。數值實驗的結果告訴我們通常在大型稀疏矩陣的情況下，我們可以藉由犧牲一些函數值的精確性來換取函數計算時間上的優勢。

In some applications of physics problems, we may want to know certain eigenvalues we''re interested in. Here our target is an arbitrarily chosen k-th eigenvalue of a generalized eigenvalue problems.

According to the Sylvester''s law of inertia in Linear algebra, we may get the eigenvalue counts before an arbitrarily given number via $LDL^T$ decomposition. If we view the eigenvalue counts before a given number as a step function f, then our problem of finding the k-th eigenvalue can be transformed into a root-finding problem ( solving $f(x)=k$ for x ). As far as our problem concerned, our goals is to find a real value $s$ such that $f(s)=k$ or $f(s)=k-1$ so that we can find the exact k-th eigenvalue by an shift-invert eigenvalue solver. In this paper, we do not emphasize the choice of a specific eigenvalue solver, but focus on developing the root-finding algorithm and introduce an approximating method which can reduce time consuming relative to $LDL^T$ decomposition.

In the literature, the Bisection method is an existing stable root-finding algorithm which can be applied to step functions. Here we will first extend the Bisection method to the Multiple-section method which can be naturally parallelized. We will also developed a sequential method based on piecewise Monotone cubic interpolation and then combine it with Multiple-section method to get a parallel algorithm. With the same number of MPI processes, it can reduce the total number of function evaluations for one time root-solving in most cases relative to the pure Multiple-section method.

On the other hand, besides computing the exact function value via $LDL^T$ decomposition, we also
took a stochastic scheme for estimating the eigenvalue counts as an alternative to achieve the aim for reducing the time consuming for one-time function evaluation. The result of numerical experiment shows that in the case of large sparse matrix, we can save function evaluation time by sacrificing some exactness of function value.

Table of Contents i
List of Figures iii
List of Tables iv
中文摘要v
Abstract vii
1 Introduction 1
2 Background 2
2.1 Sylvester’s law of inertia and its application to our problem . 2
2.2 Saad’s scheme for approximating the eigenvalue counts . . . 6
2.2.1 Polynomial expansion filtering . . . . . . . . . . . . . 6
2.2.2 Estimation of the trace with a stochastic estimator . 8
2.2.3 Estimation of the eigenvalue counts in generalized
eigenvalue problem . . . . . . . . . . . . . . . . . . . 9
2.2.4 How to choose a proper eigenvalue range in generalized
eigenvalue problem . . . . . . . . . . . . . . . . . 10
3 Methods 12
3.1 Root-Finding Algorithm . . . . . . . . . . . . . . . . . . . . 12
3.1.1 Bisection Method . . . . . . . . . . . . . . . . . . . . 13
3.1.2 Multiple-section Method . . . . . . . . . . . . . . . . 14
3.1.3 Monotone Cubic Interpolation Based Methods . . . . 15
3.1.4 A Hybrid method of Multiple-section and Monotone
Cubic Interpolation . . . . . . . . . . . . . . . . . . . 23
3.2 Combining Root-Finding Algorithm with Saad’s Scheme . . 25
4 Numerical Results 26
4.1 Number of LDLT Required between the Sequential Methods 27
4.2 Number of LDLT Required between Different Section Number
for Multiple-section Method . . . . . . . . . . . . . . . . 31
4.3 Number of LDLT Required between the Parallel Methods . . 35
4.4 Time comparison for LDLT Decomposition and Saad’s scheme
36
4.5 Error of Saad’s scheme . . . . . . . . . . . . . . . . . . . . . 37
4.6 Convergence History and Total computation time for Hybrid
Method combined with Saad’s scheme . . . . . . . . . . . . . 38
Bibliography 39

[1] J Cerda and F Soria. Accurate and transferable extended hu&;#776;ckel-type
tight-binding parameters. Physical Review B, 61(12):7965, 2000.
[2] Edoardo Di Napoli, Eric Polizzi, and Yousef Saad. Efficient estimation of
eigenvalue counts in an interval. arXiv preprint arXiv:1308.4275, 2013.
[3] Frederick N Fritsch and Ralph E Carlson. Monotone piecewise cubic
interpolation. SIAM Journal on Numerical Analysis, 17(2):238–246, 1980.
[4] Takeo Hoshi, Susumu Yamamoto, Takeo Fujiwara, Tomohiro Sogabe, and
Shao-Liang Zhang. An order-n electronic structure theory with generalized
eigenvalue equations and its application to a ten-million-atom system.
Journal of Physics: Condensed Matter, 24(16):165502, 2012.
[5] MF Hutchinson. A stochastic estimator of the trace of the influence
matrix for laplacian smoothing splines. Communications in Statistics-
Simulation and Computation, 18(3):1059–1076, 1989.
[6] Dongjin Lee, Takafumi Miyata, Tomohiro Sogabe, Takeo Hoshi, and
Shao-Liang Zhang. An interior eigenvalue problem from electronic structure
calculations. Japan Journal of Industrial and Applied Mathematics,
30(3):625–633, 2013.