数值代数思想

数值代数

第一章

舍入误差 数值解和精确解的距离。由于计算机存储理解浮点数，在进行一些计算后会丢失一些精度。很多基本运算，如加减乘除，可以保证误差控制在一定范围内。但多次执行较复杂的计算，误差就会累加，影响解的精度。因此，需要对合适的问题选择合适的算法。
浮点运算的基本模型是
fl(a.b) = (a.b)(1+delta)

第二章

解线性方程组，数值代数中最基本的问题。设 A 是 nxn矩阵，x和b是n维向量，要求解 Ax=b。在数学上易知方程有解当且仅当A非奇异，此时解为x=A-1b。但在数值上，一来，如何求解A的逆成了一个问题，二来，这是否是最快的求解方法。

已知最快的求解法是LU分解。但数值稳定性一般。算法的步骤是先将A分解成A=LU的形式，其中L是下三角矩阵，U是上三角矩阵。然后依次求解Ly=b和Ux=y。对上三角型的线性方程组，可以直接求解
xn=bn/Lnn, xn-1=(bn-1-ln-1nxn)/Ln-1n-1, …
总共需要 n(n-1)/2次减法，n(n-1)/2次乘法，n次除法，(访问了n(n+1)次元素)，因此一共需要O(n2)量级的flop。

那么自然的问题是，求解LU分解需要多少flop呢？我们先看看如何求解LU分解。我们假定L的对角元都是1，因此若LU分解存在，就是唯一的。

for i = 1 : n
uii = aii
for j = i+1 : n
lji = aji / uii
uij = aij
for j = i+1 : n
for k = i+1 : n
ajk = ajk - lji * uij

注意到在进行到第i步的时候，内层循环没有修改aij和aji的值，因此我们其实可以把L和U存储在矩阵A的内部。

for i = 1 : n
for j = i+1 : n
aji = aji / aii
for j = i+1 : n
for k = i+1 : n
ajk = ajk - aji * aik

在上述算法运行完以后，A的严格下三角区域就存储了L的值(不包括对角线)，而上三角区域存储了U的值。

LU分解的灵感来自于高斯消去法。对于方程 Ax=b，通过不断地对增广矩阵 [A b]做行变换，当A变换为单位阵时，b被变换为方程的解。

尽管想法很美好，但LU分解并不是对任何矩阵都存在。在上述算法中，要求进行到第i步时的aii不等于0。但是若A非奇异，第i列的后n-i个元素一定有非零的，因此，我们思考能不能通过行变换，将非零元放在aii的位置，使得LU分解可以继续进行。改进的算法叫做PLU分解。

p = (1:n)
for i = 1 : n
imax = max(a(i+1:n, i))
p
uii = aii
for j = i+1 : n
lji = aji / uii
uij = aij
for j = i+1 : n
for k = i+1 : n
ajk = ajk - lji * uij

现在来分析一下LU分解的复杂度。主要的计算量就在内层循环里，第i步需要 2i2次运算，一共需要大约 2/3 n3次 flop。

QR分解
速度很快的LU分解在某些情况下的确会导致精确性不足，这时就是要一些更稳定的算法，可想而知也会稍微慢一点，只要在能接受的范围内。接下来介绍的就是QR分解。QR分解的基本思路是通过左正交变换，将矩阵化为上三角矩阵。那么自然的问题是，要用什么样的正交变换？主要有两种，Housholder变换和Givens旋转。

Housholder变换形如 I-2uuT，u是单位向量。它的几何意义是以垂直于u的超平面为轴将向量镜像翻转。

Givens旋转的矩阵是两行两列的旋转矩阵，其他地方和单位阵一样。

一般利用Housholder变换做QR分解，那么将一个向量a变成e1，要怎么找Housholder变换的u呢？令v=a-e1，令u=v/|v|即可。但是在v的形成过程中，出现了大数减大数，容易造成数值不稳定，所以当a和e1接近时，令v=a+e1，将a变为-e1是更好的选择。

当QR分解做完后，矩阵Q并没有显式形成，而是一列Housholder变换，每个Housholder变换的存储需求是n个元素。做一次Housholder变换需要 3i2次flop，因此QR分解一共需要大约n3次flop。比LU分解的2/3n3次flop慢一点，但数值上稳定得多。

我们还没有用到Givens旋转，那它有啥用呢？当QR分解进行到最后一步时，只需再消去一个元素，这时用Givens旋转比Housholder变换速度更快。

如何利用分块矩阵乘法加速LU分解
LU分解的速度固然很快，但是用块矩阵乘能在计算机实际实现的时候表现更好。这是因为CPU一次只能把极少量的数据读取到cache中，cache内的数据被反复使用的速度是很快的。但如果同一个数据需要在不同的时刻被读取，那效率就不高了。我们想把LU分解变成块矩阵乘的形式，把矩阵作为算法里的每个元素。

for i = 1 : n
aii = liiuii (LU分解)
for j = i+1 : n
lji = aji uii -1
uij = lii-1 aij
for j = i+1 : n
for k = i+1 : n
ajk = ajk - lji * uij

第三章 最小二乘问题
在解完线性方程组Ax=b以后，现在面临着更复杂的情况，矩阵A的列数大于行数，因此方程没有解，然而我们还是想找一个x，在某种程度上充当最优解。这个想法有实际意义，比如线性回归问题，y=Xw+e，想找w使得|e|=|y-Xw|最小。

若A满秩，线性最小二乘问题有唯一解析解，由法方程给出。若A不满秩，解x属于一个线性空间，有唯一的最小范数解。

对于满秩的最小二乘问题，可以用QR分解。

QR分解算法LS问题
A=QR
x = R-1R-TRTQTb

其中Q不用显式地表达出来，只需用一列Housholder变换表示即可。

第四章 求一般矩阵的特征值

有代数基本定理，任何矩阵至少有一个复特征值。解特征值问题也就是求解 Ax=sx，其中 s和x都是未知的。一共有n+1个未定元。

数学上，每个方阵A都有Schur分解 A=QTQT，其中Q是正交阵，T是上三角阵。如果已求得Schur分解，那么T的对角线就是特征值。如何从特征值求特征向量呢？

事实上，可以从Schur标准型求特征向量。

我们称 将矩阵A变为QAQT的变换称为正交相似变换，其中Q是正交阵。

若矩阵的特征值是特征多项式的多重根，则对A微小的扰动会对特征值s造成很大的扰动。否则，设s是A的特征值，x和y是左右特征向量，设
(A+dA)(x+dx)=(s+ds)(x+ds)
则
dAx+Adx=ds x+sdx
y. dAx=ds y.x
ds = y.dAx / y. x

因此ds的扰动是dA量级的。
若A非奇异，s的相对扰动估计
ds/s <= kappa(A) dA/A

现在，讨论将矩阵化为Schur标准型(从而完成特征值和特征向量求解)的方法。

想要直接将矩阵化为 QTQT，可能没有直接法，因为方程 Ax=sx是一个非线性方程。但是，将A约化成稍微好一点的形式是可能的。