代码传送门:https://github.com/taifyang/optimization-method
Python实现:
import sympy
import numpy as np
from numpy import matlib as mb
#f为要求极值的函数,x0为初始位置,max_iter为最大迭代次数,epsilon为相邻两次迭代的x改变量
def Broyden_x(f, x0, max_iter, epsilon):
i = 0 # 记录迭代次数的变量
x0 = float(x0) # 浮点数计算更快
df = sympy.diff(f, x) # 定义一阶导数
beta = 0.5 #beta 0~1
delta = 0.25 #delta 0~0.5
Hk = 1 #初始对称正定矩阵H0通常取为1或单位矩阵
while i < max_iter:
gk = df.subs(x, x0)
dk = -Hk*gk
mk = 0
while mk < 10:
if f.subs(x, x0+beta**mk*dk) < f.subs(x,x0) + delta*beta**mk*gk*dk:
break
mk += 1
xnew = x0 + beta**mk*dk
sk = xnew - x0
yk = df.subs(x, xnew) - gk
phik = 0.5
vk = (yk*Hk*yk)**0.5*(sk/(yk*sk)-(Hk*yk)/(yk*Hk*yk))
Hk = Hk - (Hk*yk*yk*Hk)/(yk*Hk*yk) + (sk*sk)/(sk*yk) + phik*vk*vk
i += 1
print('迭代第%d次:%.5f' %(i, xnew))
if abs(f.subs(x, xnew)-f.subs(x, x0)) < epsilon:
break
x0 = xnew
return xnew
#f为要求极值的函数,X0为初始位置,max_iter为最大迭代次数,epsilon为相邻两次迭代的x改变量
def Broyden_x0x1(f, X0, max_iter, epsilon):
i = 0 # 记录迭代次数的变量
X0[0], X0[1] = float(X0[0]), float(X0[1]) #浮点数计算更快
df0 = sympy.diff(f, x0) #定义一阶导数
df1 = sympy.diff(f, x1)
beta = 0.5 #beta 0~1
delta = 0.25 #delta 0~0.5
Hk = mb.identity(len(X0)) #初始对称正定矩阵H0通常取为1或单位矩阵
while i < max_iter:
gk = np.mat([float(df0.subs([(x0, X0[0]), (x1, X0[1])])), float(df1.subs([(x0, X0[0]), (x1, X0[1])]))]).T #梯度矩阵
dk = -Hk*gk
mk = 0
while mk < 10:
if f.subs([(x0, X0[0]+beta**mk*dk[0,0]), (x1, X0[1]+beta**mk*dk[1,0])]) < f.subs([(x0, X0[0]), (x1, X0[1])]) + delta*beta**mk*gk.T*dk:
break
mk += 1
Xnew = [X0[0] + beta**mk*dk[0,0], X0[1] + beta**mk*dk[1,0]]
sk = np.mat([beta**mk*dk[0,0], beta**mk*dk[1,0]]).T
yk = np.mat([float(df0.subs([(x0, Xnew[0]), (x1, Xnew[1])])), float(df1.subs([(x0, Xnew[0]), (x1, Xnew[1])]))]).T - gk
phik = 0.5
vk = float(np.sqrt(yk.T*Hk*yk))*(sk/(sk.T*yk)-(Hk*yk)/(yk.T*Hk*yk))
Hk = Hk - (Hk*yk*yk.T*Hk)/(yk.T*Hk*yk) + (sk*sk.T)/(sk.T*yk) + phik*vk*vk.T
i += 1
print('迭代第%d次:[%.5f, %.5f]' %(i, Xnew[0], Xnew[1]))
if abs(f.subs([(x0, Xnew[0]), (x1, Xnew[1])])-f.subs([(x0, X0[0]), (x1, X0[1])])) < epsilon:
break
X0 = Xnew
return Xnew
if __name__ == '__main__':
x = sympy.symbols("x")
x0 = sympy.symbols("x0")
x1 = sympy.symbols("x1")
result = Broyden_x(x**4-4*x, 10, 50, 1e-5)
print('最佳迭代的位置:%.5f' %result)
result = Broyden_x0x1((x0-1)**2+(x1-1)**4, [10,10], 50, 1e-5)
print('最佳迭代位置:[%.5f, %.5f]' %(result[0], result[1]))
C++实现:
#include <iostream>
#include <vector>
#include <Eigen/Dense>
const double dx = 1e-3;
double f(double x)
{
return pow(x, 4) - 4 * x;
}
double df(double x)
{
//return 4 * pow(x, 3) - 4;
return (f(x + dx) - f(x)) / dx;
}
double f(std::vector<double> X)
{
return pow(X[0] - 1, 2) + pow(X[1] - 1, 4);
}
double df0(std::vector<double> X)
{
//return 2 * (X[0] - 1);
return (f({ X[0] + dx, X[1] }) - f(X)) / dx;
}
double df1(std::vector<double> X)
{
//return 4 * pow(X[1] - 1, 3);
return (f({ X[0], X[1] + dx }) - f(X)) / dx;
}
double Broyden_x(double x0, int max_iter, double epsilon)
{
int i = 0;
double beta = 0.5;
double delta = 0.25;
double Hk = 1;
double xnew;
while (i < max_iter)
{
double gk = df(x0);
double dk = -Hk*gk;
int mk = 0;
while (mk < 10)
{
if (f(x0 + pow(beta, mk)*dk) < f(x0) + delta*pow(beta, mk)*gk*dk)
break;
++mk;
}
xnew = x0 + pow(beta, mk)*dk;
float sk = xnew - x0;
float yk = df(xnew) - gk;
float phik = 0.5;
float vk = sqrt(yk*Hk*yk)*(sk / (yk*sk) - (Hk*yk) / (yk*Hk*yk));
Hk = Hk - (Hk*yk*yk*Hk) / (yk*Hk*yk) + (sk*sk) / (sk*yk) + phik*vk*vk;
++i;
std::cout << "迭代次数:" << i << " " << x0 << std::endl;
if (abs(f(xnew) - f(x0)) < epsilon)
break;
x0 = xnew;
}
return xnew;
}
std::vector<double> Broyden_x0x1(std::vector<double> X0, int max_iter, double epsilon)
{
int i = 0;
double beta = 0.5;
double delta = 0.25;
Eigen::Matrix2d Hk = Eigen::Matrix2d::Identity();
std::vector<double> Xnew;
while (i < max_iter)
{
Eigen::Vector2d gk;
gk << df0(X0), df1(X0);
Eigen::Vector2d dk = -Hk*gk;
int mk = 0;
while (mk < 10)
{
Xnew = { X0[0] + pow(beta, mk)*dk(0), X0[1] + pow(beta, mk)*dk(1) };
if (f(Xnew) < f(X0) + delta*pow(beta, mk)*gk.transpose()*dk)
break;
++mk;
}
Xnew = { X0[0] + pow(beta, mk)*dk(0), X0[1] + pow(beta, mk)*dk(1) };
Eigen::Vector2d sk;
sk << pow(beta, mk)*dk(0), pow(beta, mk)*dk(1);
Eigen::Vector2d yk;
yk << df0(Xnew), df1(Xnew);
yk -= gk;
float phik = 0.5;
Eigen::Vector2d vk = sqrt(yk.transpose()*Hk*yk)*(sk / (sk.transpose()*yk)- (Hk*yk) / (yk.transpose()*Hk*yk));
Hk = Hk - (Hk*yk*yk.transpose()*Hk) / (yk.transpose()*Hk*yk) + (sk*sk.transpose()) / (sk.transpose()*yk) + phik*vk*vk.transpose();
++i;
std::cout << "迭代次数:" << i << " " << X0[0] << " " << X0[1] << std::endl;
if (abs(f(Xnew) - f(X0)) < epsilon)
break;
X0 = Xnew;
}
return X0;
}
int main(int argc, char* argv[])
{
double result = Broyden_x(10, 50000, 1e-5);
std::cout << "最佳迭代位置:" << result << std::endl;
std::vector<double> results = Broyden_x0x1({ 10,10 }, 50000, 1e-5);
std::cout << "最佳迭代位置:" << results[0] << " " << results[1] << std::endl;
system("pause");
return EXIT_SUCCESS;
}