最优化第六讲练习题

使用牛顿法

def f(vec):
    x1,x2=vec[0],vec[1]
    return x1*x1/2+2*x2*x2

def first_order(vec):
    x1,x2=vec[0],vec[1]
    return np.array((x1,4*x2))

x0=np.array((2,1)) #初始点
sec=np.array([[1,0],[0,4]]) #二阶导
try:
    inv=np.linalg.inv(sec)
except:
    print("矩阵不存在逆矩阵")
res=first_order(x0) #微分
d0=inv@-res
x1=x0+d0
tmp=None
while abs(f(x1)-f(x0))>1e-3: #函数
    res=first_order(x1) #微分
    d1=inv@-res
    tmp=x1
    x1=x1+d1
    x0=tmp

最小值为0

1.1 坐标轴交替下降
令$x_i=x_{i-1}+{\alpha }_i{e}_i$，其中${\alpha }_i=arg\min f(x_{i-1}+\alpha e_i)$

def f(vec):
    x1,x2=vec[0],vec[1]
    return x1*x1+x2*x2+x1*x2+2*x1-3*x2

def first_order(vec):
    x1,x2=vec[0],vec[1]
    return np.array((2*x1+x2+2,2*x2+x1-3))

def calc_alpha(vec,idx):
    x1,x2=vec[0],vec[1]
    if idx:
        return (3-x1-2*x2)/2
    else:
        return -(2+x2+2*x1)/2
    
x0=np.array((0,0))
tmp=x0.copy()
x0[0]+=calc_alpha(x0,0)
n=1
while abs(f(x0)-f(tmp))>1e-3:
    tmp=x0.copy()
    x0[n]+=calc_alpha(x0,n)
    n=1-n

最小值为-6
1.2 最速下降法
$d_k=-\nabla f(x^k)$

class Poly:
    cons1=None
    var1=None
    cons2=None
    var2=None

def poly_weifen(vec,x):
    vec.cons1*=2
    vec.cons1+=x[1]
    vec.cons1+=2
    vec.var1*=2

    vec.cons2*=2
    vec.cons2+=x[0]
    vec.cons2-=3
    vec.var2*=2
    return vec

def calc_alpha(vec,d):
    vec.cons1*=d[0]
    vec.var1*=d[0]
    vec.cons2*=d[1]
    vec.var2*=d[1]
    return -(vec.cons1+vec.cons2)/(vec.var1+vec.var2)

x0=np.array((0,0)) #初始点
d0=-first_order(x0) #微分
x2=Poly()
x2.cons1=x0[0]
x2.var1=d0[0]
x2.cons2=x0[1]
x2.var2=d0[1]
vec=poly_weifen(x2,x0) #多项式微分
alpha=calc_alpha(vec,d0)
x1=x0+alpha*d0
tmp=None
while abs(f(x1)-f(x0))>0.2: #函数定义
    d1=-first_order(x1) #微分
    x2=Poly()
    x2.cons1=x1[0]
    x2.var1=d1[0]
    x2.cons2=x1[1]
    x2.var2=d1[1]

    vec=poly_weifen(x2,x1) #多项式微分
    alpha=calc_alpha(vec,d1)
    tmp=x1
    x1=x1+alpha*d1
    x0=tmp

最小值为-6.30859375
1.3 牛顿法
$d^k=-[{\nabla }^{2}f(x^k){]}^{-1}\nabla f(x^k)$
最优解$x=x^k+d^k$

def calc_sec(vec):
    return np.array([[2,1],[1,2]])

最小值为-6.333333333333332

2.1 坐标轴交替下降

def f(vec):
    x1,x2=vec[0],vec[1]
    return 4*x1*x1+x2*x2-x1*x1*x2

def first_order(vec):
    x1,x2=vec[0],vec[1]
    return np.array((8*x1-2*x1*x2,2*x2-x1*x1))

def calc_alpha(vec,idx):
    x1,x2=vec[0],vec[1]
    if idx:
        return (x1*x1-2*x2)/2
    else:
        return (2*x2-8)*x1/(8-2*x2)

3个初始点的情况下最小值都为0
2.2 最速下降法

def poly_weifen(vec,x):
    vec.cons1*=8
    vec.cons1-=2*x[0]*x[1]
    vec.var1*=8
    vec.var1-=2*x[1]

    vec.cons2*=2
    vec.cons2-=x[0]*x[0]
    vec.var2*=2
    return vec

初始点(1,1)的最小值为3.1338726246602966e-05，初始点(3,4)无法收敛，初始点(2,0)的最小值为3.5938810907227845e-06
1.3 牛顿法

def calc_sec(vec):
    x1,x2=vec[0],vec[1]
    return np.array([[8-2*x2,-2*x1],[-2*x1,2]])

初始点(1,1)的最小值为0.0014071388348340429；初始点(3,4)无法收敛，最小值为16；初始点(2,0)无法收敛