优化理论中的常见迭代。
Jacobi迭代
理论
对于线性方程组：Ax=b，若A非奇异且aij不等于0，i=1,2,3…n
若方程组可同解变型为：
则可以利用Jacobi迭代法计算
其中，每个迭代步的值为：
对于A可以进行分解：
例子
考虑如下方程组：
可以表示成如下形式：
根据上述理论描述，可以转化为Jacobi迭代求解：
例如对于初始点P0=(x0, y0, z0)=(1, 2, 2)，迭代一次后得到P1=(1.75, 3.375, 3.00)
多次求解后得到的迭代序列如下：
可以看出，经过多次迭代后，结果将收敛于(2, 4, 3),符合原方程组解的值。
代码(Python环境实现)
代码演示了上述例子中的迭代过程，其中利用三个函数模拟了迭代中的三个等式。设置当迭代前后值的差值小于1e-18时停止迭代。
import numpy as np
def getXValue(y, z):
return (7 + y - z) / 4
def getYValue(x, z):
return (21 + 4 * x + z) / 8
def getZValue(x, y):
return (15 + 2 * x - y) / 5
if __name__ == "__main__":
xData = []; yData = []; zData = []
xData.append(1); yData.append(2); zData.append(2)
itor = 1
print("迭代轮次 ", "X ", "Y ", "Z ")
while(True):
xTmp = getXValue(yData[-1], zData[-1]); yTmp = getYValue(xData[-1], zData[-1]); zTmp = getZValue(xData[-1], yData[-1])
if(abs(xTmp - xData[-1]) < 1e-18):
break
if(abs(yTmp - yData[-1]) < 1e-18):
break
if(abs(zTmp - zData[-1]) < 1e-18):
break
xData.append(xTmp); yData.append(yTmp); zData.append(zTmp)
print(itor, "--", xTmp, "--", yTmp, "--", zTmp)
itor += 1
代码输出结果如下：
迭代轮次 X Y Z
1 -- 1.75 -- 3.375 -- 3.0
2 -- 1.84375 -- 3.875 -- 3.025
3 -- 1.9625 -- 3.925 -- 2.9625
4 -- 1.990625 -- 3.9765625 -- 3.0
5 -- 1.994140625 -- 3.9953125 -- 3.0009375
6 -- 1.9985937500000002 -- 3.9971875 -- 2.99859375
7 -- 1.9996484374999999 -- 3.9991210937500004 -- 3.0
8 -- 1.9997802734375 -- 3.9998242187499997 -- 3.0000351562499996
9 -- 1.999947265625 -- 3.99989453125 -- 2.999947265625
10 -- 1.99998681640625 -- 3.999967041015625 -- 3.0
11 -- 1.9999917602539061 -- 3.999993408203125 -- 3.0000013183593746
12 -- 1.9999980224609375 -- 3.9999960449218745 -- 2.9999980224609373
13 -- 1.9999995056152344 -- 3.999998764038086 -- 3.0
14 -- 1.9999996910095215 -- 3.999999752807617 -- 3.0000000494384764
15 -- 1.9999999258422854 -- 3.9999998516845703 -- 2.9999999258422845
16 -- 1.9999999814605716 -- 3.9999999536514284 -- 3.0
17 -- 1.9999999884128572 -- 3.9999999907302857 -- 3.000000001853943
18 -- 1.9999999972190858 -- 3.9999999944381717 -- 2.9999999972190854
19 -- 1.9999999993047717 -- 3.9999999982619285 -- 2.9999999999999996
20 -- 1.9999999995654822 -- 3.9999999996523856 -- 3.0000000000695226
21 -- 1.9999999998957154 -- 3.9999999997914313 -- 2.999999999895716
22 -- 1.9999999999739289 -- 3.9999999999348224 -- 3.0
23 -- 1.9999999999837055 -- 3.999999999986964 -- 3.000000000002607
24 -- 1.9999999999960896 -- 3.9999999999921787 -- 2.9999999999960894
25 -- 1.9999999999990223 -- 3.9999999999975557 -- 3.0000000000000004
26 -- 1.999999999999389 -- 3.999999999999511 -- 3.000000000000098
27 -- 1.9999999999998535 -- 3.999999999999707 -- 2.9999999999998535
28 -- 1.9999999999999634 -- 3.9999999999999085 -- 3.0000000000000004
29 -- 1.999999999999977 -- 3.999999999999982 -- 3.0000000000000036
30 -- 1.9999999999999947 -- 3.999999999999989 -- 2.9999999999999942
31 -- 1.9999999999999987 -- 3.9999999999999964 -- 3.0
虽然最终结果并没有达到(2, 4 ,3)，但是就精度而言已经可以满足要求。
总结
有时，Jacobi迭代并不能收敛，可能会得到一个发散的序列：
对于上述算例，我们依旧采用之前的流程计算，得到的部分迭代序列为：
迭代矩阵M的所有特征值的模都小于1，是Jacobi迭代收敛的充要条件；若系数矩阵A按行（列）严格对角占优，则Jacobi迭代收敛。
严格对角占优矩阵：设有N×N维矩阵A，如果：
显然，严格对角占优的意思就是指对角线上元素的绝对值不小于所在行其他元素的绝对值和。
Gauss-Seidel迭代
理论
采用Jacobi迭代的算法如下：
在雅克比迭代法中，并没有对新算出的分量进行充分利用，一般来说，这些新算出计算的结果要比上一步计算的结果精确。对上式第二个方程组，第一行式子算出的x值立即投入第二行方程里，第二行式子的结果算出后投入第三行方程中，直到第n个方程。由于xk+1被认为是比xk更好的x之近似值，所以在计算yk+1时用xk+1来替换xk是合理的。同理，可用xk+1和yk+1计算zk+1。迭代过程修改为：
设迭代格式如下：
其中：
截取的一些定理如下：
例子
使用与Jacobi迭代中同样的例子。
代码
修改后的代码如下：
import numpy as np
def getXValue(y, z):
return (7 + y - z) / 4
def getYValue(x, z):
return (21 + 4 * x + z) / 8
def getZValue(x, y):
return (15 + 2 * x - y) / 5
if __name__ == "__main__":
xData = []; yData = []; zData = []
xData.append(1); yData.append(2); zData.append(2)
itor = 1
print("迭代轮次 ", "X ", "Y ", "Z ")
while(True):
xTmp = getXValue(yData[-1], zData[-1])
if (abs(xTmp - xData[-1]) < 1e-18):
break
xData.append(xTmp)
yTmp = getYValue(xData[-1], zData[-1])
if (abs(yTmp - yData[-1]) < 1e-18):
break
yData.append(yTmp)
zTmp = getZValue(xData[-1], yData[-1])
if (abs(zTmp - zData[-1]) < 1e-18):
break
zData.append(zTmp)
print(itor, "--", xTmp, "--", yTmp, "--", zTmp)
itor += 1
代码输出结果如下：
迭代轮次 X Y Z
1 -- 1.75 -- 3.75 -- 2.95
2 -- 1.95 -- 3.96875 -- 2.9862499999999996
3 -- 1.995625 -- 3.99609375 -- 2.9990312500000003
4 -- 1.999265625 -- 3.99951171875 -- 2.99980390625
5 -- 1.9999269531250001 -- 3.99993896484375 -- 2.99998298828125
6 -- 1.999988994140625 -- 3.9999923706054688 -- 2.999997123535157
7 -- 1.999998811767578 -- 3.9999990463256836 -- 2.999999715441894
8 -- 1.9999998327209474 -- 3.9999998807907104 -- 2.9999999569302367
9 -- 1.9999999809651183 -- 3.9999999850988384 -- 2.9999999953662795
10 -- 1.99999999743314 -- 3.999999998137355 -- 2.9999999993457847
11 -- 1.9999999996978925 -- 3.9999999997671694 -- 2.9999999999257234
12 -- 1.9999999999603615 -- 3.999999999970896 -- 2.999999999989966
13 -- 1.9999999999952327 -- 3.999999999996362 -- 2.9999999999988205
14 -- 1.9999999999993854 -- 3.9999999999995453 -- 2.999999999999845
15 -- 1.999999999999925 -- 3.999999999999943 -- 2.9999999999999813
16 -- 1.9999999999999905 -- 3.999999999999993 -- 2.999999999999998
17 -- 1.9999999999999987 -- 3.9999999999999987 -- 2.9999999999999996
18 -- 1.9999999999999996 -- 4.0 -- 3.0
总结
当矩阵A具有严格对角优势时，可证明高斯-塞德尔迭代法也会收敛。在大多数情况下，高斯-塞德尔迭代法比雅可比迭代法收敛得更快，因此通常会利用高斯-塞德尔迭代法。但在某些情况下，雅可比迭代会收敛，而高斯-塞德尔迭代不会收敛。
牛顿迭代
理论
设已知方程f(x)=0的近似根为x0，则在x0附近f(x)可用一阶泰勒多项式近似代替。
因此，方程f(x)=0,可以近似地表示p(x)=0。设f(x0)处的导数不为0，x1可以满足：
重复迭代，得到迭代公式：
不动点方程为：
例子
例1.求方程的根，f(x)=x*x-1,求解过程见代码1
例2.优于sqrt()的求根法，过程见代码2
例3.雷神之锤3代码黑科技(C++系)，见代码3
代码
代码1
def func(x):
# return 5 * x * x * x + 4 * x * x + 2
return x * x - 1
def func2(x):
# return 15 * x * x + 8 * x
return x * 2
def getRoot(x):
data = x
while (abs(func(data) > 1e-18)):
data = data - func(data) / func2(data)
# print(data)
print("最终结果为：", data)
if __name__=="__main__":
getRoot(1000)
getRoot(-1000)
代码2
def f(x, a):
return x * x - a
def f1(x):
return x * 2
def getRootByNowton(a):
x = 10
print("初始值为：", x)
while(abs(f(x, a)) > 1e-8):
x = (x + a/x) / 2
print("迭代值为：", x)
print(x)
if __name__=="__main__":
getRootByNowton(10)
代码3
int sqrt(float x) {
if(x == 0) return 0;
float result = x;
float xhalf = 0.5f*result;
int i = *(int*)&result;
i = 0x5f375a86- (i>>1); // what the fuck?
result = *(float*)&i;
result = result*(1.5f-xhalf*result*result); // Newton step, repeating increases accuracy
result = result*(1.5f-xhalf*result*result);
return 1.0f/result;
}
总结
是用于求解方根的黑科技！