カーネル主成分分析のC言語実装

参考文献

機械学習のためのカーネル100問 with Python
著者 鈴木 讓
発売日 2021/12/28

準備

オンラインコンパイラを使用します。

ソースコード

sample.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>  // Add this line to include the time function

#define N 100

typedef struct {
    double x, y;
} Point;

Point points[N];

double kernel(Point a, Point b, double sigma) {
    double diff = a.x - b.x;
    return exp(-diff * diff / (2 * sigma * sigma));
}

void compute_kernel_matrix(double sigma, double kernel_matrix[N][N]) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            kernel_matrix[i][j] = kernel(points[i], points[j], sigma);
        }
    }
}

void jacobi_eigenvalue_algorithm(double A[N][N], double eigenvalues[N], double eigenvectors[N][N], int max_iterations, double tolerance) {
    int i, j, k, l, m, iter;
    double Amax, theta, t, c, s, tau;
    double tmp, tmp2, tmp3;

    // Initialize eigenvectors matrix to identity matrix
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            eigenvectors[i][j] = (i == j) ? 1.0 : 0.0;
        }
    }

    for (iter = 0; iter < max_iterations; iter++) {
        // Find the largest off-diagonal element
        Amax = 0.0;
        for (i = 0; i < N - 1; i++) {
            for (j = i + 1; j < N; j++) {
                if (fabs(A[i][j]) > fabs(Amax)) {
                    Amax = A[i][j];
                    k = i;
                    l = j;
                }
            }
        }

        // If the largest off-diagonal element is below tolerance, we are done
        if (fabs(Amax) < tolerance) {
            break;
        }

        // Calculate Jacobi rotation
        theta = (A[l][l] - A[k][k]) / (2.0 * A[k][l]);
        t = ((theta >= 0) ? 1.0 : -1.0) / (fabs(theta) + sqrt(1.0 + theta * theta));
        c = 1.0 / sqrt(1.0 + t * t);
        s = t * c;
        tau = s / (1.0 + c);

        // Update the matrix A
        tmp = A[k][l];
        A[k][l] = 0.0;
        A[k][k] -= t * tmp;
        A[l][l] += t * tmp;

        for (i = 0; i < k; i++) {
            tmp = A[i][k];
            A[i][k] = tmp - s * (A[i][l] + tau * tmp);
            A[i][l] = A[i][l] + s * (tmp - tau * A[i][l]);
        }

        for (i = k + 1; i < l; i++) {
            tmp = A[k][i];
            A[k][i] = tmp - s * (A[i][l] + tau * A[k][i]);
            A[i][l] = A[i][l] + s * (tmp - tau * A[i][l]);
        }

        for (i = l + 1; i < N; i++) {
            tmp = A[k][i];
            A[k][i] = tmp - s * (A[l][i] + tau * A[k][i]);
            A[l][i] = A[l][i] + s * (tmp - tau * A[l][i]);
        }

        // Update eigenvectors
        for (i = 0; i < N; i++) {
            tmp2 = eigenvectors[i][k];
            tmp3 = eigenvectors[i][l];
            eigenvectors[i][k] = tmp2 - s * (tmp3 + tau * tmp2);
            eigenvectors[i][l] = tmp3 + s * (tmp2 - tau * tmp3);
        }
    }

    // Copy eigenvalues
    for (i = 0; i < N; i++) {
        eigenvalues[i] = A[i][i];
    }
}

void kernel_pca(int num_components, double sigma) {
    double kernel_matrix[N][N];
    double eigenvalues[N];
    double eigenvectors[N][N];
    int max_iterations = 1000;
    double tolerance = 1e-10;

    compute_kernel_matrix(sigma, kernel_matrix);

    jacobi_eigenvalue_algorithm(kernel_matrix, eigenvalues, eigenvectors, max_iterations, tolerance);

    for (int i = 0; i < num_components; i++) {
        printf("Eigenvalue: %f\n", eigenvalues[i]);
        printf("Eigenvector: ");
        for (int j = 0; j < N; j++) {
            printf("%f ", eigenvectors[j][i]);
        }
        printf("\n");
    }
}

int main() {
    srand(time(0));
    for (int i = 0; i < N; i++) {
        points[i].x = (double)rand() / RAND_MAX;
        points[i].y = (double)rand() / RAND_MAX;
    }
    kernel_pca(1, 1.0);
    return 0;
}