【SVM】c++实现

SVM使用c++实现

下面是使用c++实现一个简单的支持向量机（SVM）分类器，并提供了从文件中读取数据、训练模型和对测试数据进行预测的功能。
SVM 实现是基于线性核函数的，并采用了简化的序列最小优化（SMO）算法进行训练

SVM推导

SVM 旨在找到一个超平面，使得正负样本间的间隔最大化。

对于线性可分的数据集，SVM 的目标函数和约束条件可以写为：

目标函数（最小化）：
1. 12∥w∥2\frac{1}{2} \| \mathbf{w} \|^2 21∥w∥2 其中 www 是超平面的法向量。
约束条件：其中 yiy_iyi是第i个样本的类标签，xix_ixi是第i个样本的特征向量，bbb是偏置项。
1. yi(w⋅xi+b)≥1, ∀iy_i (\mathbf{w} \cdot \mathbf{x}_i + b) \geq 1, \; \forall i yi(w⋅xi+b)≥1,∀i

通过拉格朗日乘子法引入拉格朗日乘子 αi≥0\alpha_i \geq 0αi≥0，得到拉格朗日函数：

L(w,b,α)=12∥w∥2−∑i=1mαi[yi(w⋅xi+b)−1]L(\mathbf{w}, b, \boldsymbol{\alpha}) = \frac{1}{2} \| \mathbf{w} \|^2 - \sum_{i=1}^m \alpha_i [y_i (\mathbf{w} \cdot \mathbf{x}_i + b) - 1] L(w,b,α)=21∥w∥2−i=1∑mαi[yi(w⋅xi+b)−1]

SMO 算法

SMO 算法的目标是求解这个拉格朗日函数的对偶问题，即最大化 $L(\mathbf{w}, b, \boldsymbol{\alpha}) $关于 α\boldsymbol{\alpha}α 的值，同时满足约束条件 αi≥0\alpha_i \geq 0αi≥0 和 ∑i=1mαiyi=0\sum_{i=1}^m \alpha_i y_i = 0∑i=1mαiyi=0，关键是每次只优化两个拉格朗日乘子，因此固定其他乘子。

SMO 算法的主要步骤包括：

选择一对需要优化的拉格朗日乘子 alphaialpha_ialphai 和 αj\alpha_jαj。
固定 αi\alpha_iαi 和 αj\alpha_jαj 以外的所有乘子，只对这两个乘子进行优化。
计算新的 αi\alpha_iαi 和 αj\alpha_jαj。
更新偏置项 bbb。
重复这个过程，直到所有的拉格朗日乘子满足 KKT 条件。

KKT 条件

在 SVM 的背景下，KKT（Karush-Kuhn-Tucker）条件是确定是否找到最优解的关键。对于每个样本 iii，KKT 条件为：

$ \alpha_i = 0 \Rightarrow y_i (\mathbf{w} \cdot \mathbf{x}_i + b) \geq 1 $
$ \alpha_i > 0 \Rightarrow y_i (\mathbf{w} \cdot \mathbf{x}_i + b) = 1 $
$ \alpha_i = C \Rightarrow y_i (\mathbf{w} \cdot \mathbf{x}_i + b) \leq 1 $

其中CCC是正则化参数，在软间隔 SVM 中使用。

代码思路

SVM 类

SVM 类是代码的核心，实现了 SVM 的主要功能：

构造函数：接受惩罚参数 C、容忍度 tolerance 和最大迭代次数 max_passes 作为参数。
fit 方法：用于训练 SVM。它接受输入特征 X 和类标签 y，然后使用 SMO 算法更新拉格朗日乘子（alphas）和偏置项（b）。
predict 方法：根据训练后的模型参数对新数据进行预测。
getWeights 方法：计算并返回权重向量，这是由 alphas、y 和 X 组合而成的。
**kernel**计算线性核、f用于计算决策函数的值。
fit 方法的具体流程
1. initialize 方法负责设置初始状态。
2. mainLoop 方法控制整个算法的迭代过程，直到达到 max_passes 或找到最优解。
3. 在 mainLoop 内部，iterateThroughDataSet 负责遍历所有数据点，并对每个数据点调用 examineExample。
4. examineExample 方法检查并选择是否需要优化特定的 alpha 值。它调用 selectJ 来选择一个合适的 j 索引，并调用 updateAlphasAndB 来更新 alpha 值和偏置项 b。
5. updateAlphasAndB 方法进一步分解为更具体的操作，如计算 L 和 H 的边界（computeLH），计算 eta（computeEta），更新 alpha 值（updateAlphaJ 和 updateAlphaI），以及更新偏置项 b（updateB）。

主函数 `main`

从文件读取训练数据：使用 utils::getData 读取 data.txt，得到训练特征 X 和类标签 y。
训练 SVM：创建 SVM 实例，调用 fit 方法进行训练。
输出训练后的模型参数：打印权重和偏置项。
读取测试数据：从 test.txt 中读取测试数据（仅特征，无类标签）。
对测试数据进行预测并规范化结果：对每个测试样本进行预测，如果预测值大于 0.7 则归类为 1，小于 -0.7 则归类为 -1，否则保持原值。

SVM 训练过程

在 fit 方法中，SVM 使用 SMO 算法来找到最佳的拉格朗日乘子。SMO 算法通过迭代方式优化这些乘子，以最小化损失函数。每一步迭代中，算法选取一对乘子进行优化及计算和更新偏置项 b，以及处理乘子的上下界。

1. 构造函数

初始化 SVM 的正则化参数 $ C $、容忍度（tolerance）、和最大迭代次数（max_passes）。

SVM(double C, double tolerance, int max_passes)
: C(C), tolerance(tolerance), max_passes(max_passes) {}

2. 获取权重

对应步骤:

计算权重向量 w\mathbf{w}w。
公式: w=∑i=1mαiyixi\mathbf{w} = \sum_{i=1}^{m} \alpha_{i} y_{i} \mathbf{x}_{i}w=∑i=1mαiyixi。

vector<double> getWeights() const {
    vector<double> weights(n, 0.0);
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            weights[j] += alphas[i] * y[i] * X[i][j];
        }
    }
    return weights;
}

3. 预测方法

对应步骤:

使用训练好的模型进行预测。
公式: f(x)=∑i=1mαiyiK(x,xi)+bf(\mathbf{x}) = \sum_{i=1}^{m} \alpha_{i} y_{i} K(\mathbf{x}, \mathbf{x}_{i}) + bf(x)=∑i=1mαiyiK(x,xi)+b。

double predict(const vector<double>& x) {
    double sum = 0.0;
    for (int i = 0; i < m; i++) {
        sum += alphas[i] * y[i] * kernel(x, X[i]);
    }
    return sum + b;
}

4. 初始化操作

设置 SVM 的训练数据和参数。

void initialize(const vector<vector<double>>& X, const vector<double>& y) {
    this->X = X;
    this->y = y;
    m = X.size();
    n = X[0].size();
    alphas.resize(m, 0.0);
    b = 0.0;
}

5. 主循环

对应步骤:

SMO 算法的主循环，迭代优化 α 值。

void mainLoop() {
    int passes = 0;
    while (passes < max_passes) {
        int num_changed_alphas = iterateThroughDataSet();
        if (num_changed_alphas == 0) passes++;
        else passes = 0;
    }
}

6. 数据集迭

对应步骤:

遍历整个数据集来优化 α 值。

int iterateThroughDataSet() {
    int num_changed_alphas = 0;
    for (int i = 0; i < m; i++) {
        num_changed_alphas += examineExample(i);
    }
    return num_changed_alphas;
}

7. 检查单个样本

对应步骤:

检查和选择 α 值是否需要优化。
检查 KKT 条件。

int examineExample(int i) {
    double Ei = f(X[i]) - y[i];
    if ((y[i] * Ei < -tolerance && alphas[i] < C) || (y[i] * Ei > tolerance && alphas[i] > 0)) {
        int j = selectJ(i, Ei);
        if (j >= 0) {
            return updateAlphasAndB(i, j);
        }
    }
    return 0;
}

8. 计算 L 和 H 的边界

对应步骤:

计算 α 值的剪辑边界。

// 计算边界 L 和 H
void computeLH(int i, int j, double& L, double& H) {
    if (y[i] != y[j]) {
        L = max(0.0, alphas[j] - alphas[i]);
        H = min(C, C + alphas[j] - alphas[i]);
    } else {
        L = max(0.0, alphas[i] + alphas[j] - C);
        H = min(C, alphas[i] + alphas[j]);
    }
}

9. 计算 eta

对应步骤:

计算 α 值优化时的步长。

// 计算 eta
double computeEta(int i, int j) {
    return 2.0 * kernel(X[i], X[j]) - kernel(X[i], X[i]) - kernel(X[j], X[j]);
}

10. 更新 alphaialpha_ialphai 和 alphajalpha_jalphaj

对应步骤:

更新 α 值和偏置项 b。
根据优化问题的二次规划来更新 α 值。

// 更新 alpha_i 和 alpha_j
int updateAlphasAndB(int i, int j) {
    double Ei = f(X[i]) - y[i];
    double Ej = f(X[j]) - y[j];
    double alpha_i_old = alphas[i], alpha_j_old = alphas[j];

    double L, H;
    computeLH(i, j, L, H);
    if (L == H) return 0;

    double eta = computeEta(i, j);
    if (eta >= 0) return 0;

    int changed = 0;
    changed += updateAlphaJ(i, j, Ei, Ej, eta, L, H);
    changed += updateAlphaI(i, j, alpha_j_old);

    updateB(i, j, Ei, Ej, alpha_i_old, alpha_j_old);
    return changed;
}

11. 训练方法

对应步骤:

初始化并运行 SMO 算法来训练 SVM。

// 训练
void fit(const vector<vector<double>>& X, const vector<double>& y) {
    // 初始化操作
    initialize(X, y);
    // 主循环，迭代优化 alpha 值
    mainLoop();
}

测试

数据

数据如下（traindata.txt）

2.0 3.0 1.0
1.0 5.0 1.0
2.5 7.0 1.0
1.5 2.0 1.0
3.0 3.5 1.0
2.0 4.5 1.0
3.5 5.0 1.0
1.0 6.0 1.0
4.0 2.5 -1.0
2.5 3.0 1.0
3.0 6.0 1.0
4.5 4.0 -1.0
5.0 5.5 -1.0
3.0 1.0 -1.0
3.6 0.5 -1.0

输出结果如下

完整代码

#include <iostream>
#include <vector>
#include <cmath>
#include <limits>
#include <fstream>
#include <random>
#include <algorithm> 
using namespace std;

// 假设每个样本有2个输入特征和1个输出
const size_t INNODE = 2;
const size_t OUTNODE = 1;

struct Sample {
    vector<double> in; // 输入特征
    double out;        // 输出（类标签）
};

namespace utils {
    vector<double> getFileData(const string& filename) {
        vector<double> res;
        ifstream in(filename);
        if (in.is_open()) {
            double buffer;
            while (in >> buffer) {
                res.push_back(buffer);
            }
            in.close();
        } else {
            cout << "Error in reading " << filename << endl;
        }
        return res;
    }

    // 读取数据,includeOutData为true时，读取的数据包括输入和输出，否则只读取输入
    vector<Sample> getData(const string& filename, bool includeOutData = true) {
        vector<Sample> res;
        vector<double> buffer = getFileData(filename);

        size_t stepSize = includeOutData ? (INNODE + OUTNODE) : INNODE;

        for (size_t i = 0; i < buffer.size(); i += stepSize) {
            Sample tmp;
            for (size_t t = 0; t < INNODE; t++) {
                tmp.in.push_back(buffer[i + t]);
            }
            if (includeOutData && i + INNODE < buffer.size()) {
                tmp.out = buffer[i + INNODE];
            }
            res.push_back(tmp);
        }

        return res;
    }

}

class SVM {
public:
    SVM(double C, double tolerance, int max_passes) 
    : C(C), tolerance(tolerance), max_passes(max_passes) {}
    // 在SVM类中添加
    vector<double> getWeights() const {
        vector<double> weights(n, 0.0);
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                weights[j] += alphas[i] * y[i] * X[i][j];
            }
        }
        return weights;
    }
    // 预测方法
    double predict(const vector<double>& x) {
        double sum = 0.0;
        for (int i = 0; i < m; i++) {
            sum += alphas[i] * y[i] * kernel(x, X[i]);
        }
        return sum + b;
    }

    // 获取训练后的 alphas
    const vector<double>& getAlphas() const {
        return alphas;
    }

    // 获取训练后的偏置项 b
    double getBias() const {
        return b;
    }
    // 初始化操作
    void initialize(const vector<vector<double>>& X, const vector<double>& y) {
        this->X = X;
        this->y = y;
        m = X.size();
        n = X[0].size();
        alphas.resize(m, 0.0);
        b = 0.0;
    }
    // 主循环，迭代优化 alpha 值
    void mainLoop() {
        int passes = 0;
        while (passes < max_passes) {
            int num_changed_alphas = iterateThroughDataSet();
            if (num_changed_alphas == 0) passes++;
            else passes = 0;
        }
    }
    // 遍历数据集
    int iterateThroughDataSet() {
        int num_changed_alphas = 0;
        for (int i = 0; i < m; i++) {
            num_changed_alphas += examineExample(i);
        }
        return num_changed_alphas;
    }
    // 检查第 i 个样本是否满足 KKT 条件
    int examineExample(int i) {
        double Ei = f(X[i]) - y[i];
        if ((y[i] * Ei < -tolerance && alphas[i] < C) || (y[i] * Ei > tolerance && alphas[i] > 0)) {
            int j = selectJ(i, Ei);
            if (j >= 0) {
                return updateAlphasAndB(i, j);
            }
        }
        return 0;
    }
    // 随机选择第二个乘子 j
    int selectJ(int i, double Ei) {
        int j = rand() % m;
        while (j == i) {
            j = rand() % m;
        }
        return j;
    }
    // 计算边界 L 和 H
    void computeLH(int i, int j, double& L, double& H) {
        if (y[i] != y[j]) {
            L = max(0.0, alphas[j] - alphas[i]);
            H = min(C, C + alphas[j] - alphas[i]);
        } else {
            L = max(0.0, alphas[i] + alphas[j] - C);
            H = min(C, alphas[i] + alphas[j]);
        }
    }
    // 计算 eta
    double computeEta(int i, int j) {
        return 2.0 * kernel(X[i], X[j]) - kernel(X[i], X[i]) - kernel(X[j], X[j]);
    }
    // 更新 alpha_j
    int updateAlphaJ(int i, int j, double Ei, double Ej, double eta, double L, double H) {
        double alpha_j_old = alphas[j];
        alphas[j] -= y[j] * (Ei - Ej) / eta;
        alphas[j] = min(H, max(L, alphas[j]));

        if (abs(alphas[j] - alpha_j_old) < 1e-5) {
            return 0; // 没有足够的改变
        }
        return 1; // 发生了改变
    }
    // 更新 alpha_i
    int updateAlphaI(int i, int j, double alpha_j_old) {
        double alpha_i_old = alphas[i];
        alphas[i] += y[i] * y[j] * (alpha_j_old - alphas[j]);

        return (alphas[i] != alpha_i_old); // 如果 alpha_i 改变了返回 1，否则返回 0
    }
    // 更新阈值 b
    void updateB(int i, int j, double Ei, double Ej, double alpha_i_old, double alpha_j_old) {
        double b1 = b - Ei - y[i] * (alphas[i] - alpha_i_old) * kernel(X[i], X[i])
                        - y[j] * (alphas[j] - alpha_j_old) * kernel(X[i], X[j]);
        double b2 = b - Ej - y[i] * (alphas[i] - alpha_i_old) * kernel(X[i], X[j])
                        - y[j] * (alphas[j] - alpha_j_old) * kernel(X[j], X[j]);

        if (0 < alphas[i] && alphas[i] < C) b = b1;
        else if (0 < alphas[j] && alphas[j] < C) b = b2;
        else b = (b1 + b2) / 2.0;
    }
    // 更新 alpha_i 和 alpha_j
    int updateAlphasAndB(int i, int j) {
        double Ei = f(X[i]) - y[i];
        double Ej = f(X[j]) - y[j];
        double alpha_i_old = alphas[i], alpha_j_old = alphas[j];

        double L, H;
        computeLH(i, j, L, H);
        if (L == H) return 0;

        double eta = computeEta(i, j);
        if (eta >= 0) return 0;

        int changed = 0;
        changed += updateAlphaJ(i, j, Ei, Ej, eta, L, H);
        changed += updateAlphaI(i, j, alpha_j_old);

        updateB(i, j, Ei, Ej, alpha_i_old, alpha_j_old);
        return changed;
    }
    // 训练
    void fit(const vector<vector<double>>& X, const vector<double>& y) {
        // 初始化操作
        initialize(X, y);
        // 主循环，迭代优化 alpha 值
        mainLoop();
    }

private:
    double C;
    double tolerance;
    int max_passes;
    vector<vector<double>> X;
    vector<double> y;
    vector<double> alphas;
    double b;
    int m, n;

    double kernel(const vector<double>& x1, const vector<double>& x2) {
        // 线性核函数
        double result = 0.0;
        for (int i = 0; i < n; i++) {
            result += x1[i] * x2[i];
        }
        return result;
    }

    double f(const vector<double>& x) {
        double sum = 0.0;
        for (int i = 0; i < m; i++) {
            sum += alphas[i] * y[i] * kernel(x, X[i]);
        }
        return sum + b;
    }
};

int main() {
    // 从文件读取数据
    string filename = "traindata.txt";
    vector<Sample> data = utils::getData(filename);
    
    vector<vector<double>> X;
    vector<double> y;

    for (const auto& sample : data) {
        X.push_back(sample.in);
        y.push_back(sample.out);
    }

    // 创建 SVM 实例并使用读取的数据进行训练
    SVM svm(1.0, 0.001, 10);
    svm.fit(X, y);
    
    // 输出训练后的模型参数
    cout << "\nTrained Model Parameters:" << endl;
    cout << "Weights:" << endl;
    vector<double> weights = svm.getWeights();
    for (const auto& w : weights) {
        cout << w << " ";
    }
    cout << "\nBias (b): " << svm.getBias() << endl;

   // 读取测试数据
    string test_filename = "testdata.txt";
    vector<Sample> test_data = utils::getData(test_filename, false);
    vector<vector<double>> test_X;
    for (const auto& sample : test_data) {
        test_X.push_back(sample.in);
    }

    // 对测试数据进行预测并规范化结果
    cout << "\nTest Data Classification Results (Normalized):" << endl;
    for (const auto& x : test_X) {
        double prediction = svm.predict(x);
        // 规范化预测结果
        double normalized_prediction = (prediction > 0.7) ? 1 : (prediction <-0.7 ? -1: prediction);
        cout << "Test Data Point: (" << x[0] << ", " << x[1] << ") - Predicted Class: " << normalized_prediction << endl;
    }

    return 0;
}

SVM使用c++实现

SVM推导

SMO 算法

KKT 条件

代码思路

SVM 类

主函数 main

SVM 训练过程

1. 构造函数

2. 获取权重

3. 预测方法

4. 初始化操作

5. 主循环

6. 数据集迭

7. 检查单个样本

8. 计算 L 和 H 的边界

9. 计算 eta

10. 更新 alphaialpha_ialphai​ 和 alphajalpha_jalphaj​

11. 训练方法

测试

数据

完整代码

相关推荐

主函数 `main`

10. 更新 alphaialpha_ialphai 和 alphajalpha_jalphaj