基本问题
数值线性代数主要包括以下三个问题
- 解线性方程组的问题。给定 $Ax=b$，求 x
- 线性最小二乘问题。$\arg\min\limits_{x\in R^n} \mid\mid Ax-b \mid\mid_2$
- 求A的部分或全部 特征值、对应的 特征向量
数值线性代数的 必要性：线性代数的一些定理，证明过程非常漂亮，但对应的计算量大的惊人
设计算法的 主要技巧：矩阵分解。
敏度分析与误差分析，误差两个来源：原始数据本身的误差（可能来自测量）、计算过程产生的误差。所以需要两个理论
- 敏度分析
- 误差分析
算法复杂度和收敛速度。
- 复杂度分析：90年代之前只考虑乘除，90年代后的计算机加减乘除速度相差无几，所以都考虑
- 不能完全依据运算量来评价算法快慢，现代计算机运算速度远高于数据传输数独，所以算法的速度很大程度上取决于数据传输两。
- 对于迭代算法，还需要分析收敛速度。$\mid\mid x_k-x \mid\mid \leq c\mid\mid x_{k-1}-x \mid\mid^n$ 称为n次收敛，显然n越大越好
线性方程组消元法
Gauss消元法，适用于
- 中小规模的线性方程组，阶数不超过1000
- 系数矩阵稠密
- 没有任何特殊结构
思路：解方程 $Ax=b$，
- 把A分解为 $A=P^{-1}LU$,其中P是排列矩阵，L是下三角矩阵，U是上三角矩阵，
- 就可以把问题分解为：$Ly=Pb,Ux=y$，
- 三角矩阵对应的解法大为简化。
上面算法的问题：如果某主元不为0但很小，在计算机精度下的计算结果差很多。 解决：$PAQ=LU$
对于正定矩阵，有更好的办法：
若 A 正定，那么纯在下三角矩阵 L，使得 $A=LL^T$（Cholesky分解定理）
分块矩阵：可以分块处理，大大节省内存和计算量
计算特征值
幂法
求最大特征值的方法
- 定义要求特征值的矩阵 A，大小为 n × n。
- 设定一个 n 维的初始向量 x0，可以将其设置为随机值。
- 对 x0 进行归一化处理，即将向量除以其长度，得到单位向量 x。
-
通过以下迭代过程，更新向量 x ：x = A x / A x -
其中， A x 表示向量 A x 的长度。
-
- 计算向量 x 和 A x 的内积，即 xT A x。
- 计算特征值的近似值 λ = xT A x。
- 重复步骤 4 至步骤 6，直到 λ 收敛到某一精度。
- 将得到的特征值 λ 作为 A 的一个特征值，并将对应的单位向量 x 作为特征向量。
- 通过正交化过程，得到 A 的所有特征向量。
需要注意的是，当矩阵 A 不对称时，需要使用雅可比方法或 Householder 方法进行特征值计算。此外，在实际计算中，可能会存在计算误差、收敛速度缓慢等问题，需要进行调优。
在实际应用中，为了提高计算精度和加速收敛速度，可以对幂法进行改进，例如用反迭代法（Inverse Iteration）或QR分解迭代法（QR Iteration）来代替简单的幂法。
代码：
import numpy as np
def get_eig(A):
# 精度
loc = 1e-8
x = np.random.rand(4, 1)
lam_pre = 0
for i in range(1000):
tmp = A @ x
x = tmp / np.linalg.norm(tmp)
lam = x.T @ A @ x
if abs(lam_pre - lam) < loc:
break
lam_pre = lam
# 特征值，特征向量，迭代次数
return lam[0, 0], x, i
# 这里为了生成一个对称矩阵
A_ = np.random.rand(4, 4) * 255
A = np.sqrt(A_ @ A_.T)
# 使用数值方法
lam, x, i = get_eig(A)
# 使用 numpy 中的方法
eig, eig_vec = np.linalg.eig(A)
# 验证特征值正确
assert abs(lam - eig[0]) < 1e-7
# 验证特征向量正确。因为允许符号相反，都先调整为正数
tmp1 = eig_vec[:, 0] if eig_vec[0, 0] > 0 else -eig_vec[:, 0]
tmp2 = x[:, 0] if x[0, 0] > 0 else -x[:, 0]
assert (np.abs(tmp1 - tmp2)).max() < 1e-6
如果把循环终止的条件改为：x 不怎么变化了，性能还能提升
from numba import jit
@jit(nopython=True)
def get_eig(A):
# 精度
precision = 1e-8
x = np.random.rand(4, 1)
y_norm_old = 0
for i in range(1000):
y = A @ x
y_norm = np.linalg.norm(y)
x = y / y_norm
if -precision < y_norm - y_norm_old < precision:
# 特征值，特征向量
return (x.T @ A @ x)[0, 0], x
y_norm_old = y_norm
return (x.T @ A @ x)[0, 0], x
QR分解法
QR分解法是一种基于矩阵分解的方法，可以计算矩阵的所有特征值及其对应的特征向量。它的基本思想是将矩阵分解为一个正交矩阵和一个上三角矩阵的乘积，再不断迭代地对上三角矩阵进行QR分解，直到分解的上三角矩阵趋近于一个对角矩阵。
算法步骤：
1）给定矩阵$A$和一个初始的正交矩阵$Q^{(0)}$。
2）计算$A^{(0)}=Q^{(0)}A(Q^{(0)})^{-1}$，得到矩阵$A^{(0)}$的QR分解$A^{(0)}=Q^{(1)}R^{(1)}$。
3）继续迭代计算$A^{(k)}=(Q^{(k)})^{-1}A(Q^{(k)})$和$A^{(k)}=Q^{(k+1)}R^{(k+1)}$，直到得到的上三角矩阵$R^{(k)}$趋近于对角矩阵。
4）得到矩阵$A$的所有特征值为$R^{(k)}$的对角线元素，对应的特征向量为$Q^{(k)}$的列向量。
参考文献
《数值线性代数》北京大学出版社，徐高方
