常见优化算法实现

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

一、梯度下降

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

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(&#34;a:{}, b:{}, c:{}&#34;.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(&#34;y_i:{}, y_hat_i:{}&#34;.format(y_i, y_hat_i))
         loss_i = y_i - y_hat_i
         loss+=loss_i
         # print(&#34;loss_:{}&#34;.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(&#34;Iter: {}, a: {}, b:{}, c:{}, loss: {}&#34;.format(i+1, a, b, c, loss))
pass

if __name__ == &#39;__main__&#39;:
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 优化：
&#39;&#39;&#39;
Description: 梯度下降法优化非线性函数
Author: suyunzheng
Date: 2021-10-06 23:06:34
LastEditTime: 2021-10-07 14:57:21
LastEditors: maple
&#39;&#39;&#39;

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(&#34;a:{}, b:{}, c:{}&#34;.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(&#34;y_i:{}, y_hat_i:{}&#34;.format(y_i, y_hat_i))
         # pass
         loss_i = math.pow((y_hat_i - y_i),2)
         loss+=loss_i
         # print(&#34;loss_:{}&#34;.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(&#34;Iter: {}, a: {}, b:{}, c:{}, lr:{}, loss: {}&#34;.format(i+1, a, b, c, lr, loss))
pass

if __name__ == &#39;__main__&#39;:
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 << &#34;, &#34; << 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 << &#34;Iter: &#34; << i <<  &#34;, cost: &#34; << cost << &#34;, eps: &#34; << eps << endl;
      if(cost<=1e-10 || abs(eps)<=1e-10 && i>=30){
         cout << &#34;convergence&#34; << 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 << &#34;a: &#34; << a << &#34;, b: &#34; << b << &#34;, c: &#34; << c << endl;
}

cout << &#34;a: &#34; << a << &#34;, b: &#34; << b << &#34;, c: &#34; << 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算法

		自动登录	找回密码
密码			立即注册