最优化方法：无约束梯度问题----最速下降法理论与C++实现

背景

优化问题的关键在于每次更新时方向、步长的选择（比较简单的方式是固定步长，缺陷步长选太小收敛速度慢，太大则在靠近最优值时剧烈震荡）。

最速下降方法：方向选择的是负梯度方向，步长通过极小化 g(k) = f(x0+k*dx) 确定（这里 g(k) 是关于k这个常数的一元函数，可以通过求导等于零求解析解，也可以使用一元搜索方法：斐波那契、黄金分割法迭代计算）。

由于函数的二阶泰勒展开式（忽略无穷小项）形如：f(x) = 1/2*X'AX+X'B+C   (其中X是nx1的列向量，X’是X的转置，A是nxn矩阵，B是nx1列向量，C是常数)，所以以解二阶多项式函数最优化问题为例

如何求解梯度、确定步长参考https://www.bilibili.com/video/av14666879/?p=9

最优梯度的缺陷是不能在有限步达到最优值，共轭梯度下降在理论上可以保证N步就收敛（N是未知变量X的维度）

问题

                 min 2*x^2+2*x*y + 5*y^2      初始点 z = [2,-2]'

 

代码实现

#include<iostream>
#include<vector>
using namespace std;

vector<vector<double>> A = {{4.0,2.0},{2.0,10.0}};
vector<double> B = {0.0,0.0};
double C = 0.0;

vector<double> mat_vect(const vector<vector<double>>& arr,const vector<double>& vect){
    int m = arr.size(),n = arr[0].size();
    vector<double> ans;
    for(int i=0;i<m;++i){
        for(int j=0;j<1;++j){
            double tmp = 0;
            for(int p=0;p<n;++p) tmp += arr[i][p]*vect[p];
            ans.push_back(tmp);
        }
    }
    return ans;
}

vector<double> mat_vect(const vector<double>& vect,const vector<vector<double>>& arr){
    int a_n = arr.size(),k = arr[0].size();
    vector<double> ans;
    for(int i=0;i<1;++i){
        for(int j=0;j<k;++j){
            double tmp = 0;
            for(int p=0;p<a_n;++p) tmp += vect[p]*arr[p][j];
            ans.push_back(tmp);
        }
    }
    return ans;
}

double mat_double(const vector<double>& vect1,const vector<double>& vect2){
    double ans = 0.0;
    for(int i=0;i<vect1.size();++i) ans += vect1[i]*vect2[i];
    return ans;
}

double get_func(const vector<double>& x){
    double ans = 0;
    auto tmp = mat_vect(x,A);
    ans = mat_double(tmp,x);
    ans/=2.0;
    ans += mat_double(x,B);
    return ans+C;
}

vector<double> get_grad(const vector<double>& x){
    auto tmp = mat_vect(A,x);
    for(int i=0;i<x.size();++i) tmp[i]+=B[i];
    return tmp;
}

double get_k_star(const vector<double>& x,vector<double>& grad){
    grad = get_grad(x);
    double ans = mat_double(grad,grad);
    auto vect = mat_vect(grad,A);
    double tmp = mat_double(vect,grad);
    return -(ans/tmp);
}

void upgrade(vector<double>& x){
    vector<double> grad;
    double k = get_k_star(x,grad);
    cout<<"k_star: "<<k<<endl<<"X: ";
    for(int i=0;i<x.size();++i){
        x[i] = x[i]+k*grad[i];
        cout<<x[i]<<" ";
    }
    cout<<endl;
    cout<<"f(x): "<<get_func(x)<<endl;
}

int main(){
    vector<double> x = {2.0,-2.0};
    cout<<"f(x): "<<get_func(x)<<endl;
    for(int i=1;i<=8;++i){
        cout<<"----------------------------------"<<i<<endl;
        upgrade(x);
    }
    return 0;
}

运行结果

运行结果与视频中每次迭代求得的值存在误差，可能是计算误差也可能是代码存在问题。期待并欢迎大佬指出