Arzie100 发表于 2023-2-18 09:05

常见优化算法实现

一、梯度下降

展开的是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 =
    # print(x_list)
    y_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 =
      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+=commom_ele*x_i*x_i
            J_T+=commom_ele*x_i
            J_T+=commom_ele
      
      delta_a = 1*lr*J_T
      delta_b = 1*lr*J_T
      delta_c = 1*lr*J_T

      # 更新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 =
    # print(x_list)
    # y_list =
    y_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 =
      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+=x_i*x_i*commom_ele
            # J_T+=x_i*commom_ele
            J_T+=1*commom_ele
      
      delta_a = -1*lr*J_T
      # delta_b = -1*lr*J_T
      delta_c = -1*lr*J_T

      # 更新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;
    int b = params;
    int c = params;

    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;
    double& b = parms;
    double& c = parms;
    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算法
页: [1]
查看完整版本: 常见优化算法实现