[关闭]
@fsfzp888 2019-01-23T13:27:22.000000Z 字数 17128 阅读 1144

共轭梯度

数值优化 线性代数


本文介绍和总结共轭梯度(conjugate gradient)算法的如下内容:

共轭梯度算法的背景和由来
梯度算法的几何理解和共轭梯度算法的流程
基于SSOR预条件子的共轭梯度算法
非线性共轭梯度算法

背景

共轭梯度算法的由来,源于求解如下所示的二次型问题:

如果对于问题的求解,可以被转化为如上所示的二次型表达式(Quadratic Form)的极值,那么求解这个问题的方法也称为二次规划(Quadratic Programming)
二次规划问题,有一个极好的性质,这个性质直接导致了该问题的求解有较好的几何解释,这个性质就是如果矩阵为对称矩阵,那么把这个函数对列向量求导取极值,得到的梯度方向向量满足:

这个是经典的线性代数求解线性方程组的问题,求极值即是求解线性方程组,而矩阵的性质,直接决定了梯度方向的情况。
对一个从维到维的函数映射,对列向量求导即:


如果仅仅只看二次型中的二次项,那么

仔细观察上述表达式,可以发现,求导结果的第个元素等于矩阵的第行和第列分别与向量点积的结果,所以

如果A是对称矩阵的形式,那么,通过求导求解极值相当于求解方程组

二次规划是一些问题的最简单的建模表达形式,比如,在EDA(电子电气自动化)软件开发中,有这样的需求,就是要对各个电子元器件在电路板或者FPGA的电路平面进行自动化布局。每个元器件可以抽象为一些有面积的方型块,然后通过设定一些约束,定义损失函数,最后,通过不断迭代减少损失函数,更新各个电路元件的坐标信息,最终得到各个电子元器件在电路中的较优的全局布局结果。这些图可以有多层,而且损失函数通常也是非线性的,迭代求解通常基于各种梯度下降算法。期间需要考虑各个图层每一块面积的覆盖率,重叠程度等等参数。
如果对问题进行全面的建模,那么自然会比较复杂,如果只考虑最简单的情况,通过把元器件视为一个点,同时通过引入锚点,那么这个问题可以被简化为二次规划问题,通常二次规划问题用于电子元器件的预布局,然后才引入非线性函数进一步处理。
如果优化的目标可以被简化为让各个元器件之间的坐标之间的距离尽可能分散,那么此时损失函数等于

如果令,那么,完全可以表达成两个二次型表达式,其中后半部分的表示的是锚点到可动点的权重,实际优化该函数,锚点部分只会对的对角线,以及有影响,优化的结果就是在锚点的约束下,各个点的距离之间尽可能小。
由于问题的特殊形式,问题的求解完全也可以把左边分别分成两个部分:

也就是等价于求解线性方程组了

几何理解和共轭梯度算法的流程

共轭梯度算法来源于二次型问题的求解,二次型问题可以在几何上直观的理解,即使是最速梯度下降(Steepest Descent),也可以这样来理解。

二次型的图形

二次型的几何图解完全取决于矩阵,如果是正定的矩阵,那么图像将会是:

也就是是一个碗状朝上的结构,因为对于正定矩阵来说,任何非零向量,都有:


这个定义的本质是,所有的特征向量所对应的特征值都是正的,因为:

对于二次型来说,求导的梯度方向由决定,正定(positive definite),意味着各个向量变换到的列空间中得到的向量都是单调递增的,所以二次型表达式(如果是二维的)的几何图形就是中间凹陷,四周往上的碗状。同理,可以定义半正定,负定等,那么图形就是:

对于既非负定也非正定的情况,就是马鞍型,这个不管使用什么梯度方法,理论上都可能不会收敛。对于更高维度的情况,其实也是类似的,这里边正定或者负定的属性,决定了梯度方向是不是单调递增/递减的。对二次型,后边都是基于对称正定矩阵来分析的,那些梯度方法也源于此来进行基本的几何分析。

Steepest Descent

对于最速梯度下降(Steepest Descent)而言,选定了初始的方向向量后,后续每一次迭代所选择的方向向量由当前的梯度方向决定:


其中是方向向量,默认就是取梯度的负方向。这个过程的本质,其实就是线搜索(line search),如下图所示:

上图是二次型在变量维度为二维的时候的共势等高图,线是梯度负方向,每次迭代的本质就是在沿着这条线走,直到到达某点后停止,得到此时的向量,一般情况下,步长是可以固定的,特别是在函数是非线性,且非常复杂的时候,要想每一次都计算最好的步长,可能是非常复杂且不划算的,不过对于二次型问题而言,我们却可以设法找到最优的,通过查看上图的各个结果向量,为了更接近与最小点,应该和正交,此时,才是沿着当前的搜索方向的最优解。因为这个时候,向量已经沿着该方向尽可能往这个方向走,且已经不可能贡献再多的效用了,再走就又回去了。
对于二次型,梯度负方向等于:

利用正交关系进行推导,可以得到:

所以

这里边每次给出的时候都相当于直接求导了,如果根据迭代的形式给出下一个方向向量,那么

用迭代的形式给出每一次的梯度方向向量,虽然也许能够节省计算开销,但是限于实际机器的浮点精度,可能会产生累积误差,所以一般也需要在迭代了几十次后,重新严格计算梯度,以避免传递误差过大。
这样我们就可以得到最速梯度下降算法的流程是:

  1. import numpy as np
  2. def steepest_descent(A, b, x_initial, max_step, thresh=0.00001):
  3. assert(isinstance(A, np.matrix))
  4. assert (isinstance(b, np.matrix))
  5. assert (isinstance(x_initial, np.matrix))
  6. x = x_initial
  7. for _ in range(max_step):
  8. r = b - A * x
  9. alpha = (r.transpose() * r) / (r.transpose() * A * r)
  10. x = x + r * alpha
  11. dist = np.sqrt(np.sum(np.square(b - A * x)))
  12. if dist < thresh:
  13. break
  14. return x
  15. if __name__ == '__main__':
  16. N = 100
  17. Ar = np.mat(np.random.rand(N, N))
  18. As = Ar * Ar.transpose() # get positive definite matrix
  19. bn = np.mat(np.random.rand(N, 1))
  20. xi = np.mat(np.zeros((N, 1)))
  21. xr = steepest_descent(As, bn, xi, 1000)
  22. print('1000:', np.sqrt(np.sum(np.square(bn - As * xr))))
  23. xr = steepest_descent(As, bn, xi, 10000)
  24. print('10000:', np.sqrt(np.sum(np.square(bn - As * xr))))
  25. xr = steepest_descent(As, bn, xi, 20000)
  26. print('20000:', np.sqrt(np.sum(np.square(bn - As * xr))))

上边的Python代码体现了SD的最直接流程,梯度都是直接计算的,实际真的要应用,还要考虑稀疏矩阵,计算效率问题,否则效率太低了,也没有使用价值。当然,这个只用于二次型问题。
对一些非线性函数问题,就需要通过别的方式求导,而且步长也需要确定,或者设定个固定值!如下所示是该代码的一次运行结果:

  1. 1000: 1.9453821317304896
  2. 10000: 0.9116628408830175
  3. 20000: 0.6950563935826073

Conjugate Gradient

运行SD的例子,会发现,随机生成的一个对称矩阵,迭代了几万次仍然不收敛,可以见得,按照SD的理论公式求解效果有些不理想,收敛过慢。在数值优化里边,一个矩阵迭代求解方法是否收敛,取决于这个矩阵是否严格对角占优,即


越是严格对角占优,越是收敛得越快。
SD这个方法,使用梯度方向,每次选择梯度方向都和上一次的正交,这样迭代的结果取决于初始点的选择和,如果初始点不理想,可能一直在锯齿状的接近极值点,但是迭代多次后仍然不能命中。
对于当前的二次型而言,不是多么复杂的问题,变量的维度也许就是而已,存不存在这样的迭代方法,使得每一次迭代,都消减一个维度方向,后边的方向和前边的方向都完全正交呢?这样的话,我只需要步迭代就可以直接命中极值点了。这个思想,我觉得就是CG背后的想法。如果选择个线性无关的向量,那么依据施密特(Gram-Schmidt Process)正交化的过程,每个方向向量为:

Conjugate Direction

在CG中,引入了共轭方向的概念。这个概念本质上就是对正交的推广,任何向量在欧几里得空间内正交意味着:


如果引入矩阵,且

那么称向量正交的,这里边都是方阵,看起来就像是把变量线变换到的行列空间中,然后这些向量在的行列空间中正交。这些两两正交的向量被称为共轭方向向量。如果给个图解,那么就会如下图所示:

前后两个方向向量是正交(A-orthogonal)的,意味着,即使这两个向量在当前空间不正交,在的空间中正交也行,可以想象,当前空间中的所有向量都变换到的空间后,着两个向量在空间中就是正交的。

二次型的标准化

上边的过程容易让人联想起线性代数中的二次型标准化过程,二次型标准化,首先通过平移变换,让二次型变成了如下的形式:


然后找寻个标准正交的特征向量,然后做如下变换:

最终变成了对角矩阵的形式,这个时候对应的二次型的几何图形就是很标准的碗状了
个人觉得,正交也是希望在标准化后的空间中正交。其实标准化前后,所代表的线性空间是不变的。在几何理解里边,都以标准化后的空间中的向量来分析理解。

CG的推导过程

二次型的基本迭代流程还是不变的:


如果采用迭代求解梯度,那么下边的等式仍然成立:

对等式22向求导:

可以看到,最优的步长是使得下一个误差方向和当前的方向正交,由于

所以

最优的步长因子和当前的梯度和方向向量均相关。从上边的表达式知道,当前反向和剩余误差方向是正交的,而我们也希望:

如果下一个方向向量和当前的也正交就好了,当然了,如果下一个方向向量刚好选择到了剩余的误差向量,那么下一次迭代就会直接命中极值点了,但是一般这个是不可能的,因为是依赖于的,会有循环求解的问题,除非我们已经知道解了,否则我们是得不到的,而且在高维空间中,和向量正交的向量是有很多个的,我们希望的是找到一个满足上边表达式即可,于此同时,也满足

这样每次迭代都减少了一个维度,那么,最多步就可以收敛了。
这个过程,在变量为二维的情况下理解,大概如下图所示:

这些椭圆是三维图形在二维的等势面,其实就是相等的值所构成的一个椭圆。
一开始选择梯度作为初始方向,也就是等势椭圆的切线,下一步,由于取正交关系,是同向的,两步就命中结果了。对应到的空间中,这个椭圆就是个圆。变量是四维的情况,也可以用下图表示:

变量为三维的情况,加上因变量实际上是四个维度的,由于二次型也可以描述能量关系,所以相等能量的各个三维自变量构成的图形在三维就是个椭球体。映射到的空间中其实就是个标准的球型。每次迭代搜寻的方向向量,都是和过去所有的方向向量正交的,三维的情况可以看得比较明显,那就是,也就是方向向量在高维度的情况下,即使和过去的方向向量都正交,但是也不能说这个方向向量就是当前的错误向量。但是从事实出发,错误向量随着每次迭代,都减少了各个方向向量的那部分,所以,错误向量即,一开始,初始点到极值点的错误向量可以表示为:

而等式28中又限定了各个方向向量的关系,所以可以:

所以

上边的推到也证明,只要我们满足条件

条件1:各个方向向量之间两两正交
条件2:通过求导让每一次迭代的步长因子取得最优

那么每一次迭代求解的过程,其实就是相当于减少了这个方向向量上的误差,每一步的误差向量可以被表示为方向向量的连加的形式,这样只需要迭代步后就会直接收敛到极值点。这个和我们在几何上边的理解是完全一致的。
只是,我们该如何选择这些方向向量呢?
如果仅仅只是随便选择一组线性无关的构造基,然后如同等式17那样构造方向向量,并满足正交,那么:


通过这一种方式,也可以得到如何去构造所有的方向向量,只是,任何方向向量都和过去的个向量相关,而且每次迭代还需要求解系数,这样就直接导致了计算所需的内存和时间开销巨大,没有实用价值。所以这一组构造基不应该随便选。而应该让等式32得出的各个系数拥有良好的关系,减少计算所需的开销。
事实上,共轭方法,在一开始提出的时候,确实存在实现上的困难,不过通过仔细查看这些关系,我们还没有把梯度向量引入。
在计算过程中,我们始终要得到当前点的梯度,在上边两个条件满足的情况下,观察错误向量,梯度等之间的关系:

也就是说,任何迭代后期的梯度都是和之前的方向在当前空间正交的,进一步的

结合前边的所有关系,可以得到如下图所示的关系图:

图中所示,的投影是相等的,所以末端所在的平面和以及所在的平面共面。通过上边的各个关系可以知道,没有必要取选择别的基,直接让就是最优的选择,这个时候,各个向量的关系图将会变成如下所示:

这样做的直接结果有:


综合前边所属,CG的流程可以被表示为:


for loop until reach max_step or cost_function below threshold:





return

使用Python代码实现:

  1. import numpy as np
  2. def conjugate_gradient(A, b, x_initial, max_step, threshold=0.00001):
  3. assert(isinstance(A, np.matrix))
  4. assert(isinstance(b, np.matrix))
  5. assert(isinstance(x_initial, np.matrix))
  6. r_old = b - A * x_initial
  7. d = r_old
  8. x = x_initial
  9. for i in range(max_step):
  10. alpha = (r_old.transpose() * r_old) / (d.transpose() * A * d)
  11. x = x + d * alpha
  12. # r_new = b - A * x
  13. r_new = r_old - A * d * alpha
  14. beta = (r_new.transpose() * r_new) / (r_old.transpose() * r_old)
  15. d = r_new + d * beta
  16. r_old = r_new
  17. cf = np.sqrt(np.sum(np.square(b - A * x)))
  18. if cf < threshold:
  19. print("Using step: ", i)
  20. break
  21. return x
  22. if __name__ == '__main__':
  23. N = 200
  24. Ar = np.mat(np.random.rand(N, N))
  25. As = Ar * Ar.transpose()
  26. bn = np.mat(np.random.rand(N, 1))
  27. xi = np.mat(np.random.rand(N, 1))
  28. xr = conjugate_gradient(As, bn, xi, 1000)
  29. print('1000:', np.sqrt(np.sum(np.square(bn - As * xr))))
  30. xr = conjugate_gradient(As, bn, xi, 10000)
  31. print('10000:', np.sqrt(np.sum(np.square(bn - As * xr))))

运行上述程序的一个输出:

  1. Using step: 410
  2. 1000: 6.43707958798e-06
  3. Using step: 410
  4. 10000: 6.43707958798e-06

大约需要两倍于维度的迭代次数才可以收敛到要求的范围,虽然理论上100维的变量值需要100次迭代,但是实际上各种浮点数运算误差导致不可能实现这样理想化的结果。通过对比会发现,迭代效果要比SD要好。当然了,实际使用的时候,判断损失函数其实没必要通过那样的形式。而且可以几十次迭代后才计算一下等等。

基于SSOR预条件子的共轭梯度算法

方程组的迭代求解

上边的迭代算法,本质上其实还是在求解形如的方程。毕竟最终的结果,都是让尽可能的接近于零。一般情况下,对于的求解,通常直接通过求逆得到,这就存在求矩阵的逆可能不存在等等问题。如果不是方阵,且,那么就是在求解超定方程组,方程组可能没有解,求解的过程就会等效于最小二乘,因为,为了让误差最小,其实就是让最终的的解得到的最近,那么只能是的列空间的投影,所以它们之间的差和的列空间正交,所以有


即使是非方阵,但是却是对称方阵,所以求解逆是可行的,这个时候其实还是转换成了

所以还是可以使用上边的梯度迭代方法来求解。所以可以见得,线性系统问题,即可以表达为求解的系统问题,还是比较多的。
求解的方法,除了直接求逆,以及高斯的LU分解,其它的基本就是迭代方法了。最早的迭代有高斯-赛德尔迭代,SOR等,它们都是首先把拆分:

也就是把矩阵拆解成对角,上三角,下三角的形式,对于高斯-赛德尔方法,把变成了:

提取成对角,上下三角矩阵的原因,恐怕还是因为逆矩阵好算。上边的过程,每次迭代其实都是循环回代的过程,所以等式40右边有,其实是说要用当前解出来的那些新的变量值的意思。
对于SOR来说,则引入了松弛系数,迭代过程是

上边的等式其实是用了转换而来。
这些方法都有其收敛条件,在数值优化中有分析。
对于稀疏矩阵,或者是对称正定矩阵来说,还是一般使用本篇总结的梯度迭代方法。只不过值得留意的却是,竟然对应于一个二次型的极值求解过程

Precondition

在数值分析里边,矩阵迭代求解的收敛速度取决于矩阵条件数:


其中是矩阵范数,数值分析里边研究了矩阵迭代的收敛很大程度上取决于条件数,而预条件方法就是为了减少条件数,快速收敛。详细的分析很复杂,详情可看数值分析教材,这里只阐述用于共轭梯度算法的预条件过程。
预条件是希望找到一个,将问题转化为求解:

其中的可逆矩阵,称为预条件矩阵。在改进的共轭梯度方法中,这个矩阵是对称正定的。这个矩阵试图对逆转,从而让的条件数降低。此外,还引入了一种广义内积来取代欧几里得内积,所有内积默认都改成内积,原始的共轭梯度方法仍然成立,因为矩阵相对于新的内积仍然是对称正定矩阵:

那么,由于CG的推到基本都是基于的形式,如果把等式44的内积结果应用到所有推导的过程,就有:

这里边一个有意思的地方在于,这个方法,把的欧几里得的内积,替换成了,定义内积在的逆空间里边,因为内积仍然满足普通内积该有的性质,所以可以完全替代欧几里得内积,公式仍然不发生变化。
进一步转换可以得到:

对于SSOR(对称连续过松弛)预条件子来说:

总结上边的过程可以得到,SSOR预条件子的共轭梯度算法的流程是:


for loop until reach max_step or cost_function below threshold:






return

对应的Python代码:

  1. import numpy as np
  2. def get_ssor_precondition_matrix(A, w):
  3. UD = np.triu(A)
  4. LD = np.tril(A)
  5. dim = A.shape[0]
  6. D = np.mat(np.zeros((dim, dim)))
  7. for i in range(dim):
  8. D[i, i] = A[i, i]
  9. for i in range(dim):
  10. for j in range(i+1, dim):
  11. UD[i, j] = w * UD[i, j]
  12. for i in range(dim):
  13. for j in range(0, i):
  14. LD[i, j] = w * LD[i, j]
  15. # 对上下三角矩阵求逆矩阵,其实不必用通用的求逆方法,不停回代即可
  16. return np.mat(np.linalg.inv(UD)) * D * np.mat(np.linalg.inv(LD))
  17. def precondition_conjugate_gradient(A, b, x_initial, max_step,
  18. threshold=0.00001, w=0.2):
  19. assert(isinstance(A, np.matrix))
  20. assert(isinstance(b, np.matrix))
  21. assert(isinstance(x_initial, np.matrix))
  22. r_old = b - A * x_initial
  23. M_inv = get_ssor_precondition_matrix(A, w)
  24. z_old = M_inv * r_old
  25. d = z_old
  26. x = x_initial
  27. for i in range(max_step):
  28. alpha = (z_old.transpose() * r_old) / (d.transpose() * A * d)
  29. x = x + d * alpha
  30. # r_new = b - A * x
  31. r_new = r_old - A * d * alpha
  32. z_new = M_inv * r_new
  33. beta = (z_new.transpose() * r_new) / (z_old.transpose() * r_old)
  34. d = z_new + d * beta
  35. r_old = r_new
  36. z_old = z_new
  37. cf = np.sqrt(np.sum(np.square(b - A * x)))
  38. if cf < threshold:
  39. print("Using step: ", i)
  40. break
  41. return x
  42. if __name__ == '__main__':
  43. N = 200
  44. Ar = np.mat(np.random.rand(N, N))
  45. As = Ar * Ar.transpose()
  46. bn = np.mat(np.random.rand(N, 1))
  47. xi = np.mat(np.random.rand(N, 1))
  48. xr = precondition_conjugate_gradient(As, bn, xi, 1000, 0.00001, 0.05)
  49. print('w=0.05:', np.sqrt(np.sum(np.square(bn - As * xr))))
  50. xr = precondition_conjugate_gradient(As, bn, xi, 1000, 0.00001, 0.5)
  51. print('w=0.5:', np.sqrt(np.sum(np.square(bn - As * xr))))
  52. xr = precondition_conjugate_gradient(As, bn, xi, 1000, 0.00001, 1)
  53. print('w=1:', np.sqrt(np.sum(np.square(bn - As * xr))))

运行三次的结果:

  1. runfile('C:/Users/zczx1/.spyder-py3/temp.py', wdir='C:/Users/zczx1/.spyder-py3')
  2. Using step: 378
  3. w=0.05: 6.71578034331e-06
  4. Using step: 405
  5. w=0.5: 9.67223926338e-06
  6. Using step: 573
  7. w=1: 8.09438315554e-06
  8. runfile('C:/Users/zczx1/.spyder-py3/temp.py', wdir='C:/Users/zczx1/.spyder-py3')
  9. Using step: 371
  10. w=0.05: 6.1508910381e-06
  11. Using step: 401
  12. w=0.5: 8.51261715479e-06
  13. Using step: 602
  14. w=1: 7.57385104633e-06
  15. runfile('C:/Users/zczx1/.spyder-py3/temp.py', wdir='C:/Users/zczx1/.spyder-py3')
  16. Using step: 373
  17. w=0.05: 6.44798434081e-06
  18. Using step: 406
  19. w=0.5: 6.55294163482e-06
  20. Using step: 580
  21. w=1: 8.13783429322e-06

虽然随机生成的矩阵不太能说明问题,但是这个只能说,在松弛系数比较小的时候,PCG和CG相比,相对来说减少了迭代的次数。但是当松弛系数较大的时候,还不如CG。

非线性共轭梯度算法

算法流程

上边的总结,都是在二次型问题上进行的分析。可以看到,虽然算法本身的实现比较简单,但是算法背后的数学思路如果进行严谨地推导的话,还是有些复杂的,不过因为有了几何的理解,所以也还是比较直观。
但是如果这个问题的表达式函数是非线性的,共轭梯度的思想仍然可以被推广,只是形式上不那么精确。而事实上,就算严格按照上边的算法流程,通过上边的实践也会发现,由于round off误差,其实不能达到理论的计算效果,所以在某种程度上近似求解也是可以的。回顾共轭梯度算法和核心,其实就两个:
1. 依据步长和方向向量更新解向量,即
2. 利用梯度更新方向向量,即
其中,,我们只需要考虑怎么选择即可
,是希望找到这样的步长因子使得,本质上还是希望下一个梯度方向和当前的方向向量正交,这个时候应该就是最优的解,和SD中的线搜索是类似的。不过,在实际的实现中,可能只是设置一个固定值,然后乘以一个比例,这个比例随着迭代的进行可能会变得越来越小。所以:


而对于,不同的学者提出了使用不同的公式,比如:


也还有很多其它的计算方法,不过显然,这些只是一种近似求解的方法。因为非线性的损失函数各种各样,精确的求解不现实。这样,非线性CG的流程如下:


for loop until max step or :




非线性函数中最优的选择

寻求最优的步长因子的过程其实就是广义线搜索,根据泰勒级数展开,有:



从等式52中,可以看到二次型中求解最优的影子,如果是二次型,那么就是的形式了,不过对于非线性问题来说,就需要计算二阶导数。二阶导数为:

其实就是求解Hessian矩阵。这个计算量可谓真的大,而且随着迭代,这个矩阵的具体数值,可能还是需要不断计算的,这个开销几乎是不可忍受的,所以严格的计算步长因子基本是不可取的。何况有些时候,可能都求不出二阶导,也无怪乎实际应用中,直接给个固定步长乘以比例了事。与其算这个矩阵那么复杂,不如多迭代几次。

非线性函数求导简介

非线性CG需要对函数进行求导,而求导的方法,目前有:
1. 公式法
2. 符号微分法
3. 自动微分法
这些根据具体场景来实现,只是非线性CG迭代流程就是如上所示的过程。
值得注意的是,非线性函数的这个过程,最终收敛到的可能只是局部极值,和初始值的选择相关。如果初始值在一个类二次型的区域,那么就可以收敛到这块区域的极值点。

总结

本文总结了本人学习EDA软件中布局算法中使用到的共轭梯度相关的内容。共轭梯度算法,应用到实际工程当中,可能还是使用非线性的部分内容,只是由于也是近似的求解,所以和随机梯度下降相比也没有太多的优势,毕竟都是和初始点相关,所以引入更多的随机性,也许可以取得更好的效果!

参考资料

1. An Introduction to the Conjugate Gradient Method Without the Agonizing Pain
2. Numerical Analysis

@fsfzp888
2019 年 01月

添加新批注
在作者公开此批注前,只有你和作者可见。
回复批注