linear regression in C

linear_regression.c

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

#define NOISE_SCALE 5.0
#define EPSILON 1e-12
#define MAX_FILE_PATH_LENGTH 100
#define MAX_FEATURE_DIMENSION 2
#define MAX_DATASET_SIZE 10000
#define TEST_DATASET_SIZE_RATIO 0.1

// 데이터 타입
typedef double Feature[MAX_FEATURE_DIMENSION];
typedef double Target;
typedef struct {
    Feature feature;
    Target target;
} Data;
typedef struct {
    Data data[MAX_DATASET_SIZE];
    int size;
    int feature_dimension;
} Dataset;

// 파라미터 타입
typedef double Weight[MAX_FEATURE_DIMENSION];
typedef double Bias;
typedef struct {
    Weight weight;
    Bias bias;
    int feature_dimension;
} Parameter;

// 유틸리티
double get_mean(const double array[], const int size) {
    double mean = 0.0;
    for(int i = 0; i < size; i++) {
        mean += array[i] / size;
    }
    return mean;
}
double get_std(const double array[], const int size) {
    double mean = get_mean(array, size);
    double std = 0.0;
    for(int i = 0; i < size; i++) {
        std += (array[i] - mean) * (array[i] - mean) / (size - 1);
    }
    return sqrt(std);
}
double get_random(void) {
    return (double)rand() / RAND_MAX;
}
double get_noise(const double scale) {
    return get_random() * scale - scale / 2;
}
double oracle(double x1, double x2) {
    return 3 * x1 * x1 - 5 * x2 * x2 + 1;
}

// 파라미터 초기화
void initialize_parameter(Parameter *parameter, const int feature_dimension) {
    parameter->feature_dimension = feature_dimension;
    for(int j = 0; j < feature_dimension; j++) {
        parameter->weight[j] = 0.0;
    }
    parameter->bias = 0.0;
}

// 모델 1 정의
double model1(const Parameter *parameter, const Feature *feature) {
    double prediction = 0.0;
    for(int j = 0; j < parameter->feature_dimension; j++) {
        prediction += parameter->weight[j] * (*feature)[j];
    }
    prediction += parameter->bias;
    return prediction;
}
void optimizer1(
    const Parameter *parameter,
    const Data *data,
    const double error,
    const int dataset_size,
    const double regularization_strength,
    Parameter *gradient
) {
    for(int j = 0; j < parameter->feature_dimension; j++) {
        gradient->weight[j] += 2 * error * data->feature[j] / dataset_size;
        gradient->weight[j] += 2 * regularization_strength * parameter->weight[j] / dataset_size;
    }
    gradient->bias += 2 * error / dataset_size;
    gradient->bias += 2 * regularization_strength * parameter->bias / dataset_size;
}

// 모델 2 정의
double model2(const Parameter *parameter, const Feature *feature) {
    double prediction = 0.0;
    for(int j = 0; j < parameter->feature_dimension; j++) {
        prediction += parameter->weight[j] * (*feature)[j] * (*feature)[j];
    }
    prediction += parameter->bias;
    return prediction;
}
void optimizer2(
    const Parameter *parameter,
    const Data *data,
    const double error,
    const int dataset_size,
    const double regularization_strength,
    Parameter *gradient
) {
    for(int j = 0; j < parameter->feature_dimension; j++) {
        gradient->weight[j] += 2 * error * data->feature[j] * data->feature[j] / dataset_size;
        gradient->weight[j] += 2 * regularization_strength * parameter->weight[j] / dataset_size;
    }
    gradient->bias += 2 * error / dataset_size;
    gradient->bias += 2 * regularization_strength * parameter->bias / dataset_size;
}

// 학습기
Parameter train(
    double (*model)(const Parameter *, const Feature *),
    void (*optimizer)(
        const Parameter *,
        const Data *,
        const double,
        const int,
        const double,
        Parameter *
    ),
    const Dataset *dataset_train,
    const Dataset *dataset_validation,
    const int number_of_epochs,
    const int print_interval,
    const double learning_rate,
    const double regularization_strength
) {

    // 파라미터 초기화
    Parameter best_parameter;
    Parameter current_parameter;
    initialize_parameter(&current_parameter, dataset_train->feature_dimension);

    // 학습 변수 초기화
    int initial_flag = 0;
    double best_loss_validation;
    double previous_loss_train;

    for(int epoch = 1; epoch <= number_of_epochs; epoch++) {

        // 그래디언트 초기화
        Parameter gradient;
        initialize_parameter(&gradient, dataset_train->feature_dimension);

        // 손실 계산 (학습 데이터셋)
        double loss_train = 0.0;
        for(int i = 0; i < dataset_train->size; i++) {
            double prediction = model(&current_parameter, &dataset_train->data[i].feature);
            double error = prediction - dataset_train->data[i].target;
            loss_train += error * error / dataset_train->size;
            optimizer(&current_parameter, &dataset_train->data[i], error, dataset_train->size, regularization_strength, &gradient);
        }

        // 파라미터 업데이트
        for(int j = 0; j < current_parameter.feature_dimension; j++) {
            current_parameter.weight[j] -= gradient.weight[j] * learning_rate;
        }
        current_parameter.bias -= gradient.bias * learning_rate;

        // 손실 계산 (검증 데이터셋)
        double loss_validation = 0.0;
        for(int i = 0; i < dataset_validation->size; i++) {
            double prediction = model(&current_parameter, &dataset_validation->data[i].feature);
            double error = prediction - dataset_validation->data[i].target;
            loss_validation += error * error / dataset_validation->size;
        }

        // 에폭 끝
        if(initial_flag == 0) {
            best_loss_validation = loss_validation;
            best_parameter = current_parameter;
            previous_loss_train = loss_train;
            initial_flag = 1;
        } else {
            if(loss_validation < best_loss_validation) {
                best_loss_validation = loss_validation;
                best_parameter = current_parameter;
            }
            if(fabs(loss_train - previous_loss_train) < EPSILON) {
                printf("Early stopping at epoch %d\n", epoch);
                break;
            }
            previous_loss_train = loss_train;
        }

        // 출력
        if(epoch % print_interval == 0) {
            printf("[Epoch %d] train: %lf, validation: %lf\n", epoch, loss_train, loss_validation);
        }
    }
    return best_parameter;
}

int main(int argc, char *argv[]) {

    // 시드값 고정
    srand(42);

    // 학습 데이터셋 생성
    printf("Generating train dataset\n");
    FILE *fp_train_generator = fopen("train.data", "w");
    for(int i = 0; i < MAX_DATASET_SIZE; i++) {
        double x1 = get_random() * 10 - 5;
        double x2 = get_random() * 10 - 5;
        double y = oracle(x1, x2) + get_noise(NOISE_SCALE);
        fprintf(fp_train_generator, "%lf %lf %lf\n", x1, x2, y);
    }
    fclose(fp_train_generator);

    // 테스트 데이터셋 생성
    printf("Generating test dataset\n");
    FILE *fp_test_generator = fopen("test.data", "w");
    for(int i = 0; i < (int)(MAX_DATASET_SIZE * TEST_DATASET_SIZE_RATIO); i++) {
        double x1 = get_random() * 10 - 5;
        double x2 = get_random() * 10 - 5;
        double y = oracle(x1, x2);
        fprintf(fp_test_generator, "%lf %lf %lf\n", x1, x2, y);
    }
    fclose(fp_test_generator);

    printf("Done!\n");

    printf("\n");
    printf("----------------------\n");
    printf("\n");

    // 설정값 세팅
    char train_data_path[MAX_FILE_PATH_LENGTH + 1];
    char test_data_path[MAX_FILE_PATH_LENGTH + 1];
    int feature_dimension;
    double train_validation_ratio;
    int number_of_epochs;
    int warm_up_epochs;
    int print_interval;
    double learning_rate;
    double regularization_strength;
    sscanf(argv[1], "%s", train_data_path);  // train.data
    sscanf(argv[2], "%s", test_data_path);  // test.data
    sscanf(argv[3], "%d", &feature_dimension);  // 2
    sscanf(argv[4], "%lf", &train_validation_ratio);  // 0.9
    sscanf(argv[5], "%d", &number_of_epochs);  // 10000
    sscanf(argv[6], "%d", &print_interval);  // 1000
    sscanf(argv[7], "%lf", &learning_rate);  // 0.001
    sscanf(argv[8], "%lf", &regularization_strength);  // 0.0
    printf("Config:\n");
    printf("- train_data_path: %s\n", train_data_path);
    printf("- test_data_path: %s\n", test_data_path);
    printf("- feature_dimension: %d\n", feature_dimension);
    printf("- train_validation_ratio: %.1lf\n", train_validation_ratio);
    printf("- number_of_epochs: %d\n", number_of_epochs);
    printf("- print_interval: %d\n", print_interval);
    printf("- learning_rate: %.04lf\n", learning_rate);
    printf("- regularization_strength: %.04lf\n", regularization_strength);

    printf("\n");
    printf("----------------------\n");
    printf("\n");

    // 학습 데이터셋 로드
    printf("Loading train dataset\n");
    Dataset dataset_train;
    dataset_train.feature_dimension = feature_dimension;
    FILE *fp_train = fopen(train_data_path, "r");
    int i_train = 0;
    while(fscanf(
        fp_train,
        "%lf %lf %lf",
        &dataset_train.data[i_train].feature[0],
        &dataset_train.data[i_train].feature[1],
        &dataset_train.data[i_train].target
    ) != EOF) {
        i_train++;
    }
    dataset_train.size = i_train;
    fclose(fp_train);

    // 검증 데이터셋 분리
    printf("Splitting validation dataset\n");
    Dataset dataset_validation;
    dataset_validation.feature_dimension = feature_dimension;
    int size_train = (int)(dataset_train.size * train_validation_ratio);
    dataset_train.size = size_train;
    for(int i = size_train; i < i_train; i++) {
        dataset_validation.data[i - size_train] = dataset_train.data[i];
    }
    dataset_validation.size = i_train - size_train;

    // 학습 데이터셋 정규화
    Feature feature_mean;
    Feature feature_std;
    Target target_mean;
    Target target_std;
    printf("Standardizing train dataset\n");
    for(int j = 0; j < feature_dimension; j++) {
        double array[MAX_DATASET_SIZE];
        for(int i = 0; i < dataset_train.size; i++) {
            array[i] = dataset_train.data[i].feature[j];
        }
        feature_mean[j] = get_mean(array, dataset_train.size);
        feature_std[j] = get_std(array, dataset_train.size);
    }
    double array[MAX_DATASET_SIZE];
    for(int i = 0; i < dataset_train.size; i++) {
        array[i] = dataset_train.data[i].target;
    }
    target_mean = get_mean(array, dataset_train.size);
    target_std = get_std(array, dataset_train.size);
    for(int j = 0; j < feature_dimension; j++) {
        for(int i = 0; i < dataset_train.size; i++) {
            dataset_train.data[i].feature[j] = (dataset_train.data[i].feature[j] - feature_mean[j]) / feature_std[j];
        }
    }
    for(int i = 0; i < dataset_train.size; i++) {
        dataset_train.data[i].target = (dataset_train.data[i].target - target_mean) / target_std;
    }

    // 검증 데이터셋 정규화
    printf("Standardizing validation dataset\n");
    for(int j = 0; j < feature_dimension; j++) {
        for(int i = 0; i < dataset_validation.size; i++) {
            dataset_validation.data[i].feature[j] = (dataset_validation.data[i].feature[j] - feature_mean[j]) / feature_std[j];
        }
    }
    for(int i = 0; i < dataset_validation.size; i++) {
        dataset_validation.data[i].target = (dataset_validation.data[i].target - target_mean) / target_std;
    }

    // 테스트 데이터셋 로드
    printf("Loading test dataset\n");
    Dataset dataset_test;
    dataset_test.feature_dimension = feature_dimension;
    FILE *fp_test = fopen(test_data_path, "r");
    int i_test = 0;
    while(fscanf(
        fp_test,
        "%lf %lf %lf",
        &dataset_test.data[i_test].feature[0],
        &dataset_test.data[i_test].feature[1],
        &dataset_test.data[i_test].target
    ) != EOF) {
        i_test++;
    }
    dataset_test.size = i_test;
    fclose(fp_test);

    // 테스트 데이터셋 정규화
    printf("Standardizing test dataset\n");
    for(int j = 0; j < feature_dimension; j++) {
        for(int i = 0; i < dataset_test.size; i++) {
            dataset_test.data[i].feature[j] = (dataset_test.data[i].feature[j] - feature_mean[j]) / feature_std[j];
        }
    }
    for(int i = 0; i < dataset_test.size; i++) {
        dataset_test.data[i].target = (dataset_test.data[i].target - target_mean) / target_std;
    }

    printf("Done!\n");

    printf("\n");
    printf("----------------------\n");
    printf("\n");

    // 모델 1 학습
    printf("Training model 1\n");
    Parameter parameter1 = train(
        model1,
        optimizer1,
        &dataset_train,
        &dataset_validation,
        number_of_epochs,
        print_interval,
        learning_rate,
        regularization_strength
    );
    printf("Done!\n");

    printf("\n");
    printf("----------------------\n");
    printf("\n");

    // 모델 1 테스트
    printf("Testing model 1:\n");
    double loss_test_model1 = 0.0;
    for(int i = 0; i < dataset_test.size; i++) {
        double prediction = model1(&parameter1, &dataset_test.data[i].feature);
        double error = prediction - dataset_test.data[i].target;
        loss_test_model1 += error * error / dataset_test.size;
    }
    printf("- MSE: %lf\n", loss_test_model1);

    printf("\n");
    printf("----------------------\n");
    printf("\n");

    // 모델 2 학습
    printf("Training model 2\n");
    Parameter parameter2 = train(
        model2,
        optimizer2,
        &dataset_train,
        &dataset_validation,
        number_of_epochs,
        print_interval,
        learning_rate,
        regularization_strength
    );
    printf("Done!\n");

    printf("\n");
    printf("----------------------\n");
    printf("\n");

    // 모델 2 테스트
    printf("Testing model 2:\n");
    double loss_test_model2 = 0.0;
    for(int i = 0; i < dataset_test.size; i++) {
        double prediction = model2(&parameter2, &dataset_test.data[i].feature);
        double error = prediction - dataset_test.data[i].target;
        loss_test_model2 += error * error / dataset_test.size;
    }
    printf("- MSE: %lf\n", loss_test_model2);

    printf("\n");
    printf("----------------------\n");
    printf("\n");

    return 0;
}