%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright (C) 2020 Zaiwen Wen, Haoyang Liu, Jiang Hu
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see .
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Change log:
%
% 2020.4.18 (Jiang Hu):
% First version
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 利用牛顿-共轭梯度法解优化问题
% 该算法利用非精确牛顿法（牛顿-共轭梯度法）求解无约束优化问题：
%
% $$ \min_{x} \;\; f(x)$$
%
% 在第 $k$ 步迭代，下降方向 $d^k$ 通过求解下面的牛顿方程
% $\displaystyle (\nabla^2 f(x^k)) d^k = -\nabla f(x^k)$
% 得到。选取合适的步长 $\alpha_k$，牛顿法的迭代格式为
%
% $$ x^{k+1} = x^{k} + \alpha_k d_k.$$
%
% 对于规模较大的问题，精确求解牛顿方程组的代价比较高。事实上，牛顿方程求解等价于无约束二次优化问题：
%
% $$ \displaystyle\min_{d^k} \frac{1}{2}(d^k)^\top \nabla^2 f(x^k) d^k +
% (\nabla f(x^k))^\top d^k, $$
%
% 其可以通过共轭梯度法来进行求解。
%
%% 初始化和迭代准备
% 输入信息：迭代初始点 |x| ，目标函数 |fun| （依次返回给定 |x| 处的函数值和梯度），
% 海瑟矩阵 |hess| 以及提供参数的结构体 |opts| 。
%
% 输出信息：迭代得到的解 |x| 和迭代信息结构体 |out| 。
%
% * |out.msg| ：标记是否达到收敛
% * |out.nrmG| ：迭代退出时的梯度范数
% * |out.iter| ：迭代退出时的迭代步数、函数值
% * |out.f| ：迭代退出时的目标函数值
% * |out.nfe| ：调用原函数的次数
function [x, out] = fminNewton(x, fun, hess, opts, varargin)
%%%
% 检查输入信息，MATLAB 以 |nargin| 表示函数输入的参数个数。当参数量不足 3 时报错，
% 等于 3 时认为 |opts| 为空结构体，即全部采用默认参数。
%
% 从输入的结构体 |opts| 中读取参数或采取默认参数。
%
% * |opts.xtol| ：针对优化变量的停机准则
% * |opts.gtol| ：针对梯度范数的停机准则
% * |opts.ftol| ：针对函数值的停机准则
% * |opts.rho1| ：线搜索参数 $\rho_1$
% * |opts.rho2| ：线搜索参数 $\rho_2$
% * |opts.maxit| ：主循环的最大迭代次数
% * |opts.verbose| ： |>=1| 时输出迭代信息
% * |opts.itPrint| ：每隔几步输出一次迭代信息
if (nargin < 3); error('[x, out] = fminNewton(x, fun, hess, opts)'); end
if (nargin < 4); opts = []; end
if ~isfield(opts,'gtol'); opts.gtol = 1e-6; end
if ~isfield(opts,'xtol'); opts.xtol = 1e-6; end
if ~isfield(opts,'ftol'); opts.ftol = 1e-12; end
if ~isfield(opts, 'rho1'); opts.rho1 = 1e-4; end
if ~isfield(opts, 'rho2'); opts.rho2 = 0.9; end
if ~isfield(opts,'maxit'); opts.maxit = 200; end
if ~isfield(opts,'verbose'); opts.verbose = 0; end
if ~isfield(opts,'itPrint'); opts.itPrint = 1; end
%%%
% 从结构体 |opts| 中复制参数。
maxit = opts.maxit; verbose = opts.verbose; itPrint = opts.itPrint;
xtol = opts.xtol; gtol = opts.gtol; ftol = opts.ftol;
%%%
% |parsls| 为由线搜索参数构成的结构体。
parsls.ftol = opts.rho1; parsls.gtol = opts.rho2;
%%%
% 迭代准备，计算初始点 $x$ 处的函数值和梯度值。
[f,g] = fun(x);
nrmg = norm(g,2);
nrmx = norm(x,2);
%%%
% 预先设定输出信息，在没有达到收敛条件时，记为超出最大迭代次数退出。
out = struct();
out.msg = 'MaxIter';
out.nfe = 1;
out.nrmG = nrmg;
%%%
% 当需要详细输出时，设定输出格式。
if verbose >= 1
if ispc; str1 = ' %10s'; str2 = ' %8s';
else str1 = ' %10s'; str2 = ' %8s'; end
stra = ['%5s', str1, str2, str2, str2, '\n'];
str_head = sprintf(stra, ...
'iter', 'F', 'nrmg', 'fdiff', 'xdiff');
fprintf('%s', str_head);
str_num = ['%4d %+8.7e %+2.1e %+2.1e %+2.1e\n'];
end
%%%
% 非精确牛顿法使用共轭梯度法求解牛顿方程， |opts_tCG| 为共轭梯度法提供参数。
opts_tCG = struct();
%% 迭代主循环
for iter = 1:maxit
%%%
% 记录上一步迭代的信息（上一步的优化变量和对应的目标函数值、梯度、范数）。
fp = f;
gp = g;
xp = x;
nrmxp = nrmx;
%%%
% 调用截断共轭梯度法 |tCG| （见参考页面）不精确地求解牛顿方程得到牛顿方向 $d$。
% 这里函数 |hess(x,d)| 对应矩阵向量乘 $\nabla^2 f(x) d$。
opts_tCG.kappa = 0.1;
opts_tCG.theta = 1;
hess_tCG = @(d) hess(x,d);
[d, ~, ~] = tCG(x, gp, hess_tCG, inf, opts_tCG);
%%%
% 沿方向 $d$ 做线搜索。调用函数 |ls_csrch| 进行线搜索，具体步骤可以参考
% MINPACK-2。
%
% 首先初始化线搜索标记 |workls.task| 为 1， |deriv| 为当前线搜索方向上的导数。
% 线搜索函数寻找合适的步长 $\alpha_k$，在对应的搜索方向上满足条件
%
% $$ \begin{array}{rl}
% f(x^k+\alpha_k d^k)&\hspace{-0.5em}\le
% f(x^k)+\rho_1\cdot\alpha_kg(x^k), \\
% |g(x^k+\alpha_kd^k)|&\hspace{-0.5em}\le \rho_2\cdot |g(x^k)|.
% \end{array} $$
%
% |ls_csrch| 每次调用只执行线搜索的一步，并用 |workls.task| 指示下一步应当执行的操作。
% 此处 |workls.task==2| 意味着需要重新计算当前点函数值和梯度等。具体的步长寻找过程可以参考源文件。
%
% 当线搜索条件满足时，退出线搜索循环，得到更新的迭代点 $x$。
workls.task = 1;
deriv = d'*g;
normd = norm(d);
stp = 1;
while 1
[stp, f, deriv, parsls, workls] = ....
ls_csrch(stp, f, deriv , parsls , workls);
if (workls.task == 2)
x = xp + stp*d;
[f, g] = feval(fun, x, varargin{:});
out.nfe = out.nfe + 1;
deriv = g'*d;
else
break
end
end
%%%
% 计算一步更新后的点 $x$处的相关数据以判断是否停机。
% |nrms| 表示 $\|x^{k+1}-x^k\|_2$，
% |xdiff| 表示 $\|x^{k+1}-x^k\|/\max\{1,\|x^k\|\}$ 为迭代点 $x$的相对变化量，
% |nrmg| 表示 $x^k$ 处梯度范数， |nrmx| 为点 $x$ 的范数，
% |fdiff| 为函数值相对变化量。 |out.nfe| 记录函数值计算的总次数。
nrms = stp*normd;
xdiff = nrms/max(nrmxp,1);
nrmg = norm(g,2);
out.nrmG = [out.nrmG; nrmg];
nrmx = norm(x,2);
out.nfe = out.nfe + 1;
fdiff = abs(fp-f)/(abs(fp)+1);
%%%
% 停机判断：当梯度范数小于阈值，或者函数值的相对变化量和优化变量的相对变化量均小于阈值时，
% 迭代终止。
cstop = nrmg <= gtol || (abs(fdiff) <= ftol && abs(xdiff) <= xtol);
%%%
% 当需要详细输出时，在（1）开始迭代时（2）达到收敛时（3）达到最大迭代次数退出迭代时
% （4）每若干步，打印详细结果。
% 当满足收敛条件时，标记为达到最优值。退出迭代，记录退出信息为达到最优值退出。
if verbose>=1 && (cstop || iter == 1 || iter==maxit || mod(iter,itPrint)==0)
if mod(iter,20*itPrint) == 0 && iter ~= maxit && ~cstop
fprintf('\n%s', str_head);
end
fprintf(str_num, ...
iter, f, nrmg, fdiff, xdiff);
end
if cstop
out.msg = 'Optimal';
break;
end
end
%%%
% 迭代退出时，记录输出信息。
out.iter = iter;
out.f = f;
out.nrmg = nrmg;
end
%% 参考页面
% 算法中的共轭梯度法部分，详见 。在 中，我们展示该算法的一个应用。
%
% 代码中利用了翻译为 MATLAB 代码的 MINPACK-2 的线搜索函数 |ls_csrch| ，详见
%
% ，或参考 Fortran 版本的官方代码 。
%
% 此页面的源代码请见： <../download_code/newton/fminNewton.m fminNewton.m>。
%% 版权声明
% 此页面为《最优化：建模、算法与理论》、《最优化计算方法》配套代码。
% 代码作者：文再文、刘浩洋、户将，代码整理与页面制作：杨昊桐。
%
% 著作权所有 (C) 2020 文再文、刘浩洋、户将