C++ 逻辑回归模型求解 (基于Eigen)
1. 目标函数:sigmoid函数
f_wb=1/[1+exp(-z)]
z=wx+b,w和x是n维向量
/**
* @author @还下着雨ZG
* @brief sigmoid函数,用于实现逻辑回归
* @param[in] z
* @return 1/(1+exp(-z))
*/
double sigmoid(const double &z){
return 1/(1+exp(-z));
}
2. 损失函数:
3. 代价函数:
计算模型的代价:
/**
* @author @还下着雨ZG
* @brief 计算逻辑回归模型的代价cost
* @param[in] X
* @param[in] y
* @param[in] w_in
* @param[in] b_in
* @return cost
*/
double compute_cost(MatrixXd X,VectorXi &y,VectorXd &w_in,double b_in){
int m=X.rows();
double cost=0;
for(int i=0;i<m;++i){
double z=w_in.dot(X.row(i))+b_in;
double f_wb=sigmoid(z);
cost+=-y[i]*log(f_wb)-(1-y[i])*log(1-f_wb);
}
return cost/m;
}
4. 梯度下降法优化模型得到最佳参数
(1)计算代价函数的梯度:
/**
* @author @还下着雨ZG
* @brief 计算逻辑回归模型的梯度
* @param[in] X
* @param[in] y
* @param[in] w_in
* @param[in] b_in
* @return pair(diff_w,diff_b) 可以携带两组参数
*/
std::pair<VectorXd,double> compute_gradient(MatrixXd &X,VectorXi &y,
VectorXd &w_in,double &b_in){
int m=X.rows();int n=X.cols();
VectorXd diff_w=VectorXd::Zero(n);
VectorXd w=w_in;double b=b_in;
double diff_b;
double f_wb=0.0,z=0.0;
for(int i=0;i<m;++i){
z=w.dot(X.row(i))+b;
f_wb=sigmoid(z);
double error=f_wb-y[i];
for(int j=0;j<n;j++) diff_w[j]+=error*X(i,j);
diff_b+=error;
}
diff_w=diff_w/m;diff_b=diff_b/m;
return std::make_pair(diff_w,diff_b);
}
(2) 梯度下降法求解参数:
/**
* @brief 梯度下降求解目标函数sigmoid的最佳参数
* @param[in] X
* @param[in] y
* @param[in] w_init
* @param[in] b_init
* @param[in] learning_rate
* @param[in] iterations 迭代次数
* @return w_out,b_out
*/
std::pair<VectorXd,double> gradient_descent(MatrixXd &X,VectorXi &y,VectorXd &w_init,
double &b_init,double &alpha,int &iterations){
int m=X.rows();int n=X.cols();
VectorXd w(w_init);double b=b_init;
VectorXd diff_w(n);double diff_b=0.0;
for(int i=0;i<iterations;++i){
auto gradient_res=compute_gradient(X,y,w,b);
diff_w=gradient_res.first;diff_b=gradient_res.second;
w=w-alpha*diff_w;
b=b-alpha*diff_b;
}
return std::make_pair(w,b);
}
可直接运行代码实例:
#include <iostream>
#include <Eigen/Dense>
#include <string>
using namespace Eigen;
/**
* @author @还下着雨ZG
* @brief sigmoid函数,用于实现逻辑回归
* @param[in] z
* @return 1/(1+exp(-z))
*/
double sigmoid(const double &z){
return 1/(1+exp(-z));
}
/**
* @author @还下着雨ZG
* @brief 计算逻辑回归模型的代价cost
* @param[in] X
* @param[in] y
* @param[in] w_in
* @param[in] b_in
* @return cost
*/
double compute_cost(MatrixXd X,VectorXi &y,VectorXd &w_in,double b_in){
int m=X.rows();
double cost=0;
for(int i=0;i<m;++i){
double z=w_in.dot(X.row(i))+b_in;
double f_wb=sigmoid(z);
cost+=-y[i]*log(f_wb)-(1-y[i])*log(1-f_wb);
}
return cost/m;
}
/**
* @author @还下着雨ZG
* @brief 计算逻辑回归模型的梯度
* @param[in] X
* @param[in] y
* @param[in] w_in
* @param[in] b_in
* @return pair(diff_w,diff_b) 可以携带两组参数
*/
std::pair<VectorXd,double> compute_gradient(MatrixXd &X,VectorXi &y,
VectorXd &w_in,double &b_in){
int m=X.rows();int n=X.cols();
VectorXd diff_w=VectorXd::Zero(n);
VectorXd w=w_in;double b=b_in;
double diff_b;
double f_wb=0.0,z=0.0;
for(int i=0;i<m;++i){
z=w.dot(X.row(i))+b;
f_wb=sigmoid(z);
double error=f_wb-y[i];
for(int j=0;j<n;j++) diff_w[j]+=error*X(i,j);
diff_b+=error;
}
diff_w=diff_w/m;diff_b=diff_b/m;
return std::make_pair(diff_w,diff_b);
}
/**
* @brief 梯度下降求解目标函数sigmoid的最佳参数
* @param[in] X
* @param[in] y
* @param[in] w_init
* @param[in] b_init
* @param[in] learning_rate
* @param[in] iterations 迭代次数
* @return w_out,b_out
*/
std::pair<VectorXd,double> gradient_descent(MatrixXd &X,VectorXi &y,VectorXd &w_init,
double &b_init,double &alpha,int &iterations){
int m=X.rows();int n=X.cols();
VectorXd w(w_init);double b=b_init;
VectorXd diff_w(n);double diff_b=0.0;
for(int i=0;i<iterations;++i){
auto gradient_res=compute_gradient(X,y,w,b);
diff_w=gradient_res.first;diff_b=gradient_res.second;
w=w-alpha*diff_w;
b=b-alpha*diff_b;
}
return std::make_pair(w,b);
}
int main(){
using std::cout;using std::endl;using std::string;
MatrixXd X_train{{0.5, 1.5}, {1,1}, {1.5, 0.5}, {3, 0.5}, {2, 2}, {1, 2.5}};
int m=X_train.rows();
int n=X_train.cols();
VectorXi y_train(m);
y_train<<0, 0, 0, 1, 1, 1;
VectorXd w_init(n);w_init<<2.0,3.0;
double b_init=1.0;
double learning_rate=0.1;
int iterations=10000;
auto res=gradient_descent(X_train,y_train,w_init,b_init,learning_rate,iterations);
VectorXd w=res.first;
double b=res.second;
cout<<w<<endl;
cout<<b<<endl;
return 0;
}