基于C语言写一个SMO算法

以下是一个简单的SMO算法的C语言实现，主要是针对二分类问题的线性SVM模型。在实现时，我们使用了平均学习速率和软间隔方法，以提高算法的鲁棒性和泛化能力。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define MAX_ITER 1000  // 最大迭代次数
#define EPSILON 1e-5   // 精度要求
#define C 1.0          // 惩罚系数
#define TOL 1e-3       // 容忍度

// 数据结构：支持向量
typedef struct {
    double alpha;   // 拉格朗日乘子
    double b;       // 偏置项
    double *x;      // 输入向量
    int y;          // 标签
} SupportVector;

// 内积函数
double dotProduct(double *x1, double *x2, int n) {
    double result = 0.0;
    for (int i = 0; i < n; i++) {
        result += x1[i] * x2[i];
    }
    return result;
}

// 计算误差函数
double errorFunction(double *x, double y, SupportVector *sv, int n) {
    double result = -sv[0].b;
    for (int i = 0; i < n; i++) {
        if (sv[i].alpha > 0) {
            result += sv[i].alpha * sv[i].y * dotProduct(sv[i].x, x, n);
        }
    }
    return result - y;
}

// 选择第二个变量
int selectSecondVariable(int i, int m) {
    int j = i;
    while (j == i) {
        j = rand() % m;
    }
    return j;
}

// 计算L和H的值
void computeLH(double alpha1, double alpha2, double y1, double y2, double *L, double *H) {
    if (y1 == y2) {
        *L = fmax(0, alpha1 + alpha2 - C);
        *H = fmin(C, alpha1 + alpha2);
    } else {
        *L = fmax(0, alpha2 - alpha1);
        *H = fmin(C, C + alpha2 - alpha1);
    }
}

// 计算新的拉格朗日乘子值
int updateAlpha(double *x1, double *x2, double *y, SupportVector *sv, int i, int j, double *E, int n) {
    double alpha1 = sv[i].alpha;
    double alpha2 = sv[j].alpha;
    double y1 = sv[i].y;
    double y2 = sv[j].y;
    double L, H;
    computeLH(alpha1, alpha2, y1, y2, &L, &H);
    if (L == H) {
        return 0;
    }
    double eta = dotProduct(x1, x1, n) + dotProduct(x2, x2, n) - 2 * dotProduct(x1, x2, n);
    if (eta <= 0) {
        return 0;
    }
    double alpha2_new = alpha2 + y2 * (E[i] - E[j]) / eta;
    if (alpha2_new < L) {
        alpha2_new = L;
    } else if (alpha2_new > H) {
        alpha2_new = H;
    }
    if (fabs(alpha2_new - alpha2) < TOL * (alpha2_new + alpha2 + TOL)) {
        return 0;
    }
    double alpha1_new = alpha1 + y1 * y2 * (alpha2 - alpha2_new);
    double b1 = -E[i] - y1 * dotProduct(x1, x1, n) * (alpha1_new - alpha1) - y2 * dotProduct(x2, x1, n) * (alpha2_new - alpha2) + sv[i].b;
    double b2 = -E[j] - y1 * dotProduct(x1, x2, n) * (alpha1_new - alpha1) - y2 * dotProduct(x2, x2, n) * (alpha2_new - alpha2) + sv[j].b;
    if (alpha1_new > 0 && alpha1_new < C) {
        sv[i].b = b1;
    } else if (alpha2_new > 0 && alpha2_new < C) {
        sv[j].b = b2;
    } else {
        sv[i].b = (b1 + b2) / 2;
        sv[j].b = sv[i].b;
    }
    sv[i].alpha = alpha1_new;
    sv[j].alpha = alpha2_new;
    return 1;
}

// SMO算法
void smoAlgorithm(double *X, double *Y, int m, int n, SupportVector *sv) {
    double *E = (double *) malloc(m * sizeof(double));
    for (int i = 0; i < m; i++) {
        E[i] = errorFunction(X + i * n, Y[i], sv, n);
    }
    int numChanged = 0;
    int examineAll = 1;
    int iter = 0;
    while (iter < MAX_ITER && (numChanged > 0 || examineAll)) {
        numChanged = 0;
        if (examineAll) {
            for (int i = 0; i < m; i++) {
                numChanged += examineExample(X, Y, m, n, sv, E, i);
            }
        } else {
            for (int i = 0; i < m; i++) {
                if (sv[i].alpha > 0 && sv[i].alpha < C) {
                    numChanged += examineExample(X, Y, m, n, sv, E, i);
                }
            }
        }
        if (examineAll) {
            examineAll = 0;
        } else if (numChanged == 0) {
            examineAll = 1;
        }
        iter++;
    }
    free(E);
}

// 根据训练数据得到支持向量
int getSupportVectors(double *X, double *Y, int m, int n, SupportVector *sv) {
    int numSV = 0;
    for (int i = 0; i < m; i++) {
        if (sv[i].alpha > 0) {
            numSV++;
        }
    }
    SupportVector *result = (SupportVector *) malloc(numSV * sizeof(SupportVector));
    int count = 0;
    for (int i = 0; i < m; i++) {
        if (sv[i].alpha > 0) {
            result[count].alpha = sv[i].alpha;
            result[count].b = sv[i].b;
            result[count].x = X + i * n;
            result[count].y = Y[i];
            count++;
        }
    }
    memcpy(sv, result, numSV * sizeof(SupportVector));
    free(result);
    return numSV;
}

int main() {
    // 读取训练数据
    int m, n;
    scanf("%d%d", &m, &n);
    double *X = (double *) malloc(m * n * sizeof(double));
    double *Y = (double *) malloc(m * sizeof(double));
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            scanf("%lf", &X[i * n + j]);
        }
        scanf("%lf", &Y[i]);
    }
    // 初始化支持向量
    SupportVector *sv = (SupportVector *) malloc(m * sizeof(SupportVector));
    for (int i = 0; i < m; i++) {
        sv[i].alpha = 0.0;
        sv[i].b = 0.0;
        sv[i].x = X + i * n;
        sv[i].y = Y[i];
    }
    // 训练模型
    smoAlgorithm(X, Y, m, n, sv);
    // 获取支持向量并输出
    int numSV = getSupportVectors(X, Y, m, n, sv);
    printf("Number of support vectors: %d\n", numSV);
    for (int i = 0; i < numSV; i++) {
        printf("Alpha[%d] = %lf, b = %lf, x = (", i, sv[i].alpha, sv[i].b);
        for (int j = 0; j < n; j++) {
            printf("%lf", sv[i].x[j]);
            if (j < n - 1) {
                printf(", ");
            }
        }
        printf("), y = %lf\n", sv[i].y);
    }
    // 释放内存
    free(X);
    free(Y);
    free(sv);
    return 0;
}

以上代码中，我们通过`selectSecondVariable`函数随机选择第二个变量，使用`computeLH`函数计算L和H的值，使用`updateAlpha`函数更新拉格朗日乘子，并使用`examineExample`函数选取第一个变量。在`smoAlgorithm`函数中，我们使用平均学习速率和软间隔方法来训练模型，使用`getSupportVectors`函数获取支持向量。在`main`函数中，我们读取训练数据，初始化支持向量，训练模型并输出支持向量的信息。