Skip to content

Latest commit

 

History

History
160 lines (124 loc) · 7.8 KB

cg.md

File metadata and controls

160 lines (124 loc) · 7.8 KB

Conjugate Gradient/共轭梯度

背景

前面提到了最速下降法,最速下降法很靠人品,如果某一轮的梯度能够和A的任一特征向量共线,那么立刻收敛,如果没有,可能要走很长的zig-zag路线,其中会有多轮优化选择近似的梯度。

针对这个问题,我们一定会想,能不能设计一种算法,一次性把某一个方向的残差走完?这就是共轭梯度要解决的问题。

共轭

共轭这个词真的很难从字面理解,所以我们慢慢引入这个概念。

前面的Steepest Descent里面,每一轮的梯度和上一轮的梯度正交,现在根据上面的想法,我们不希望同一个方向重复走多次,那么经过某一轮更新的方向d[i], 应该与后面的残差正交,也就是说:

# 目标函数,求x
A * x = b
# 前面的一些函数
def r(iteration):
	return b - A * x[iteration]

def e(iteration):
	return x[iteration] - x # x是ground truth

d[i] * e(i+1) = 0

根据这条,我们可以进行一定的推导,得到关于步长的公式:

d[i] * e(i+1) = 0
=> d[i] * (e(i) + step[i] * d[i]) = 0
=> step[i] = -(d[i] * e[i]) / (d[i] * d[i])

步长公式求出来了,但是这没有什么用,因为e[i]是没有办法知道的。所以利用正交的方法不可行。

于是乎,我们采用了另外一个方法,就是关于A正交,公式是

# 如果说x,y关于A正交,那么
x * A * y == 0 # 这里不考虑向量的转置问题

这是什么意思呢?过去的正交是两个向量的点乘为0,现在相当于先把向量y根据A进行线性变换,变成A的column space下的线性组合,再与x做点乘,结果正交。从几何意义上,有点像把y投影到A上,然后与x正交。这种正交就叫共轭

有了这个性质会有什么好处么?看看上面的式子:

step[i] = -(d[i] * A * e(i)) / (d[i] * A * d[i])
=> step[i] = (d[i] * r(i)) / (d[i] * A * d[i])

这里用到了一个变换,不熟悉的可以看Steepest那章。可以看到,当利用了共轭之后,e(i)被换成了可以求解的r(i),问题被解决了。

共轭的性质

性质1: 如果d[i]和e(i+1)共轭,那么求得的x[i+1]在d[i]的方向上是最小值

这个性质说了什么呢?在Steepest Descent中,我们证明了每一步的梯度都是与上一步正交,因为每一步都是选择了当前最优的点,也就是说沿着梯度方向,找到了这个方向上的最小点。那么我们刚才提到了共轭梯度的方法,能不能做到每一步找到当前最小的点呢?答案是肯定的,下面证明:

# 如果x[i+1]是d[i]方向上的最小点,那么盖点关于步长的导数应该为0f(x[i+1])/∂(step[i]) = 0 # 使用链式法则f(x[i+1])/∂(x[i+1]) * ∂(x[i+1])/∂(step[i]) = 0 
# 两部分分别变换,前面是f的导数,也就是残差,
# 后面x[i+1]=x[i] + step[i] * d[i], 对其求导
r(i+1) * d[i] = 0 # 继续利用残差公式
d[i] * A * e(i+1) = 0

最终得出了前面的共轭条件,证毕。

了解了这个性质,我们可以对比一下最速下降法和共轭方向法(现在还叫共轭方向),他们的区别:

最速下降法 共轭方向法
方向 梯度 共轭的方向
步长 使方向取得极值的位置 使方向取得极值的位置

这么一比较发现,两者的区别主要在方向上了,其他的原理似乎比较一致。

性质2: 共轭梯度能够保证同一个方向只走一次

# step1:构造误差,我们把误差看作是所有共轭方向方法的一个线性组合,有
e(0) = x(0) - x = sum([delta[j] * d[j] for j in range(size(iterations))])
# 其中delta表示了每个共轭方向的线性组合系数

# step2:对上面的公式左右乘以d[k] * A,有
d[k] * A * e(0) = sum([delta[j] * d[k] * A * d[j] for j in range(size(iterations))])
# 根据共轭的概念,除自身以外所有的共轭方向相互之间应该是共轭的(一句废话),那么有如果i != j, d[i] * A * d[j] = 0,所以
d[k] * A * e(0) = delta[k] * d[k] * A * d[k]
=> delta[k] = (d[k] * A * e(0)) / (d[k] * A * d[k])
# 继续根据共轭的性质,对于i < k, d[i] * A * d[k] = 0,
# 所以分子部分可以加一项: d[k] * A * sum([d[i] for i in range(0,k)])
# 那么两项合并,得到
=> delta[k] = (d[k] * A * e(k)) / (d[k] * A * d[k])
=> delta[k] = -(d[k] * r[k]) / (d[k] * A * d[k])
# 根据前面的结论,得到step[k] = -delta[k],说明了这组线性组合的系数和步长可以一一对应的,说明我们的方向在满足共轭时,和它们相应的步长组合可以构成整个误差,也说明了每一步解决一个方向的正确性

(证毕)

这个性质让我们对共轭方向的方法充满了信心,他确实会很快,现在就剩下一个问题,共轭方向d要怎么确定?

共轭方向

对于共轭方向,我们可以采用一个经典的套路进行构建,那就是Gram-Schmidt方法,这是对于一个满秩矩阵求正交向量的经典方法。

它的思路是这样

def gramSchidt(u): # u是向量的数组
	d = [] # d保存结果
	for i in range(len(u)):
		if i == 0: # 第一向量就是它本身
			d.append(u[i])
		else:
			# 初始化一个向量d_cur 
			for j in range(len(d)):
				d_cur = beta(u[i], d[j]) * d[j]
			d.append(u[i] + d_cur)
	return d

也就是说,每一个要新加入的正交(共轭)方向的基,必须把它在前面基的投影部分减掉,剩下的就是与前面的正交(共轭)基的正交(共轭)部分。这样就保证了性质的正确性,因为前提条件是u这组向量线性无关,所以在正交(共轭)化的过程中不会出现一个向量被变成了0向量的悲剧事情。

从上面的代码中,我们发现了一个叫beta的东西,它是前面共轭基的分量,表示了新的向量在这个基上投影后的大小,我们可以用下面的推导求证:

# 根据上面的代码,可以知道
d[i] = u[i] + sum([beta(u[i],d[j]) * d[k]] for k in range(0,i)) # 后面简写beta(...)为beta(i,j)
# 两边同时乘以 A* d[j],这是正交系证明的惯用手法
=> d[i] * A * d[j] = u[i] * A * d[j] + sum(beta(i,j) * d[k] * A * d[j] for k in range(0,i))
=> 0 = u[i] * A * d[j] + beta(i,j) * d[j] * A * d[j]
=> beta(i,j) = -(u[i] * A * d[j]) / (d[j] * A * d[j])

这就是beta的公式。

至于u怎么来呢?我们知道对于一个n维空间,最多构成n个基的向量空间,也就是说我们可以人为构建一组基(比方说全部平行于坐标轴的一组基),还有利用他们计算共轭方向,不就可以了?

可以是可以,就是有点慢……

上面这个方法需要保存优化中每一步的方向,而且计算还要牵扯其中,比较麻烦。所以,得像个办法才行。

对了上面这个方法叫共轭方向法,和共轭梯度还有些差别,那么是时候把共轭梯度拿出来了……

共轭梯度

共轭梯度法,全程是共轭后的梯度——法,十分清楚地表达了,这里的正交基来自梯度。因为在最速下降法中,梯度拥有一个很好的性质,那就是当前的梯度和前一个方向正交,原因是前面一步选择当前方向的最小值,导致下一步的梯度不会出现该方向的分量,那么这个正交性保证了每一轮出现的梯度都不会与前面的共轭方向线性有关,这是个非常好的性质。

TODO

整体代码

def conjugateGradient(A,b, x0, max_iter):
	d = []
	step = []
	r = []
	x = [x0]
	beta = []
	d[0] = b - A * x[0]
	for i in range(size(max_iter)):
		step = (r(i) * r(i)) / (d[i] * A * d[i])
		x[i+1] = x[i] + step[i] * d[i]
		r[i+1] = r[i] - step[i] * A * d[i]
		beta[i+1] = (r(i+1) * r(i+1))/(r(i)*r(i))
		d[i+1] = r(i+1) + beta[i+1] * d[i]