【非线性方程（组）求解】牛顿-拉夫逊算法应用浅析

【非线性方程求解】割线法
【非线性方程（组）求解】牛顿-拉夫逊算法 多输入单输出

1. 概述

通过以上两篇博客，你已经体会到了牛顿-拉夫逊算法在求解含单/多变量非线性方程的便利，牛顿-拉夫逊算法的另一个广泛应用是无约束优化，一般用于求解无约束目标函数的极小值，本文具体来阐述这个过程，怎么从解非线性方程组变成求极小值的。

2. 原理步骤

2.1 解多元非线性方程组

假设你看过之前内容，知道问题描述和符号定义，了解基本牛顿迭代流程。我把主要公式抄过来：

$$\begin{array}{rl}X=& [x_1,\dots ,x_{nx}{]}^T\\ & \\ Y=& [y_1,\dots ,y_{ny}{]}^T\\ & \\ Y^t=& [y_1^t,\dots ,y_{ny}^t{]}^T\\ & \\ f(X)=& \phi (X)-Y^t=0\end{array}$$

$$\begin{array}{rl}& \to f(X)=f\left(X^k\right)+f^{\prime }\left(X^k\right)\left(X-X^k\right)+f^{\prime \prime }(\xi ){\left(X-X^k\right)}^2/(2!)\\ & \\ & \to f\left(X^k\right)+f^{\prime }\left(X^k\right)\left(X-X^k\right)=0\\ & \\ & \to f^{\prime }\left(X^k\right)\left(X-X^k\right)=-f\left(X^k\right)\\ & \\ & \to \cancel{X-X^k=-\frac{f\left(X^k\right)}{f^{\prime }\left(X^k\right)}}\\ & \\ & \to \cancel{X=X^k-\frac{f\left(X^k\right)}{f^{\prime }\left(X^k\right)}}\end{array}$$

多元非线性方程组含有多个未知量，且方程数不止一个，即 $ny>1$，为了将公式（1）继续推导下去，方程组还需满足条件 $nx=ny$，此时 $f^{\prime }\left(X^k\right)$ 是一个方阵，实际应用需保证可逆，公式（1）变为：

$$\begin{array}{rl}& \to inv[f^{\prime }\left(X^k\right)]\cdot f^{\prime }\left(X^k\right)(X-X^k)=-inv[f^{\prime }\left(X^k\right)]\cdot f\left(X^k\right)\\ & \\ & \to (X-X^k)=-inv[f^{\prime }\left(X^k\right)]\cdot f\left(X^k\right)\\ & \\ & \to dX=-inv[f^{\prime }\left(X^k\right)]\cdot f\left(X^k\right)\end{array}$$
$inv[f^{\prime }\left(X^k\right)]$表示对 $f(X^k)$ 求偏导数后再求逆矩阵。得到增量的表达式后，牛顿-拉夫逊算法求解非线性方程组最终迭代表达式为：
$$\begin{array}{r}{X}^k=X^{k-1}-inv[f^{\prime }\left(X^{k-1}\right)]\cdot f\left(X^{k-1}\right)\end{array}$$
迭代式（2）要求多元非线性方程组的未知量个数等于方程数，换一种说法，系统的输入个数与输出个数一致。

2.2 求目标函数极小值

假设有无约束多维极值问题：
$$\begin{array}{r}minF(X),X\in R^{nx}\end{array}$$
为区别于方程组 $f(X)=0$，目标函数记为 $F(X)$ ，它是一个数。无约束多维极值的一般问题需要求全局最小点，但是大多数优化算法做不到这一点，大多只能找到局部最优点，但是实际问题都有一定 的 应用背景，局部最优点不多，有时甚至局部最优点就是全局最优点，所以多半可以凭借经验来判断结果的可用性¹。于是大家喜欢令 $F(X)$ 的偏导数 $\frac{\partial F(X)}{\partial X}$ 为0：
$$\begin{array}{r}\frac{\partial F(X)}{\partial X}=▽F(X)=\left[\begin{array}{c}\frac{\partial F}{\partial x_1}\\ \frac{\partial F}{\partial x_2}\\ ⋮\\ \frac{\partial F}{\partial x_{nx}}\end{array}\right]=0\end{array}$$
这不巧了吗！公式（4）正好是含有 $nx$ 个未知数和 $nx$ 个方程的非线性方程组。直接套用公式（2）就得到迭代公式：
$$\begin{array}{r}{X}^k=X^{k-1}-[{▽}^{2}F(X^{k-1}){]}^{-1}▽F(X^{k-1})\end{array}$$
其中 ${▽}^{2}F(X^{k-1})$ 是二阶偏导数构成的方阵，著名的黑塞（Hessian）矩阵，计算思路可以参照之前文章。

3. 实例测试

$$\begin{array}{rl}X& =[x_1,x_2{]}^T\\ & \\ f(X)& =\left[\begin{array}{c}(x_1-3{)}^2+(x_2-4{)}^2-25\\ x_1^2+(x_2-4{)}^2-25\end{array}\right]\\ & \\ X^0& =[-0.9,-0.9]\\ & \\ △X& =[0.0001,0.0001]\\ & \\ ϵ& =0.001\end{array}$$

结果如下：

优化历程	优化历程
两个目标值变化历程	变量优化历程

你会发现有趣的现象：迭代初值点处于两个不同曲面上，之后向曲面的相交线移动，逐渐重合，最后沿着相交线迭代到为零的点。

4. 总结

介绍了牛顿-拉夫逊法在求极值问题中的应用。由于实际问题中的求极小值可以被转换成求偏导为零的非线性方程组，而牛顿-拉夫逊法刚好适于求解未知数与方程数相同的方程组，简单又高效，基本上迭代几次即可得解，大家都喜欢用。此外在最优化中优化参数过多时可以将参数分块求解以降低优化的复杂度，需要一些主成分/灵敏度分析，把有些参数和等式约束转换成求解方程组，进而用上牛顿-拉夫逊法，有空会讲到。

5. 详细代码

ExampleNewtonLaphson.m

%
% Auther:  sanfan66
%
close all
clear all
x1=-2:0.5:2;
x2=-2:0.5:2;
[x,y]=meshgrid(x1,x2);
for i=1:size(x,1)
   for j=1:size(x,2) 
    f(i,j,1:2)=TestFunction([x(i,j),y(i,j)]');
   end
end

X0=[-0.9; -0.9];
deltaX0=[0.0001; 0.0001];
epsF=0.001;
try
    [index,X,fx,dx] = NewtonLaphsonMethod(X0,deltaX0,epsF,@TestFunction);
catch ME
     disp(ME.message);
end

figure(1)
mesh(x,y,f(:,:,1))
hold on;
mesh(x,y,f(:,:,2))
hold on;
mesh(x,y,f(:,:,1)*0)
% hold on;
% plot3(X(1,1),X(2,1),fx(1),'b+',X(1,index),X(2,index),fx(index),'g+')
hold on;
plot3(X(1,:),X(2,:),fx(1,:),'r-o')
hold on;
plot3(X(1,:),X(2,:),fx(2,:),'r-o')
hold on;
xCircle=-2:0.1:8;
plot3(xCircle,sqrt(25-(xCircle-3).^2)+4,xCircle*0,'r')
hold on;
plot3(xCircle,-sqrt(25-(xCircle-3).^2)+4,xCircle*0,'r')
xlabel("x1")
ylabel("x2")
zlabel("f")
xlim([-2 3])
ylim([-2 4])

figure
plot(1:index,fx(1,1:index),'-o',1:index,fx(2,1:index),'-+')
xlabel("迭代次数")
legend('fx1','fx2');
grid on;grid minor;
figure
plot(X(1,:),X(2,:),'-o')
xlabel("x1")
ylabel("x2")
grid on;grid minor;

fprintf("%5s%10s%10s%10s%10s\n","i","f0","ff","x1","x2");
for i=1:index
fprintf("%5.f%10f%10f%10f%10f\n",i,fx(1),fx(i),X(1,i),X(2,i));
end

function f=TestFunction(x)
f(1,:)=(x(1,:)-3).^2+(x(2,:)-4).^2-25;
f(2,:)=x(1,:).^2+(x(2,:)-4).^2-25;
end

NewtonLaphsonMethod.m

function [i,X,fx,dx] = NewtonLaphsonMethod(X0,deltaX0,epsF,objFun)
nx=size(X0,1);%自变量的个数
X(1:nx,1)=X0;%每次迭代保存自变量
fx(:,1)=objFun(X(1:nx,1));
ny=size(fx,1);%应变量的个数
diffF(1:ny,1:nx,1)=SolveDiff(X0,deltaX0,objFun);%计算偏导数,矩阵的存放方式与python不同
dx(1:nx,1)=inv(diffF(1:ny,1:nx,1))*(-fx(1:ny,1));%将总变化量均分给两部分,((nx*ny)*(ny*1)).*(nx*1)=nx*1
maxIter=10;
for i=2:maxIter
    X(1:nx,i)=X(1:nx,i-1)+dx(1:nx,i-1);
    fx(1:ny,i)=objFun(X(1:nx,i));
    diffF(1:ny,1:nx,i)=SolveDiff(X(1:nx,i),deltaX0,objFun);
    dx(1:nx,i)=inv(diffF(1:ny,1:nx,i))*(-fx(1:ny,i));%将总变化量均分给两部分
    if(norm(fx(1:ny,i))<=epsF)
        break;
    end
    
end

% if(i==maxIter && norm(fx(:,i))>epsF)
%     ME=MException('NewtonLaphsonError:maxIter', '超过最大迭代次数：i=%d，maxIter=%d',i+1,maxIter);
%     throw(ME);
% end

end

SolveDiff.m

function diffF = SolveDiff(X0,deltaX0,objFun)
%X0列向量
fx0=objFun(X0);
nx=size(X0,1);
ny=size(fx0,1);
deltaX=diag(deltaX0);%增量对角矩阵
XPlus=deltaX+repmat(X0,[1,nx]);%多列X0组成的方阵，并使对角线元素加上小增量
XMinus=-deltaX+repmat(X0,[1,nx]);
fxPlus=objFun(XPlus);%维度变为ny*nx
fxMinus=objFun(XMinus);
diffF=(fxPlus-fxMinus)./repmat(deltaX0',[length(fx0),1])/2;%返回维度ny*nx

end

---

1. 龚纯. 精通MATLAB最优化计算. 2009:电子工业出版社. ↩︎