多维非线性最优化牛顿下山法
对于多元函数极值问题描述如下
\begin{equation}
\min f\left( {\bf{x}} \right)
\end{equation}
其中
\begin{equation}
{\bf{x}} = \left( \begin{array}{l}
{x_1}\\
{x_2}\\
\vdots \\
{x_n}
\end{array} \right)
\end{equation}
多元函数极值问题有多种数值方法求解，本节先介绍一个基础且重要的方法-牛顿法。
牛顿法的基本思想是把函数在局部用泰勒级数的二次展开取近似原函数，
而二次函数的极值问题是容易求得的，故而逐步迭代，最终得到满足精度需求的方法。
对目标函数进行泰勒展开，保留至二次项
\[f\left( {\bf{x}} \right) \approx \phi \left( {\bf{x}} \right) = f\left( {{{\bf{x}}_k}} \right) + {\left[ {\nabla f\left( {{{\bf{x}}_k}} \right)} \right]^T}\left( {{\bf{x}} - {{\bf{x}}_k}} \right) + \frac{1}{2}{\left( {{\bf{x}} - {{\bf{x}}_k}} \right)^T}{\nabla ^2}f\left( {{{\bf{x}}_k}} \right)\left( {{\bf{x}} - {{\bf{x}}_k}} \right)\]
求其驻值
\[\nabla \phi \left( {\bf{x}} \right) = {\bf{0}}\]
即
\[\nabla f\left( {{{\bf{x}}_k}} \right) + {\nabla ^2}f\left( {{{\bf{x}}_k}} \right)\left( {{\bf{x}} - {{\bf{x}}_k}} \right) = {\bf{0}}\]
从而有下面的迭代格式
\[\left\{ \begin{array}{l}
{\bf{H}}\left( {{{\bf{x}}_k}} \right){\bf{v = - }}\nabla {\bf{f}}\left( {{{\bf{x}}_k}} \right)\\
{{\bf{x}}_{k + 1}} = {{\bf{x}}_k} + {\bf{v}}
\end{array} \right.\]
其中$$\nabla {\bf{f}}\left( {{{\bf{x}}_k}} \right)$$为梯度向量
\[\nabla {\bf{f}}\left( {{{\bf{x}}_k}} \right) = \left[ \begin{array}{l}
\frac{{\partial f}}{{\partial {x_1}}}\left( {{x_1}, \cdots ,{x_n}} \right)\\
\quad \quad \vdots \\
\frac{{\partial f}}{{\partial {x_n}}}\left( {{x_1}, \cdots ,{x_n}} \right)
\end{array} \right]\]
$${\bf{H}}\left( {{{\bf{x}}_k}} \right)$$
为Hessian矩阵
\[{\bf{H}} = \left[ {\begin{array}{*{20}{c}}
{\frac{{{\partial ^2}f}}{{\partial {x_1}\partial {x_1}}}}& \cdots &{\frac{{{\partial ^2}f}}{{\partial {x_1}\partial {x_n}}}}\\
\vdots & \ddots & \vdots \\
{\frac{{{\partial ^2}f}}{{\partial {x_n}\partial {x_1}}}}& \cdots &{\frac{{{\partial ^2}f}}{{\partial {x_n}\partial {x_n}}}}
\end{array}} \right]\]
c#版本代码如下。
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;
class cnewton
/*--------------------------------------class comment
Version : V1.0
Coded by : syz
Date : 2011-07-11 08:46:04 *星期一*
----------------------------------------------------
Desciption : 牛顿下山法计算最优化问题
Methods :
* newton--------牛顿下山方法函数
* pivot---------选主元消去法解线性方程
* uptri---------上三角方程回代方法
* func----------目标函数
* grad----------梯度函数
* hess----------hessian矩阵
--------------------------------------------------*/
{
public StreamWriter fout = new StreamWriter("fout.txt");
//输出文件
public void newton(out VEC x, out double fx,VEC x0)
/*------------------------------------------ comment
Author : syz
Date : 2011-07-11 08:52:58
---------------------------------------------------
Desciption : 牛顿下山法计算非线性最优化
* 注意：牛顿法形式上需要计算Hessian矩阵的逆，
* 实际计算时直接解算线性方程
paramtegers :
* x------极小值点
* fx-----极值点的函数值
* x0-----初值向量
*
* ~~~~以下量在程序开头部分控制，也可以把控制量放在参数中，由外部控制
* imax------最大迭代次数
* tol-------误差容限
-------------------------------------------------*/
{
int imax = 50;
//最大允许迭代次数
double tol = 1e-7;
//误差容限
VEC dx, x1, df;
MAT H;
//这里皆不指定维数大小
for (int k = 0; k < imax; k++)
{
grad(out df, x0);
//求的梯度向量
hess(out H, x0);
//求Hessian矩阵
//实际方程为 H dx= -df
//这里转一次负号，重载算符
//df = 0 - df;
pivot(H, df, out dx);
x1 = x0 - dx;
//------------------------
//向量的模，算符重载
double err = ~dx;
//满足精度时退出循环
if (err < tol)
break;
fout.WriteLine("{0,5:D3}{1,12:F6}{2,12:F6}",
k, x1[0], x1[1]);
x0 = x1;
}
x = x0;
func(out fx, x);
}
public void pivot(MAT A, VEC b, out VEC x)
{
/*------------------------------------------ comment
Author : syz
Date : 2011-06-20 10:59:59
---------------------------------------------------
Desciption :
* 消去法解线性方程组
Post Script :
* 要用到上三角方程的回代方法
paramtegers :
* A---方程组系数矩阵
* b---右向量
* x---计算结果
-------------------------------------------------*/
int n, i, j, k;
n = b.dim;
//方程维数
double tmp;
//增广矩阵
MAT Ab = new MAT(n, n + 1);
//赋值给增广矩阵
for (i = 0; i < n; i++)
{
for (j = 0; j < n; j++)
Ab.ele[i, j] = A.ele[i, j];
Ab.ele[i, n] = b.ele[i];
}
double elmax;
int id_max;
double tmp1, tmp2;
//这段是高斯消去法的核心
for (k = 0; k <= n - 2; k++)
{
//主元素初值
elmax = Math.Abs(Ab.ele[k, k]);
//主元素标号
id_max = k;
//下面一段为找出主元素所在的行
//即获得主元素标号的值
for (i = k + 1; i <= n - 1; i++)
{
if (Math.Abs(Ab.ele[i, k]) > elmax)
{
elmax = Math.Abs(Ab.ele[i, k]);
//注意这里把该元素的绝对值复制给最大元素
id_max = i;
}
}
//n+1个元素，该循环为交换两行元素
for (i = 0; i < n + 1; i++)
{
tmp1 = Ab.ele[k, i];
tmp2 = Ab.ele[id_max, i];
Ab.ele[k, i] = tmp2;
Ab.ele[id_max, i] = tmp1;
}
//交换后实施消去
for (i = k + 1; i <= n - 1; i++)
{
tmp = Ab.ele[i, k] / Ab.ele[k, k];
for (j = 0; j < n + 1; j++)
Ab.ele[i, j] = Ab.ele[i, j] - tmp * Ab.ele[k, j];
}
}
//在C#中使用传值方法，不会改变原始变量的值
for (i = 0; i < n; i++)
{
for (j = 0; j < n; j++)
A.ele[i, j] = Ab.ele[i, j];
b.ele[i] = Ab.ele[i, n];
}
uptri(A, b, out x);
}
public void uptri(MAT A, VEC b, out VEC x)
{
/*------------------------------------------ comment
Author : syz
Date : 2011-06-17 10:04:27
---------------------------------------------------
Desciption :
* 上三角矩阵方程解法
Post Script :
* 要用到MAT 与 VEC 类
paramtegers :
* A-----系数矩阵
* b-----右向量
* x-----方程的解
-------------------------------------------------*/
int N;
N = A.dim1;
//获取方程维数
x = new VEC(N);
x.ele[N - 1] = b.ele[N - 1] / A.ele[N - 1, N - 1];
for (int i = N - 2; i >= 0; i--)
{
x.ele[i] = b.ele[i];
for (int j = i + 1; j <= N - 1; j++)
x.ele[i] = x.ele[i] - A.ele[i, j] * x.ele[j];
x.ele[i] = x.ele[i] / A.ele[i, i];
}
}
}