使用牛顿法
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
+
α
i
e
i
x_i=x_{i-1}+\alpha_ie_i
xi=xi−1+αiei,其中
α
i
=
a
r
g
min
f
(
x
i
−
1
+
α
e
i
)
\alpha_i=arg\min f(x_{i-1}+\alpha e_i)
αi=argminf(xi−1+αei)
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
=
−
∇
f
(
x
k
)
d_k=-\nabla f(x^k)
dk=−∇f(xk)
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
=
−
[
∇
2
f
(
x
k
)
]
−
1
∇
f
(
x
k
)
d^k=-[\nabla^2f(x^k)]^{-1}\nabla f(x^k)
dk=−[∇2f(xk)]−1∇f(xk)
最优解
x
=
x
k
+
d
k
x=x^k+d^k
x=xk+dk
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)无法收敛