常见优化算法实现

Arzie100

一、梯度下降

展开的是loss^2的导数
不知道这个函数： y = exp(1*x^{2}+2*x+1) 为什么会不能优化。。。
'''
Description: 梯度下降法优化非线性函数
Author: suyunzheng
Date: 2021-10-06 23:06:34
LastEditTime: 2021-10-06 23:38:21
LastEditors: maple
'''

import os
import sys
import numpy as np
import math
import random

def generator_sample(a, b, c):
    x_list = [i/100 for i in range(0, 100)]
    # print(x_list)
    y_list = [math.exp(a*x*x+b*x+c)+random.random()/10 for x in x_list]
    print(y_list)

    return x_list, y_list
def GD_optimizer(x_list, y_list, epoch, a, b, c, lr):

    print("a:{}, b:{}, c:{}".format(a, b, c))
    for i in range(epoch):
        J_T = [0, 0, 0]
        loss = 0
        # pass
        for (x_i, y_i) in zip(x_list, y_list):
            y_hat_i = math.exp(a*x_i*x_i+b*x_i+c)
            # print("y_i:{}, y_hat_i:{}".format(y_i, y_hat_i))
            loss_i = y_i - y_hat_i
            loss+=loss_i
            # print("loss_:{}".format(loss_i))
            commom_ele = -2*loss_i*y_hat_i
            J_T[0]+=commom_ele*x_i*x_i
            J_T[1]+=commom_ele*x_i
            J_T[2]+=commom_ele
        
        delta_a = 1*lr*J_T[0]
        delta_b = 1*lr*J_T[1]
        delta_c = 1*lr*J_T[2]

        # 更新a b c
        a+=delta_a
        b+=delta_b
        c+=delta_c
        print("Iter: {}, a: {}, b:{}, c:{}, loss: {}".format(i+1, a, b, c, loss))
    pass

   


if __name__ == '__main__':
    x_list, y_list = generator_sample(1,2,1)
    a = 1
    b = 1
    c = 1
    GD_optimizer(x_list, y_list, 100, a, b, c, 0.1)下面对 y=x^{2}+1 优化：
'''
Description: 梯度下降法优化非线性函数
Author: suyunzheng
Date: 2021-10-06 23:06:34
LastEditTime: 2021-10-07 14:57:21
LastEditors: maple
'''

import os
import sys
import numpy as np
import math
import random
import matplotlib.pyplot as plt

def show_x_y(x_list, y_list):
    plt.scatter(x_list, y_list, s = 0.5)
    plt.show()

def generator_sample(a, b, c):
    x_list = [i/100 for i in range(0, 100)]
    # print(x_list)
    # y_list = [math.exp(a*x*x+b*x+c)+random.random() for x in x_list]
    y_list = [a*x*x+b*x+c+random.random()/100 for x in x_list]
   
    # print(y_list)

    show_x_y(x_list, y_list)


    return x_list, y_list
def GD_optimizer(x_list, y_list, epoch, a, b, c, lr):

    print("a:{}, b:{}, c:{}".format(a, b, c))
    for i in range(epoch):
        J_T = [0, 0, 0]
        loss = 0
        # pass
        for (x_i, y_i) in zip(x_list, y_list):
            # y_hat_i = math.exp(a*x_i*x_i+b*x_i+c)
            y_hat_i = a*x_i*x_i+0*x_i+c
            # print("y_i:{}, y_hat_i:{}".format(y_i, y_hat_i))
            # pass
            loss_i = math.pow((y_hat_i - y_i),2)
            loss+=loss_i
            # print("loss_:{}".format(loss_i))
            # commom_ele = 2*(y_hat_i - y_i)*y_hat_i
            commom_ele = 2*(y_hat_i - y_i)
            # commom_ele = y_hat_i
            J_T[0]+=x_i*x_i*commom_ele
            # J_T[1]+=x_i*commom_ele
            J_T[2]+=1*commom_ele
        
        delta_a = -1*lr*J_T[0]
        # delta_b = -1*lr*J_T[1]
        delta_c = -1*lr*J_T[2]

        # 更新a b c
        a+=delta_a
        # b+=delta_b
        c+=delta_c
        print("Iter: {}, a: {}, b:{}, c:{}, lr:{}, loss: {}".format(i+1, a, b, c, lr, loss))
    pass

   


if __name__ == '__main__':
    x_list, y_list = generator_sample(1,0,1)
    a = 0.1
    b = 0.0
    c = 0.3
    GD_optimizer(x_list, y_list, 1000, a, b, c, 0.0001)二、牛顿法

三、高斯牛顿法

参考：
/*
* @Description: 常见优化算法的实现
* @Author: suyunzheng
* @Date: 2021-10-05 21:50:27
* @LastEditTime: 2021-10-05 22:55:12
* @LastEditors: maple
*/

#include <iostream>
#include <vector>
#include <eigen3/Eigen/Core>
#include <eigen3/Eigen/Dense>


using namespace std;

constexpr int N = 100;
void generator_sample(const vector<int>& params, vector<double>& x_data, vector<double>& y_data){
    int a = params[0];
    int b = params[1];
    int c = params[2];

    for(int i = 0; i<N; i++){
        double x = i/100.0;
        x_data.push_back(x);
        y_data.push_back(exp(a*x*x+b*x+c)+(rand()%100)/100);
    }
   
    for(int i = 0; i<N; i++){
        cout << x_data << ", " << y_data << endl;
    }
    cout << endl;
   
}

int GN(vector<double>& parms, vector<double>& x_data, vector<double>& y_data){
    int iter_num = 100;
    double& a = parms[0];
    double& b = parms[1];
    double& c = parms[2];
    double pre_cost = 0;
    double cost = INT_MAX;

    for(int i = 0; i<iter_num; i++){
        double eps = pre_cost - cost;
        cout << "Iter: " << i <<  ", cost: " << cost << ", eps: " << eps << endl;
        if(cost<=1e-10 || abs(eps)<=1e-10 && i>=30){
            cout << "convergence" << endl;
            break;
        }
        pre_cost = cost;
        cost = 0;           // cost清零
        Eigen::Matrix3d H_ = Eigen::Matrix3d::Zero();
        Eigen::Vector3d b_ = Eigen::Vector3d::Zero();
        for(int i = 0; i<x_data.size(); i++){
            Eigen::Matrix3d tmp_H = Eigen::Matrix3d::Zero();
            Eigen::Vector3d tmp_b = Eigen::Vector3d::Zero();
            Eigen::Vector3d J_T = Eigen::Vector3d::Zero();
            double x = x_data;
            double y = y_data;
      
            double exp_ele = exp(a*x*x+ b*x + c );
            J_T.x() = -exp_ele*x*x;
            J_T.y() = -exp_ele*x;
            J_T.z() = -exp_ele;
            
            // H和b求和
            H_+=J_T*J_T.transpose();  ;      
            b_+=-1*J_T*(y-exp_ele);
            cost+=pow(y-exp_ele, 2);
            
        }
        
        // cout << H_ << endl;
        // Eigen::Vector3d delta = H_.inverse()*b_;
        // Eigen::Vector3d delta = H_*b_;
        Eigen::Vector3d delta = H_.ldlt().solve(b_);

        a+=delta.x();
        b+=delta.y();
        c+=delta.z();
        // cout << "a: " << a << ", b: " << b << ", c: " << c << endl;
    }

    cout << "a: " << a << ", b: " << b << ", c: " << c << endl;
   
}


int main(){

    vector<int> params = {1,2,1};
    vector<double> x_data;
    vector<double> y_data;
   
    generator_sample(params, x_data, y_data);

    vector<double> guess_params = {0, 0 , 0};
   
    // 优化
    GN(guess_params, x_data, y_data);


    return 0;
}



四、L-M算法