跳转至

优化方法笔记

凸优化概括

变量选择之LASSO—(一)凸正则化方法1

least-square

参考:

https://textbooks.math.gatech.edu/ila/least-squares.html

https://zhuanlan.zhihu.com/p/31341436

https://en.wikipedia.org/wiki/Ordinary_least_squares

对于方程

Ax=b
Ax=b

最小二乘法的解为

ˆx=(ATA)1ATb
^x=(ATA)1ATb

不适用于不满秩的情况。

Tikhonov regularization/Ridge Regression(RR)

参考:

https://en.wikipedia.org/wiki/Ridge_regression

也称岭回归,岭回归是在线性回归模型中存在多共线性(高度相关)自变量时,作为解决最小二乘估计量不精确问题的一种可能方法.

通过在对角线上添加正元素来缓解,从而减少其条件数。类似于普通的最小二乘估计器,λλ为正则化阻尼Regularization dampings

ˆβR=(XTX+λI)1XTy

在最小二乘最小化过程中加入正则化项,防止过拟合,获得更加理想的解,类似于L2正则化。

Axb22+Γx22

解为

ˆx=(AA+ΓΓ)1Ab.

Tikhonov matrix Γ=αI

Generalized Tikhonov regularization

x=(APA+Q)1(APb+Qx0)

其中P为b的协方差矩阵,x0x的期望值

Lavrentyev regularization

如果A是对称正定的,即有A=A>0, 所以x2P=xA1x ????

故最优解为

x=(A+Q)1(b+Qx0)

可以看作Tikhonov regularization中A=A=P1的情况

Bregman method

参考:

如何理解Bregman divergence?

Bregman 散度(Bregman divergence)和Bregman信息(Bregman information)

机器学习中的数学——距离定义(二十五):布雷格曼散度(Bregman Divergence)

Split-Bregman迭代方式

定义Bregman散度为

DF(p,q)=F(p)F(q)F(q),pq

几何意义如下图

选择不同的F会有不同的损失函数

例如最小均方误差

d2(x,y)=xy2=xy,xy=x2y22y,xy

Bregman Iteration

uk+1=\argminuDpkJ(u,uk)+λHH(u),pkJ(uk)

等价于

uk+1=argminuJ(u)p,uuk+λH(u),pkJ(uk)

Result,证明见文档《Notes on Bregman Iteration》

HH(uk+1)HH(uk+1)+1λDpkJ(uk+1,uk)HH(uk)
{uk+1=\argminuDpkJ(u,uk)+λH(u)pk+1=pkλH(uk+1)

假设H(u,f)=12||Auf||22,则p0=0时的迭代过程等价于 $$ \left{ \begin{array}{l} u_{k+1}=\arg \min {u} J(u)+\lambda H\left(u, f{k}\right) \ f_{k+1}=f_{k}+\left(f-A u_{k+1}\right), \quad f_{0}=f \end{array}\right. $$

Linearized Bregman

假设J(u)=||u||1,将H(u)进行矩阵泰勒展开得

H(u)H(uk)+H(uk)(uku)

但是这种近似展开只有在uuk十分接近的时候上式才准确,所以加入惩罚项12δuuk22,最终如下

uk+1=argminuDpkJ(u,uk)+λH(uk)+λH(uk),uuk+12δuuk22

化简得 $$ u_{k+1}=\arg \min {u} J(u)+\left\langle\lambda \nabla H\left(u{k}\right)-p_{k}, u\right\rangle+\frac{1}{2 \delta}\left|u-u_{k}\right|_{2}^{2} $$

Split Bregman

假设ΦE(u)均为凸函数

minuΦ(u)1+E(u)

通过运算符分裂

minu,dd1+E(u) subject to Φ(u)=d

J(u,d)=||d||1+E(u),H(u,d)=12||dΦ(u)||22,带入Bregman方法

{(uk+1,dk+1)=minu,dJ(u,d)pku,uukpkd,ddk+λH(u,d)pk+1u=pkuλuH(uk+1,dk+1)pk+1d=pkdλdH(uk+1,dk+1)

lsqr方法

最小二乘法原理详解

最小平方QR分解法

min[AλI]x[b0]2

基础方法即代码实现

对于问题

min[AλI]x[b0]2

在之后的代码实现中,令G=[AλI],obs=[b0]

其中λI为正则化项,用于为反问题增加约束,减小多解性

梯度下降

通过梯度下降最小化目标泛函f(m)的迭代流程为

mk+1=mkαkf(mk)

对于Ax=b的问题,可以求出

minf(m)=||Gmd||22
f(m)=(Gmd)T(Gmd)=mTGTGm2(GTd)Tm+dTd
f(m)=2GTGm2GTd
mk+1=mkαk(2GTGmk2GTd)

步长的选取

αk2λmax(GTG)

代码实现

def GD_run(self, init_guess, G, obs, Iter, show=False):
    pre = init_guess
    loss = []
    for k in range(Iter):
        step = 1
        grad = G.T @ G @ pre - G.T @ obs
        residual = self.calculate_residual(pre, G, obs)
        rmse = np.linalg.norm(residual)
        # 原定损失小于沿学习率迭代后的损失,说明步长较大,减小步长
        while rmse < np.linalg.norm(self.calculate_residual(pre - step * grad, G, obs)):
            step = step * 0.5
        pre = pre - step * grad
        loss.append(rmse)

    if show:
        plt.figure()
        plt.plot(loss, linewidth=1, linestyle="solid", label="loss")
        plt.legend()
        plt.title('GD Loss curve(final loss = ' + str(round(loss[-1], 4)) + ')')
        plt.xlabel('iters')
        plt.ylabel('res')

    return pre, loss

随机梯度下降

类似梯度下降,只不过不是求全部列的梯度,而是求其随机一列的梯度进行下降

minf(m)=||Gmd||22=mi=1fi(m)=mi=1(Gmd)2if(m)i=2GTiGim2diGTi

代码实现(固定梯度) TODO 变梯度

def SGD_run(self, init_guess, G, obs, iter, show=False):
    self.init_guess = init_guess
    pre = self.init_guess
    rmse_prev = np.inf
    loss = []
    for k in range(iter):
        step = 1  # 步长 
        lossres = self.calculate_residual(pre, self.A, self.b)  # 损失,如果未加正则化,self.A,self.b即为G,obs
        index = round(np.random.random() * (G.shape[0] - 1))  # 随机取行号
        grad = G[index:index + 1, :].T @ G[index:index + 1, :] @ pre - obs[index] * G[index, :].T  # 随机一行的梯度
        pre = pre - step * grad  # 梯度下降
        rmse = np.linalg.norm(lossres)
        loss.append(rmse)

    if show:
        plt.figure()
        plt.plot(loss, linewidth=1, linestyle="solid", label="loss")
        plt.legend()
        plt.title('SGD Loss curve(final loss = ' + str(round(loss[-1], 4)) + ')')
        plt.xlabel('iters')
        plt.ylabel('res')

    return pre, loss

高斯牛顿法

参考链接:

Algorithms from scratch: Gauss-Newton

最小二乘问题的四种解法——牛顿法,梯度下降法,高斯牛顿法和列文伯格-马夸特法的区别和联系

对于最小二乘问题

minx12||f(x)||22

f(x+Δx)进行泰勒展开得到

f(x+Δx)f(x)+J(x)TΔx+o(Δx)

所以最小二乘问题化为

Δx=argmin12||f(x+Δx)||22argmin12||f(x)+J(x)TΔx||22

化简

m(x)=12||f(x)+J(x)TΔx||2=12(f(x)+J(x)TΔx)T(f(x)+J(x)TΔx)

求导令导数等于零可得

m(x)=0J(x)J(x)TΔx=J(x)f(x)Δx=J(x)f(x)J(x)J(x)T

代码实现

def GN_run(self, init_guess, G, obs, iter, show=False):
    pre = init_guess
    loss = []
    for k in range(iter):
        residual = self.calculate_residual(pre, G, obs)
        lossres = self.calculate_residual(pre, self.A, self.b)
        jacobian = self.calculate_jacobian(pre, G, obs, step=10 ** (-6))
        pre = pre - np.linalg.pinv(jacobian.T @ jacobian) @ (jacobian.T @ residual)
        rmse = np.linalg.norm(lossres)
        loss.append(rmse)
    if show:
        plt.figure()
        plt.plot(loss, linewidth=1, linestyle="solid", label="loss")
        plt.legend()
        plt.title('Gauss-Newton Loss curve(final loss = ' + str(round(loss[-1], 4)) + ')')
        plt.xlabel('iters')
        plt.ylabel('res')

    return pre, loss

列文伯格-马夸特方法(Levenberg–Marquardt)

基本原理与高斯牛顿法类似,在高斯牛顿法中,可以看到,损失曲线在后期会进入锯齿状,迭代的次数较长。所以,为了避免其步长过大导致的问题,该方法提出了信赖区域,动态调整步长减小迭代次数

在迭代过程中设定一个评判标准用于评判估计的好坏,从而动态调整步长

ρ=f(x+Δx)f(x)J(x)TΔx

分为以下三种情况

  • ρ 接近1,近似是好的,不需要更改
  • ρ 太小,则实际减少的值小于近似减少的值,近似较大,需要缩小近似的范围
  • ρ 太大,则实际减少的值大于近似减少的值,近似较小,需要扩大近似的范围。

反问题为

min12||f(x)+J(x)TΔx||s.t||DΔx<μ||2

其中D为系数矩阵,μ为信赖半径

可以构建拉格朗日函数,λ为系数因子

L(Δx,λ)=12||f(x)+J(x)TΔx||2+λ2(||DΔx||2μ)

求导并化简得

J(x)f(x)+J(x)JT(x)Δx+λDTDΔx=0(JJT+λDTD)Δx=Jf

在实际使用中,通常将I代替DTD

代码实现

def LM_run(self, init_guess, G, obs, iter, show=False):
    pre = init_guess
    lam = 0.01  # todo 如何设定lambda
    loss = []
    for k in range(iter):
        residual = self.calculate_residual(pre, G, obs)
        lossres = self.calculate_residual(pre, self.A, self.b)
        jacobian = self.calculate_jacobian(pre, G, obs, step=10 ** (-6))
        pre = pre - np.linalg.pinv(jacobian.T @ jacobian + lam * np.identity(jacobian.shape[1])) @ (
                jacobian.T @ residual)
        rmse = np.linalg.norm(lossres)
        loss.append(rmse)
    if show:
        plt.figure()
        plt.plot(loss, linewidth=1, linestyle="solid", label="loss")
        plt.legend()
        plt.title('Levenberg-Marquardt Loss curve(final loss = ' + str(round(loss[-1], 4)) + ')')
        plt.xlabel('iters')
        plt.ylabel('res')

    return pre, loss

LM方法是由牛顿法的基础上得到的,他们之间的联系主要取决于λ的取值

对于迭代过程

(JJT+λI)Δxk=Jf
  • λ较大时,更新方程中λI占主要,迭代近似为λIΔxk=Jf,此时为梯度下降法
  • λ较大时,更新方程中JJT占主要,迭代近似为JJTΔxk=Jf,此时为高斯牛顿法

IRLS

正则化反演

参考链接

变量选择之LASSO—(一)凸正则化方法1

对于最小二乘问题min||Axb||2进行反演时,因为约束条件不够,可能会出现多个不同的x满足min||Axb||2,当并不是我们预期的x,也就是问题的多解性,这时我们需要在求解反问题的方程中添加其他约束项,从数学上来讲,就是在损失函数中加个正则项,来防止过拟合。目前主要有L1正则化和L2正则化,即LASSO回归和岭回归

参考链接:l1正则与l2正则的特点是什么,各有什么优势?

  • L1-Regularization(一范数) λki=0|wi|

  • L2-Regularization(二范数) <spanclass="arithmatex"><spanclass="MathJaxPreview">λki=0w2i</span><scripttype="math/tex">λki=0w2i

对于地震反演,常用的就是TV正则化,即对模型的差分的一范数,提高模型的不连续性,可以形成稀疏解

min||Gmd||22+α||Lm||1

其中L为差分矩阵

这里介绍一种求解带有一范数正则化反问题的求解方法 IRLS

对于上式求梯度得到

f(m)=2GTGm2GTd+αmi=0|yi|

由于yi不为零,故

|yi|mk=Li,ksgn(yi)

因为|Lm|=mi=1|yi|,故

|Lm|mk=mi=1Li,kyi|yi|

所以我们令W表示一个对角阵,且

Wi,i=1|yi|

|Lm|=LTWLm

梯度变为

f(m)=2GTGm2GTd+αLTWLm

因为W依赖于Lm,所以这是一个非线性方程组,且W在任意Lm取值为0的点无定义

为解决不可微问题,我们设定一个容许偏差ϵ,并将W改写如下

Wi,i={1/|y(k)i||y(k)i|ϵ 1/ϵ|y(k)i|<ϵ.

迭代过程为

m(k+1)=argmin[Gα2WL]m[d0]2

代码实现,示例为高斯牛顿法优化策略

def GN_run(self, init_guess, G, obs, iter, show=False, IRLS=None):  # IRLS为一范数参数[alpha, eps, L]
    pre = init_guess
    loss = []
    with tqdm(total=iter) as t:
        for i in range(iter):
            t.set_description('solver GN')
            if IRLS is not None:
                L = IRLS[-1]
                alpha = IRLS[0]
                eps = IRLS[1]
                Reg = self._IRLSupdate(L,pre,alpha,eps)
                G = np.vstack((self.A, Reg))
                obs = np.hstack((self.b, np.zeros((Reg.shape[0]))))  # TODO 一维数据适用,二维数据要改!!!
            residual = self.calculate_residual(pre, G, obs)
            lossres = self.calculate_residual(pre, self.A, self.b)
            jacobian = self.calculate_jacobian(pre, G, obs, step=10 ** (-6))
            pre = pre - np.linalg.pinv(jacobian.T @ jacobian) @ (jacobian.T @ residual)
            rmse = np.linalg.norm(lossres)
            loss.append(rmse)
            t.set_postfix(loss=rmse)
            t.update(1)
    if show:
        plt.figure()
        plt.plot(loss, linewidth=1, linestyle="solid", label="loss")
        plt.legend()
        plt.title('Gauss-Newton Loss curve(final loss = ' + str(round(loss[-1], 4)) + ')')
        plt.xlabel('iters')
        plt.ylabel('res')

    return pre, loss

ISTA方法

参考链接

软阈值迭代算法(ISTA算法)

对于反问题

min||Axb||22+λ||x||1

通过推导得到(具体参照上面的参考链接)

xk+1=Tλt(xk2tAT(Axkb))

其中,软阈值函数为:

参考链接:硬阈值 & 软阈值

Tω(x)i=(|xi|ω)+sgn(xi)

[Tω(x)]i={xiω,xi>ω0,|xi|ωxi+ω,xi<ω

当A为单位矩阵时,即

min||xb||22+λ||x||1

解为

x=Tλ2(b)

代码实现

def Ista_run(self, init_guess, G, obs, iter, show=False):
    pre = init_guess
    lamb = 0.01
    tk = 1 / np.linalg.norm(G.T @ G)  # tk最大值

    loss = []
    for k in range(iter):
        residual = self.calculate_residual(pre, G, obs)
        rmse = np.linalg.norm(residual)
        loss.append(rmse)
        grad = G.T @ G @ pre - G.T @ obs
        grad = np.squeeze(grad)
        zk = pre - tk * grad
        pre = self._softThresh(lamb * tk, zk)
        tk = tk / 1.1
    if show:
        plt.figure()
        plt.plot(loss, linewidth=1, linestyle="solid", label="loss")
        plt.legend()
        plt.title('ISTA Loss curve')
        plt.xlabel('iters')
        plt.ylabel('res')
    return pre, loss

FISTA方法

Amir Beck等人将Nesterov加速算法引入ISTA算法中,并称之为FISTA算法,其将复杂度从O(1/k)降低到了O(1/k2),减小了ISTA的计算复杂度,提高了计算速度,使损失可以得到更快收敛

ISTA方法是直接迭代xk,FISTA是迭代ykykxk计算得到

yk+1=xk+(tk1tk+1)(xkxk1)

其中

tk+1=1+1+4t2k2
def Fista_run(self, init_guess, G, obs, iter, show=False):
    pre = init_guess
    tk = 1
    t = 1 / np.linalg.norm(G.T @ G)
    lamb = 0.01
    loss = []
    xk = init_guess
    xk_1 = xk
    for k in range(iter):
        tk1 = 0.5 * (1 + np.sqrt(1 + 4 * tk ** 2))
        yk1 = xk + (tk - 1) / tk1 * (xk - xk_1)
        grad = G.T @ G @ yk1 - G.T @ obs
        grad = np.squeeze(grad)
        zk = yk1 - t * grad

        # update iterative
        xk_1 = xk
        xk = self._softThresh(lamb * t, zk)

        tk = tk1
        # err = np.linalg.norm(self.pre - yk1)
        # if k!=0 and err<0.0001:
        #     print('第%d轮已经收敛,增量为%f' %(k,err))
        #     break
        pre = yk1
        t = t / 1.1

        # loss
        residual = self.calculate_residual(pre, G, obs)
        rmse = np.linalg.norm(residual)
        loss.append(rmse)

    if show:
        plt.figure()
        plt.plot(loss, linewidth=1, linestyle="solid", label="loss")
        plt.legend()
        plt.title('FISTA Loss curve')
        plt.xlabel('iters')
        plt.ylabel('res')
    return pre, loss

Split-Bregman方法

参考见《反问题的迭代求解算法》课件

对于一阶TV反问题

minm||dataGm||22+α||Lm||1

dk=Lm,增加约束项,得到

(mk+1,dk+1)=minm,d||dataGm||22+α||dk||1+λ2dkLmbk22bk+1=bk+(dataGmdk+1)

故,上式可拆解成两个子问题,岭回归问题与LASSO回归问题

子问题一(岭回归)

mk+1=minm||dataGm||22+λ2||dkbkLm||22

可看作

min||[Gλ2L]m[dataλ2(db)]||2

子问题二(LASSO回归)

dk+1=mindλ2||dkbkLm||22+α||dk||1

用软阈值函数可求解

d={b+Lmλ4,b+Lm>λ40,b+Lmλ4b+Lm+λ4,b+Lm<λ4
d={b+Lm2αλ,b+Lm>2αλ0,|b+Lm|2αλb+Lm+2αλ,b+Lm<2αλ

代码实现

def run(self, show=True):
    iiter = 0
    flag = 1
    loss = []
    lamb = 0.01
    x = self.x0
    for _ in tqdm(range(self.Iter)):
        # regularized problem
        dataregs = self.d - self.b_
        # subquestion1
        Sub1 = RegularizedInv(self.Op, self.obs, self.x0, lamb, self.Reg, dataregs)  # 初始化岭回归反演
        x, sub1loss = Sub1.run(self.x0 if flag else x, 10, 'GN', True)  # 选择优化方法迭代求解
        flag = 0
        # subquestion2
        dataregs2 = self.Reg @ x + self.b_
        # self.d, sub2loss = self.Ista_run(self.d,np.eye(self.Reg.shape[0]),dataregs2,50,None) # ISTA方法
        self.d = self._softThresh(lamb / 4, dataregs2)  # TODO lambda 应该怎么确定

        residual = np.linalg.norm(self.calculate_residual(x, self.A, self.b))
        loss.append(residual)

        self.b_ = self.b_ + self.tau * (self.Reg @ x - self.d)
        iiter = iiter + 1

    if show:
        plt.figure()
        plt.plot(loss, linewidth=1, linestyle="solid", label="loss")
        plt.legend()
        plt.title('Split Bregman Loss curve(final loss = ' + str(round(loss[-1], 4)) + ')')
        plt.xlabel('iters')
        plt.ylabel('res')

    return x, loss

对于二阶TV反问题

minmi(xm)2i+(ym)2i+μ2||Gmd||22

dx=xm,dy=ym

进行分裂L1,L2项后,反问题化为

minu,dk,dy||(dx,dy)||2+μ2||Gmd||22+λ2||dxxmbx||22+λ2||dyymby||22

其中

||(dx,dy)||2=i,j|dx,i,j|2+|dy,i,j|2

两个子问题分别为

子问题一

mk+1=minmμ2||Gmd||22+λ2||dxxmbx||22+λ2||dyymby||22

化为岭回归问题

min||[Gλ2xλ2y]m[dλ2(dxbx)λ2(dyby)]||2

[Gλ2xI]m=[dλ2(dxbx)λ2(dyby)1y]

子问题二

(dk+1x,dk+1y)=mindx,dy||(dx,dy)||2+λ2||dxxmbx||22+λ2||dyymby||22

可以通过广义阈值函数求解,参考论文A Fast Algorithm for Image Deblurring with Total Variation Regularization

迭代过程为

dk+1x=max(sk1λ,0)xmk+bkxskdk+1y=max(sk1λ,0)ymk+bkysksk=|xmk+bkx|2+|ymk+bky|2

拟牛顿法

深入机器学习系列17-BFGS & L-BFGS

数值优化(六)——拟牛顿法之LBFGS理论与实战

机器学习基础·L-BFGS算法

数值优化(6)——拟牛顿法:SR1,BFGS,DFP,DM条件

BFGS算法

L-BFGS算法

蒙特卡洛

目标泛函构建

贝叶斯理论

正则化项构建

CV交叉验证

广义交叉验证

多道反演

以上的优化方法都只针对于单道,但地震数据往往需要多道进行处理

目前本人实现的方法只是循环进行不同道的单道反演,未完待续。。。