行列の総合レポート

0. 読み方・記法

行列を「数表」としてではなく、線形写像・基底変換・二次形式・データ変換・計算手続きの統一言語として扱います。Jordan 標準形は既知として、数値線形代数で実際に安定に使われる分解へ重点を移します。

1. 行列とは何か：線形写像・座標・合成

1.1 線形写像の座標表示

有限次元ベクトル空間 \(V,W\) と基底 \(\mathcal B=(v_1,\dots,v_n)\)、\(\mathcal C=(w_1,\dots,w_m)\) を固定すると、線形写像 \(T:V\to W\) は一意な行列 \([T]_{\mathcal C\leftarrow\mathcal B}\in\mathbb F^{m\times n}\) で表されます。

この意味で、行列の列は「基底ベクトルの像」を並べたものです。第 \(j\) 列は \(T(v_j)\) の \(\mathcal C\) 座標です。

1.2 行列積は写像の合成

\(T:V\to W\)、\(S:W\to U\) なら、

行列積の非可換性は、写像の合成順序が一般に可換でないことの座標表示です。

1.3 基底変換と相似

同じ線形作用素 \(T:V\to V\) を基底 \(\mathcal B\) と \(\mathcal B'\) で表すと、基底変換行列 \(P=[\operatorname{id}]_{\mathcal B\leftarrow\mathcal B'}\) により

この関係を相似といいます。固有値、特性多項式、最小多項式、Jordan ブロック構造、trace、determinant は相似不変です。一方、個々の成分、行和、列和、対角成分そのものは一般には不変ではありません。

1.4 双線形形式・二次形式としての行列

行列 \(A\) は \(x^*Ay\) により双線形または半双線形形式を定めます。Hermitian 行列 \(A=A^*\) の二次形式 \(x^*Ax\) は実数値です。Rayleigh 商

は固有値、最適化、正定値性、反復法の解析で中心的役割を持ちます。

2. 基本構造：rank、核、像、ブロック、射影

2.1 四つの基本部分空間

\(A\in\mathbb F^{m\times n}\) の rank を \(r\) とすると、四つの基本部分空間があります。

rank-nullity は

SVD はこの四部分空間を最も透明に表します。非零特異値に対応する右特異ベクトルは行空間を張り、零特異値に対応する右特異ベクトルは零空間を張ります。

2.2 ブロック行列と Schur 補

ブロック行列

で \(B\) が正則なら、\(B\) に関する Schur 補を

と定義します。このとき

2.3 射影行列

\(P^2=P\) を射影といいます。さらに \(P=P^*\) なら直交射影です。列フルランク \(A\in\mathbb F^{m\times n}\) に対し、\(\mathcal R(A)\) への直交射影は

ただし数値計算では \(A^*A\) を直接作るより、QR 分解 \(A=QR\) から \(P_A=QQ^*\) と扱う方が安定です。

3. ノルム・条件数・浮動小数点

3.1 ベクトルノルムと誘導ノルム

代表的なベクトルノルムは

行列ノルムがベクトルノルムから誘導されるとは

であることです。代表例は

Frobenius ノルムは誘導ノルムではありませんが、特異値と trace で

と書けます。

3.2 条件数

正則行列 \(A\) の線形方程式 \(Ax=b\) に対する条件数は

2-ノルムでは SVD により

摂動 \(A+\Delta A\), \(b+\Delta b\) に対し、\(\|A^{-1}\Delta A\|<1\) なら概略として

3.3 浮動小数点と後退安定性

丸めモデルでは、基本演算 \(\circ\in\{+,-,\times,/\}\) について

と書きます。\(u\) は unit roundoff です。

後退安定とは、計算結果 \(\widehat x\) が少し摂動した問題の厳密解になっていることです。線形方程式なら

Householder QR や Cholesky は典型的に安定です。部分ピボット付き LU も実用上は非常に安定ですが、最悪ケースでは成長因子が大きくなり得ます。

4. 行列クラスの分類

行列クラスを認識できると、分解・固有値・計算量・安定性が一気に見通せます。

4.1 対角・スカラー・置換行列

対角行列 \(D=\operatorname{diag}(d_1,\dots,d_n)\) は成分ごとに作用します。\(D^{-1}\) は全 \(d_i\ne0\) のとき存在し、\(D^{-1}=\operatorname{diag}(d_i^{-1})\) です。置換行列 \(P\) は標準基底を並べ替える行列で、\(P^{-1}=P^T\) です。

4.2 三角行列

上三角行列 \(U\) は \(i>j\Rightarrow u_{ij}=0\) を満たします。積、逆行列、固有値に関して閉じています。固有値は対角成分です。

三角行列の連立方程式は \(O(n^2)\) で解けるため、LU、QR、Schur の最終形として頻出します。

4.3 Hermitian・正規・ユニタリ・正定値

4.4 帯行列・三重対角行列・Hessenberg 行列

帯行列では非零成分が対角近傍に限定されます。半帯幅 \(p,q\) を用いると \(a_{ij}=0\) for \(i-j>p\) or \(j-i>q\)。三重対角行列は \(p=q=1\) です。上 Hessenberg 行列は \(i>j+1\Rightarrow h_{ij}=0\)。一般行列の QR 固有値アルゴリズムでは、まず Hessenberg 形へ縮約します。Hermitian 行列なら三重対角形へ縮約できます。

4.5 Toeplitz・巡回行列

Toeplitz 行列は \(a_{ij}=t_{i-j}\)、つまり各対角線上で成分が一定です。巡回行列はさらに各行が循環シフトであり、離散 Fourier 行列で対角化されます。

この構造は畳み込み、信号処理、FFT に直結します。

5. 分解の地図：何を、なぜ、どの順番で使うか

行列分解は、複雑な線形変換を「解きやすい」「直交的」「対角的」「低ランク的」な因子に分ける技術です。

6. LU 分解：数理・計算・安定性

6.1 Gaussian 消去と LU

正方行列 \(A\in\mathbb F^{n\times n}\) に対し、下三角行列 \(L\) と上三角行列 \(U\) で

と分解するのが LU 分解です。通常は \(L\) の対角を 1 とする Doolittle 型を使います。Gaussian 消去の乗数

を \(L\) の下三角成分に保存し、消去後の上三角部分が \(U\) になります。

6.2 Doolittle 公式

\(L\) は単位下三角、\(U\) は上三角とします。成分比較から

したがって、\(i\le j\) では

\(i>j\) では

for k = 1..n:
    for j = k..n:
        U[k,j] = A[k,j] - sum_{s=1}^{k-1} L[k,s] U[s,j]
    L[k,k] = 1
    for i = k+1..n:
        L[i,k] = (A[i,k] - sum_{s=1}^{k-1} L[i,s] U[s,k]) / U[k,k]

6.3 PLU：部分ピボット

ピボット \(a^{(k)}_{kk}\) が零または小さいと、消去が不可能または不安定になります。そこで列 \(k\) の下側から絶対値最大の成分を探して行交換します。

この結果は

と表されます。\(P\) は置換行列です。

Input: A in F^{n x n}
P = I
for k = 1..n-1:
    p = argmax_{i=k..n} |A[i,k]|
    if A[p,k] == 0: stop; matrix is singular to working precision
    swap rows k and p of A
    swap rows k and p of P
    swap L multipliers already stored in columns 1..k-1
    for i = k+1..n:
        A[i,k] = A[i,k] / A[k,k]              # multiplier l_ik
        for j = k+1..n:
            A[i,j] = A[i,j] - A[i,k] A[k,j]
Output: P, L = unit lower part of A, U = upper part of A

分解後、\(Ax=b\) は \(PAx=Pb\)、すなわち \(LUx=Pb\) として、前進代入 \(Ly=Pb\) と後退代入 \(Ux=y\) で解きます。

6.4 具体例

6.5 完全ピボット・rook pivoting・成長因子

部分ピボットは行交換のみです。完全ピボットは残り部分行列から最大絶対値成分を選び、行と列を交換します。

完全ピボットは理論的により安全ですが、探索コストと列交換管理が重く、通常は部分ピボットが標準です。対称不定値行列には Bunch–Kaufman 型の対称ピボットがよく使われます。

LU の誤差解析で重要なのが成長因子

後退誤差評価は概ね \(g\) に比例します。部分ピボットでは最悪 \(g\) が指数的に大きくなり得ますが、実用上は多くの問題で小さいままです。

6.6 ブロック LU

現代計算機では BLAS-3、つまり行列-行列積を活用するためブロック化が重要です。

すると

最後に Schur 補 \(S\) を再帰的に分解します。大きな更新 \(A_{22}\leftarrow A_{22}-L_{21}U_{12}\) が性能の鍵です。

6.7 行列式・逆行列・複数右辺

\(PA=LU\) なら

逆行列を明示的に作ることは多くの場合避けます。\(A^{-1}B\) が必要なら、各列について \(AX=B\) を LU の三角解法で解く方が良いです。右辺が複数ある場合、分解 \(PA=LU\) は一度だけ行い、三角解法を使い回します。

7. Cholesky 分解・LDL* 分解

7.1 正定値行列と Cholesky

Hermitian 正定値行列 \(A\succ0\) に対して、一意な下三角行列 \(L\) が存在し、対角成分は正で

と分解できます。LU の約半分の計算量で、ピボット不要かつ安定です。

7.2 成分公式

\(A=LL^*\) の成分比較から

for k = 1..n:
    s = A[k,k] - sum_{j=1}^{k-1} |L[k,j]|^2
    if s <= 0: stop; not positive definite in working precision
    L[k,k] = sqrt(s)
    for i = k+1..n:
        L[i,k] = (A[i,k] - sum_{j=1}^{k-1} L[i,j] conj(L[k,j])) / L[k,k]

7.3 LDL* 分解

平方根を避けたい場合、Hermitian 行列を

と分解します。\(L\) は単位下三角、\(D\) は対角またはブロック対角です。正定値なら \(D\) は正対角で済みます。不定値なら安定性のため 1×1 と 2×2 ピボットを使う Bunch–Kaufman 型分解が標準です。

7.4 正定値判定

Hermitian 行列 \(A\) について、次は同値です。

1. \(A\succ0\)。
2. 全固有値が正。
3. すべての先頭主座小行列式が正。
4. Cholesky 分解が失敗せず、対角成分が全て正。

7.5 例

8. QR 分解：直交性を使う分解

8.1 QR 分解の定義

\(A\in\mathbb F^{m\times n}\)、\(m\ge n\) に対し、列フルランクなら

と分解できます。ここで \(Q\in\mathbb F^{m\times n}\) は \(Q^*Q=I_n\) を満たす列直交行列、\(R\in\mathbb F^{n\times n}\) は上三角行列です。完全 QR では \(Q\in\mathbb F^{m\times m}\)、\(R\in\mathbb F^{m\times n}\)。

QR は最小二乗問題

を安定に解く基本手段です。

8.2 Gram–Schmidt

古典 Gram–Schmidt は列ベクトル \(a_j\) から直交基底を作ります。

丸め誤差のため、古典 Gram–Schmidt は列がほぼ線形従属だと直交性を失いやすいです。修正 Gram–Schmidt は投影を逐次的に更新し、より安定です。さらに必要なら再直交化を行います。

for j = 1..n:
    v = A[:,j]
    for i = 1..j-1:
        R[i,j] = Q[:,i]^* v
        v = v - Q[:,i] R[i,j]
    R[j,j] = ||v||_2
    Q[:,j] = v / R[j,j]

8.3 Householder QR

実務上の標準は Householder 反射です。ベクトル \(x\in\mathbb F^k\) に対し、\(Hx=\alpha e_1\) となるユニタリ反射

を作ります。実数なら \(\tau=2/(v^Tv)\)。桁落ちを避けるため、典型的には

を選びます。

R = A
for k = 1..n:
    x = R[k:m, k]
    construct Householder vector v so that (I - tau v v*) x = alpha e_1
    R[k:m, k:n] = R[k:m, k:n] - tau v (v* R[k:m, k:n])
    store v below diagonal of R
Q is represented implicitly as H_1 H_2 ... H_n

Householder QR は後退安定で、密行列では計算量 \(2mn^2-2n^3/3\) flops です。

8.4 Givens 回転

Givens 回転は平面 \((i,j)\) 内の回転で一つの成分を零にします。

疎行列、Hessenberg 行列、逐次更新 QR では Givens が便利です。

8.5 列ピボット付き QR と rank revealing

数値ランク判定には列ピボット付き QR を使います。

列ノルムが大きい列を先に選び、\(|r_{11}|\ge |r_{22}|\ge\cdots\) に近い構造を得ます。\(|r_{kk}|\) の急落は数値ランクの手がかりです。ただし最良低ランク近似の精度診断では SVD がより信頼できます。

8.6 QR と最小二乗

列フルランク \(A\in\mathbb F^{m\times n}\)、\(m\ge n\) とし、薄い QR \(A=QR\) を使います。すると

したがって解は三角方程式 \(Rx=Q^*b\) から得られます。正規方程式 \(A^*Ax=A^*b\) は \(\kappa_2(A^*A)=\kappa_2(A)^2\) と条件数を二乗するため、数値的には QR が標準です。

9. Schur 分解・固有値計算・Jordan 形との関係

9.1 Schur 分解

任意の正方複素行列 \(A\in\mathbb C^{n\times n}\) はユニタリ行列 \(Q\) と上三角行列 \(T\) により

と分解できます。\(T\) の対角成分が \(A\) の固有値です。実行列では実 Schur 分解として、\(T\) は 1×1 と 2×2 ブロックをもつ準上三角行列になります。

9.2 Hessenberg 縮約

QR 固有値アルゴリズムの前処理として、一般行列を上 Hessenberg 形へユニタリ相似変換します。

Hessenberg 形は QR 反復で保存され、1回の反復を \(O(n^2)\) にできます。Hermitian 行列では三重対角形に縮約できます。

9.3 QR アルゴリズム

基本 QR 反復は

各ステップはユニタリ相似なので固有値を保ちます。実装ではシフト \(\mu_k\) を使い、

Hessenberg 構造を保ちながら「bulge chasing」で実装します。

1. Householder 相似変換で A を Hessenberg H に縮約する。
2. H の右下 1x1 または 2x2 ブロックからシフトを選ぶ。
3. implicit QR step を行い、Hessenberg 構造を保つ。
4. 副対角成分が十分小さくなったら deflation する。
5. 残ったブロックに対して反復する。

9.4 Hermitian 固有値問題

\(A=A^*\) なら固有値は実数、固有ベクトルは直交基底に選べます。まず三重対角化

を行い、三重対角 \(T\) の固有値を QR、二分法、dqds、分割統治、MRRR などで求めます。三重対角化は \(O(n^3)\)、三重対角固有値計算は多くの場合それより低コストです。

9.5 固有値の摂動と pseudospectrum

正規行列では固有値摂動がよく制御されます。非正規行列では固有値が非常に敏感になり得ます。これを理解するには pseudospectrum

が有効です。固有値だけを見ると安定そうに見えても、非正規性が強いと \(e^{tA}\) や反復過程で大きな過渡増幅が起こります。

10. 特異値分解 SVD：理論・幾何・アルゴリズム

10.1 SVD 定理

任意の \(A\in\mathbb C^{m\times n}\) に対し、ユニタリ行列 \(U\in\mathbb C^{m\times m}\)、\(V\in\mathbb C^{n\times n}\)、非負対角成分をもつ \(\Sigma\in\mathbb R^{m\times n}\) が存在して

\(\Sigma\) の対角成分 \(\sigma_1\ge\sigma_2\ge\cdots\ge0\) を特異値といいます。rank \(r\) なら \(\sigma_1\ge\cdots\ge\sigma_r>0=\sigma_{r+1}=\cdots\)。薄い SVD は

10.2 幾何学的意味

SVD は線形写像を

に分解します。単位球は楕円体に写り、その半軸長が特異値です。

10.3 SVD と四つの基本部分空間

薄い SVD を \(A=U_r\Sigma_rV_r^*\) とします。さらに \(V=[V_r,V_0]\)、\(U=[U_r,U_0]\) と拡張すると

これが SVD の診断力の本質です。

10.4 低ランク近似：Eckart–Young–Mirsky

\(A=U\Sigma V^*\) とし、

すると rank \(k\) 以下の行列全体で、\(A_k\) は最良近似です。

画像圧縮、PCA、潜在意味解析、モデル縮約の中心定理です。

10.5 疑似逆

Moore–Penrose 疑似逆は

で定義されます。\(\Sigma^+\) は非零特異値を逆数にして転置した対角行列です。最小二乗解のうち最小ノルム解は

特異値が小さいと \(1/\sigma_i\) が大きくなりノイズを増幅します。このため切断 SVD や Tikhonov 正則化が使われます。

10.6 SVD の計算アルゴリズム概観

SVD を \(A^*A\) の固有値問題として解く方法は理論的に簡単ですが、条件数を二乗し、小さい特異値の精度を悪化させます。実用的な密行列 SVD は次の二段階が標準です。

1. Householder 変換で \(A\) を双対角行列 \(B\) へ縮約する。
2. 双対角 \(B\) の SVD を QR 型、dqds、分割統治などで求める。

10.7 Golub–Kahan 双対角化

\(m\ge n\) とします。左と右から Householder 変換をかけ、

となる上双対角行列

を作ります。

for k = 1..n:
    choose left Householder H_k to zero A[k+1:m, k]
    A[k:m, k:n] = H_k A[k:m, k:n]
    if k < n:
        choose right Householder G_k to zero A[k, k+2:n]
        A[k:m, k+1:n] = A[k:m, k+1:n] G_k
Result: bidiagonal B stored in A, with accumulated U and V if needed

密行列でのコストは概ね \(4mn^2-4n^3/3\) flops 程度です。特異ベクトルも必要なら左右の Householder を蓄積します。

10.8 双対角 SVD・Jacobi SVD・Lanczos・ランダム化 SVD

双対角 \(B\) の特異値は対称三重対角 \(T=B^*B\) の固有値の平方根です。ただし明示的に \(B^*B\) を作ると相対精度に問題が出るため、実装では双対角構造上で implicit QR や dqds を行います。

Jacobi 型 SVD は列対の直交化を反復します。非常に高精度ですが、大規模密行列では標準法より遅いことがあります。巨大疎行列や行列ベクトル積だけ可能な場合は Lanczos 双対角化を使います。

Input: A, target rank k, oversampling p, power iterations q
1. Draw random Ω in F^{n x (k+p)}
2. Y = A Ω
3. for t = 1..q: Y = A (A* Y), with orthonormalization
4. Orthonormalize Y = Q R
5. Form small B = Q* A
6. Compute SVD B = Ũ Σ V*
7. U = Q Ũ
Output: U[:,1:k], Σ[1:k], V[:,1:k]

10.9 具体例

11. 特殊行列クラス別の理論とアルゴリズム

11.1 Hermitian 行列

Rayleigh 商は

なので、\(\lambda_{\min}\le R_A(x)\le\lambda_{\max}\)。Courant–Fischer の min-max 原理は

Hermitian 行列のアルゴリズムでは、三重対角化、Lanczos 法、共役勾配法、Cholesky が基本です。

11.2 ユニタリ行列

\(U^*U=I\) は内積を保存します。

したがって \(\|Ux\|_2=\|x\|_2\)。数値計算でユニタリ変換が好まれる理由は、ノルムと条件を悪化させないためです。Householder、Givens、QR、Schur、SVD の中心にユニタリ行列があります。

ユニタリ行列の固有値は \(|\lambda|=1\) を満たします。特異値はすべて 1 です。

11.3 上三角行列

上三角行列 \(U\) の固有値は対角成分です。関数 \(f(U)\) も上三角で、対角成分は \(f(u_{ii})\) になります。三角方程式 \(Ux=b\) は後ろから

と解けます。計算量は \(O(n^2)\)。上三角は Schur 分解と LU の \(U\) に現れます。

11.4 三重対角行列

三重対角行列

に対する連立方程式は Thomas 法で \(O(n)\) です。これはピボットなし LU を三重対角構造に特殊化したものです。

Given lower a[2..n], diagonal b[1..n], upper c[1..n-1], rhs d[1..n]
for i = 2..n:
    m = a[i] / b[i-1]
    b[i] = b[i] - m c[i-1]
    d[i] = d[i] - m d[i-1]
x[n] = d[n] / b[n]
for i = n-1 down to 1:
    x[i] = (d[i] - c[i] x[i+1]) / b[i]

対称三重対角行列の固有値は Sturm 列で個数判定できます。\(T-\mu I\) の LDL^T 分解の符号変化数により、\(\mu\) 未満の固有値数を数えられます。これにより二分法で固有値を高信頼に求められます。

11.5 正規行列

正規行列 \(A^*A=AA^*\) はユニタリ対角化可能です。Hermitian、歪 Hermitian、ユニタリはいずれも正規です。正規行列では

が成り立ちます。非正規行列では固有値だけからノルムや過渡挙動を判断できません。

11.6 冪零行列・Jordan ブロック

\(N^k=0\) となる行列を冪零行列といいます。Jordan ブロック

は固有値 \(\lambda\) まわりの非対角化可能性を表します。行列関数では

しかし数値計算で Jordan 構造を直接求めるのは避けます。

12. 低ランク近似・疑似逆・最小二乗

12.1 最小二乗の三つの解法

12.2 Tikhonov 正則化

ノイズ増幅を抑えるため

を解くと、SVD 表現で

小さい特異値方向の係数が抑制されます。切断 SVD は \(\sigma_i\) が閾値以下の項を捨てる方法です。

12.3 PCA

データ行列 \(X\in\mathbb R^{m\times n}\) の列を中心化したとき、共分散行列は \(C=(1/(m-1))X^TX\)。SVD \(X=U\Sigma V^T\) により主成分方向は \(V\)、分散は \(\sigma_i^2/(m-1)\) です。

rank \(k\) 近似 \(X_k=U_k\Sigma_kV_k^T\) は Frobenius ノルムで最良なので、PCA は最小二乗的な次元削減です。

13. 行列関数・指数関数・極分解

13.1 行列多項式と解析関数

多項式 \(p(z)=\sum_k c_k z^k\) に対し

解析関数 \(f\) も冪級数や Schur 分解を通じて定義できます。\(A=QTQ^*\) なら

対角化可能なら \(A=X\Lambda X^{-1}\) により \(f(A)=Xf(\Lambda)X^{-1}\)。Jordan 形を使う場合は導関数が現れます。

13.2 行列指数関数

行列指数関数は

で定義され、線形微分方程式 \(x'(t)=Ax(t)\) の解は \(x(t)=e^{tA}x(0)\) です。計算には scaling and squaring + Padé 近似、Krylov 法、対角化・Schur 法が使われます。

13.3 極分解

任意の \(A\in\mathbb C^{m\times n}\) は

と分解できます。\(H=(A^*A)^{1/2}\) は Hermitian 半正定値、\(Q\) は部分等長写像です。列フルランクなら

SVD \(A=U\Sigma V^*\) から \(Q=UV^*\)、\(H=V\Sigma V^*\) と得られます。

14. 疎行列・グラフ・帯行列

14.1 疎行列の本質

疎行列とは非零成分が少ない行列です。ただし「少ない」だけでは不十分で、非零パターン、fill-in、メモリアクセス、並列性が重要です。疎行列を密行列として扱うと計算量・メモリが破綻します。

14.2 行列とグラフ

対称疎行列の非零パターンは無向グラフで表せます。\(a_{ij}\ne0\) なら頂点 \(i,j\) に辺を張ります。Cholesky 分解で新たに非零が発生する fill-in は、グラフの消去過程で辺を追加することに対応します。

14.3 並べ替え

疎 Cholesky/LU の fill-in を減らすため、対称問題では Approximate Minimum Degree、Nested Dissection などの並べ替えを使います。非対称 LU では列・行の並べ替えと数値ピボットの両立が課題です。

14.4 反復法と前処理

大規模疎行列では直接法より反復法が有効なことがあります。基本は Krylov 部分空間

正定値 Hermitian なら共役勾配法、非対称なら GMRES、BiCGSTAB などを使います。収束には前処理 \(M^{-1}A\) が決定的です。

Jacobi、SSOR、不完全 LU/Cholesky、代数的マルチグリッドなどがあります。良い前処理はスペクトルを集中させ、Krylov 法の反復回数を減らします。

15. 対話的デモ

以下はブラウザ上で動く簡易デモです。外部ライブラリに依存しない部分はオフラインでも動きます。Graphviz 図は d3-graphviz の読み込み後にレンダリングされます。

15.1 2×2 線形変換の可視化

15.2 小行列計算機

15.3 LU 分解ステップ

15.4 Power iteration

15.5 SVD 特異値デモ

16. 例題集

17. 演習問題

各問題の下にヒントを折りたたみで置いています。

18. 実装チェックリスト

19. 参考文献・次に読むもの

さらに深めるには以下の古典的文献が有用です。

• Golub and Van Loan, Matrix Computations.
• Trefethen and Bau, Numerical Linear Algebra.
• Horn and Johnson, Matrix Analysis.
• Demmel, Applied Numerical Linear Algebra.
• Higham, Accuracy and Stability of Numerical Algorithms.
• Saad, Iterative Methods for Sparse Linear Systems.