【Matlab学习手记】拟牛顿型信頼域方法求解函数极值
信頼域子问题求解过程,包含Hessian矩阵计算啊,这里采用拟牛顿法(BFGS)来计算近似的Hessian矩阵。
例子比较简单:
f = 100*(x(1)^2 - x(2))^2 + (x(1) - 1)^2;
很容易看出,此函数有唯一一个全局极小值点[1, 1],下面给出拟牛顿信頼域法的计算过程,给一个偏离很远的初始值[-100, 100],经过205次迭代收敛到[1, 1]。
- 源码
function x = Opt_TrustRegionBFGS(x0, err, MaxIter)
%{
函数功能:BFGS信赖域方法求解无约束问题: min f(x);
输入:x0:初始点;err:梯度误差阈值;MaxIter:最大迭代次数;
输出:x:最小值点;
调用格式:
clear; clc;
x = Opt_TrustRegionBFGS([-100 100]', 1e-5, 1000)
%}
eta1 = 0.25;
eta2 = 0.75;
tau1 = 0.5;
tau2 = 2;
delta = 1;
x = x0;
n = length(x0);
Bk = eye(n); % 初始 Hessian 矩阵
k = 0;
while k < MaxIter fk = fun(x);gk = gfun(x);if norm(gk) < err disp(k);break;end% 调用子程序 Trust_q[d, val] = Trust_q(fk, gk, Bk, delta);dq = fk - val;df = fun(x) - fun(x + d);rk = df/dq;if rk <= eta1 delta = tau1*delta;elseif rk >= eta2 && norm(d) == delta delta = tau2*delta;endif rk > eta1 x = x + d;dx = d; % 迭代点差值dg = gfun(x) - gk; % 梯度差值if dg'*dx > 0 Bk = Bk - (Bk*(dx*dx')*Bk)/(dx'*Bk*dx) + (dg*dg')/(dg'*dx); % Hesse矩阵更新endendk = k+1;
endfunction [d, val] = Trust_q(fk, gk, Bk, deltak)
%{函数功能: 求解信赖域子问题:min qk(d)=fk+gk'*d+0.5*d'*Bk*d, s.t.||d|| <= delta
输入: fk:xk处的目标函数值;gk:xk处的梯度;Bk:近似Hesse阵;delta:当前信赖域半径
输出: d:最优值;val:最优值
%}
n = length(gk);
rho = 0.6;
sigma = 0.2;
mu0 = 0.05;
lam0 = 0.05;
gamma = 0.05;
epsilon = 1e-6;
d0 = ones(n, 1);
zbar = [mu0, zeros(1, n + 1)]';
i = 0;
mu = mu0;
lam = lam0;
d = d0;
while i <= 150H = dah(mu, lam, d, gk, Bk, deltak);if norm(H) <= epsilonbreak;endJ = JacobiH(mu, lam, d, Bk, deltak);b = psi(mu, lam, d, gk, Bk, deltak, gamma)*zbar - H;dz = J\b;dmu = dz(1); dlam = dz(2); dd = dz(3 : n + 2);m = 0; mi = 0;while m < 20t1 = rho^m;Hnew = dah(mu + t1*dmu, lam + t1*dlam, d + t1*dd, gk, Bk, deltak);if norm(Hnew) <= (1 - sigma*(1 - gamma*mu0)*rho^m)*norm(H) mi = m; break;endm = m+1;endalpha = rho^mi;mu = mu + alpha*dmu;lam = lam + alpha*dlam;d = d + alpha*dd;i = i + 1;
end
val = fk + gk'*d + 0.5*d'*Bk*d;function p = phi(mu, a, b)
p = a + b - sqrt((a - b)^2 + 4*mu^2);function H = dah(mu, lam, d, gk, Bk, deltak)
n = length(d);
H = zeros(n + 2, 1);
H(1) = mu;
H(2) = phi(mu, lam, deltak^2 - norm(d)^2);
H(3 : n + 2) = (Bk + lam*eye(n))*d + gk;function J = JacobiH(mu, lam, d, Bk, deltak)
n = length(d);
t2 = sqrt((lam + norm(d)^2 - deltak^2)^2 + 4*mu^2);
pmu = -4*mu/t2;
thetak = (lam + norm(d)^2 - deltak^2)/t2;
J=[1, 0, zeros(1, n);pmu, 1 - thetak, -2*(1 + thetak)*d';zeros(n, 1), d, Bk + lam*eye(n)];function si = psi(mu, lam, d, gk, Bk, deltak, gamma)
H = dah(mu, lam, d, gk, Bk, deltak);
si = gamma*norm(H)*min(1, norm(H));function f = fun(x)
f = 100*(x(1)^2 - x(2))^2 + (x(1) - 1)^2;function gf = gfun(x)
gf = [400*x(1)*(x(1)^2 - x(2)) + 2*(x(1) - 1);-200*(x(1)^2 - x(2))];
【Matlab学习手记】拟牛顿型信頼域方法求解函数极值相关推荐
- MATLAB求解函数极值及函数图像
MATLAB具有求解函数极值以及函数图像的功能,简单举一个例子. 求解上述函数极值与图像: 1.驻点求解 syms x >> y = (3*x^2 + 4*x +4)/(x^2 + x + ...
- 爬山算法求解函数极值(matlab实现)
爬山算法: 爬山算法是一种简单的贪心搜索算法,在算法迭代的过程中,会从当前解的临近空间中随机选取下一个点,如果比当前结果好则会选取这个点作为新的最优解,否则再次进行选取.因为不是遍历得到的最优解,而是 ...
- 数值优化(Numerical Optimization)学习系列-拟牛顿方法(Quasi-Newton)
概述 拟牛顿方法类似于最速下降法,在每一步迭代过程中仅仅利用梯度信息,但是通过度量梯度之间的变化,能够产生超线性的收敛效果.本节主要学习一下知识点: 1. 拟牛顿方程推导 2. 几个常见的拟牛顿方法 ...
- Matlab学习手记——非线性方程组求解:牛顿下山法
功能:牛顿下山法求解非线性方程组. 牛顿下山法 function [x, n] = NonLinearEquations_NewtonDown(x0, err) %{ 函数功能:牛顿下山法求解非线性方 ...
- 【Matlab学习手记】Matlab积分问题
一个程序彻底搞懂Matlab的数值积分.符号积分问题. 数值积分问题,给定被积分函数和积分上下限,使用 integral 函数得到积分值: 符号积分问题,通常结果是解析解,即需要知道被积分函数的原函数 ...
- Matlab学习手记——非线性拟合方法:压缩因子粒子群算法
目的:采用压缩因子粒子群算法实现双指数拟合. function x_opt = PSO_ExpFit2(t, Et) %{ 函数功能:压缩因子粒子群算法实现指数拟合:y = a1*exp(-x/b1) ...
- 【Matlab学习手记】标签显示在刻度之间
问题:Matlab标签和刻度线默认是对齐的,如何将标签设置到刻度线之间? 三个实例. plot类型 clear; clc; x = 0:0.1:2*pi; y = sin(x); plot(x, y) ...
- Matlab学习手记——制作GIF动图
目的:利用Matlab制作GIF动图. 结果图 测试代码 clear;clc; filename = '页岩碎屑.gif'; % 保存文件名 Iters = [1:9 10*(1:9) 100*(1: ...
- 【Matlab学习手记】sym8小波滤波
提供sym8小波,四层全局软阈值滤波源代码,采用Matlab语言编写,可移植性强. 源代码 clear;clc; load leleccum; indx = 1:3450; noisez = lele ...
最新文章
- R语言威布尔分布函数F Distribution(dweibull, pweibull, qweibull rweibull )实战
- jquery如何判断div是否隐藏
- Navicat设置unique报错
- 自己动手制作chm格式开源文档
- linux下php启动实例,linux下实现定时执行php脚本_php实例
- 模板方法源码解析(jdk+servlet+mybatis)
- Hive的10种常用优化总结,再也不怕MapReduce分配不均了
- python3 asyncio 不阻塞_Python中的并发处理之asyncio包使用的详解
- Chrome启动后打开第一个网页很慢的解决方案
- 面向对象程序设计之封装性、继承性、多态性
- 查看Linux内核版本命令
- 什么是主数据,如何做好主数据管理?
- EUI插件服务器负载显示不兼容,EUI - 魔兽世界最贴心的插件
- 用云服务器架设好服务器显示无法连接
- 开机广告页面2017流行样式 dialogTheme的popuwindow版本
- 大佬云集的在线少儿英语市场,谁才是那匹冲出重围的黑马?
- html、css 实现一个漂亮的表格
- Wireshark实践总结
- linux中mysql的安装与卸载_linux的mysql安装与卸载
- Linux · 教程
热门文章
- 电脑之间快速传输超大文件(100GB以上)的方法
- 加盐密码哈希:如何正确使用 (转)
- 查看计算机内存过高,物理内存过高怎么办,小编教你电脑物理内存过高怎么办...
- 智慧景区“数字孪生“三维可视化运营管理平台-景区“元宇宙”的数字
- 零基础入门数据挖掘之金融风控-贷款违约预测
- 小米手机5X获得Root权限的方法
- 程序员面试闪充--iOS密码学
- 计算机毕业设计选题、开题、答辩、模板大全(有源码)
- 小猿圈Linux零基础自学之路
- 最新小漫画Android下载,迷妹漫画安卓app2021最新版