高斯-牛顿法

Gauss-Newton Method
考虑
$$\underset{\mathbf{x}\in {\mathbb{R}}^n}{\min }\left\{g(\mathbf{x})\equiv \sum _{i=1}^m{\left(f_i(\mathbf{x})-c_i\right)}^2\right\}$$
假设$f_1,\cdots ,f_m$在${\mathbb{R}}^n$上连续可微
$c_1,\cdots ,c_m\in \mathbb{R}$
有时候会向量化
$$F\left(\mathbf{x}\right)=\left(\begin{array}{c}{f}_1\left(\mathbf{x}\right)-c_1\\ f_2\left(\mathbf{x}\right)-c_2\\ ⋮\\ f_m\left(\mathbf{x}\right)-c_m\end{array}\right)$$
于是问题就是
$$\min \Vert F\left(\mathbf{x}\right){\Vert }^2$$
每一次迭代
$${\mathbf{x}}_{k+1}=\arg \underset{\mathbf{x}\in {\mathbb{R}}^n}{\min }\left\{\sum _{i=1}^m{\left[f_i\left({\mathbf{x}}_k\right)+\nabla f_i{\left({\mathbf{x}}_k\right)}^T\left(\mathbf{x}-{\mathbf{x}}_k\right)-c_i\right]}^2\right\}$$
这个问题可以转化为
$$\underset{\mathbf{x}\in {\mathbb{R}}^n}{\min }\Vert {\mathbf{A}}_k\mathbf{x}-{\mathbf{b}}_k{\Vert }^2$$
其中
$${\mathbf{A}}_k=\left(\begin{array}{c}\nabla f_1{\left({\mathbf{x}}_k\right)}^T\\ \nabla f_2{\left({\mathbf{x}}_k\right)}^T\\ ⋮\\ \nabla f_m{\left({\mathbf{x}}_k\right)}^T\end{array}\right)=J\left({\mathbf{x}}_k\right)$$
也叫做雅可比矩阵
$${\mathbf{b}}_k=\left(\begin{array}{c}\nabla f_1{\left({\mathbf{x}}_k\right)}^T{\mathbf{x}}_k-f_1\left({\mathbf{x}}_k\right)+c_1\\ \nabla f_2{\left({\mathbf{x}}_k\right)}^T{\mathbf{x}}_k-f_2\left({\mathbf{x}}_k\right)+c_2\\ ⋮\\ \nabla f_m{\left({\mathbf{x}}_k\right)}^T{\mathbf{x}}_k-f_m\left({\mathbf{x}}_k\right)+c_m\end{array}\right)=J\left({\mathbf{x}}_k\right){\mathbf{x}}_k-F\left({\mathbf{x}}_k\right)$$
然后最小二乘法一下
$$\begin{array}{rl}{\mathbf{x}}_{k+1}& ={\left(J{\left({\mathbf{x}}_k\right)}^{T}J\left({\mathbf{x}}_k\right)\right)}^{-1}J\left({\mathbf{x}}_k\right){\mathbf{b}}_k\\ & ={\left(J{\left({\mathbf{x}}_k\right)}^{T}J\left({\mathbf{x}}_k\right)\right)}^{-1}J\left({\mathbf{x}}_k\right)\left(J\left({\mathbf{x}}_k\right){\mathbf{x}}_k-F\left({\mathbf{x}}_k\right)\right)\\ & ={\mathbf{x}}_k-{\left(J{\left({\mathbf{x}}_k\right)}^{T}J\left({\mathbf{x}}_k\right)\right)}^{-1}J{\left({\mathbf{x}}_k\right)}^{T}F\left({\mathbf{x}}_k\right)\end{array}$$
于是下降方向${\mathbf{d}}_k={\left(J{\left({\mathbf{x}}_k\right)}^{T}J\left({\mathbf{x}}_k\right)\right)}^{-1}J{\left({\mathbf{x}}_k\right)}^{T}F\left({\mathbf{x}}_k\right)$
注意到$\nabla g\left(\mathbf{x}\right)=2J{\left(\mathbf{x}\right)}^{T}F\left(\mathbf{x}\right)$
$${\mathbf{d}}_k=\frac{1}{2}{\left(J{\left({\mathbf{x}}_k\right)}^{T}J\left({\mathbf{x}}_k\right)\right)}^{-1}\nabla g\left({\mathbf{x}}_k\right)$$
这个迭代时没有步长的，如果用了步长就是阻尼高斯-牛顿法（Damped Gauss-Newton Method）

具体步骤：
输入：$ϵ>0$
初始化：选择任意${\mathbf{x}}_0\in {\mathbb{R}}^n$
步骤：
(a)${\mathbf{d}}_k={\left(J{\left({\mathbf{x}}_k\right)}^{T}J\left({\mathbf{x}}_k\right)\right)}^{-1}J{\left({\mathbf{x}}_k\right)}^{T}F\left({\mathbf{x}}_k\right)$
(b)选择步长$t_k$
$$h\left(t\right)=g\left({\mathbf{x}}_k-t{\mathbf{d}}_k\right)$$
$(c)$${\mathbf{x}}_{k+1}={\mathbf{x}}_k-t_k{\mathbf{d}}_k$
(d)如果$\Vert \nabla g\left({\mathbf{x}}_{k+1}\right)\Vert \le ϵ$，就停止，并输出${\mathbf{x}}_{k+1}$

代码：这里步长用回溯法选的

function [x,fun_val]=damped_Gauss_Newtow(g,grad,J,F,x0,s,alpha,...
beta,epsilon)
% Gradient method with backtracking stepsize rule
%
% INPUT
%=======================================
% g ......... objective function
% grad ...... gradient of the objective function
% J ......... Jacobian matrix
% F ......... vector-valued function
% x0......... initial point
% s ......... initial choice of stepsize
% alpha ..... tolerance parameter for the stepsize selection
% beta ...... the constant in which the stepsize is multiplied
%             at each backtracking step (0<beta<1)
% epsilon ... tolerance parameter for stopping rule
% OUTPUT
%=======================================
% x ......... optimal solution (up to a tolerance)
%             of min f(x)
% fun_val ... optimal function value
x=x0;
J_val=J(x);
F_val=F(x);
d=(J_val'*J_val)\(J_val'*F_val);
fun_val=g(x);
gval=grad(x);
iter=0;
while (norm(gval)>epsilon&&(iter<10000))
    iter=iter+1;
    t=s;
    while (fun_val-g(x-t*d)<alpha*t*norm(d)^2)
        t=beta*t;
    end
    x=x-t*d;
    J_val=J(x);
    F_val=F(x);
    d=(J_val'*J_val)\(J_val'*F_val);
    fun_val=g(x);
    gval=grad(x);
    fprintf('iter_number = %3d norm_grad = %2.6f fun_val = %2.6f \n',...
        iter,norm(gval),fun_val);
end
if (iter==10000)
    fprintf('did not converge\n')
end