[函数] Simple Linear Regression Cholesky

使用Cholesky分解求解一元线性回归。

Input

• 自变量序列 $\boldsymbol{X}$
• 因变量序列 $\boldsymbol{y}$

Execution

1. 设计矩阵，第一列为常数项，第二列为自变量的值 $X = [1, \boldsymbol{X}]$

2. 构造正规方程

   $\text{XtX} = X^T * X$

   和

   $\text{Xty} = X^T * y$

   则待求系数 $θ$ 满足： $\text{XtX} * θ = \text{Xty}$

3. 对 $\text{XtX}$ 做 Cholesky 分解： $XtX = L * L^T$ , 其中 $L$ 为下三角矩阵

4. 依次解方程： $L * z = \text{Xty}$

   和

   $L^T * θ = z$

5. 得到回归系数 $α = θ[0]$ # 截距 $β = θ[1]$ # 斜率

Demo Code (C++)

#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>

// Cholesky分解函数
bool cholesky_decomposition(const std::vector<std::vector<double>>& A,
    std::vector<std::vector<double>>& L) {
    int n = A.size();
    L = std::vector<std::vector<double>>(n, std::vector<double>(n, 0.0));

    for (int i = 0; i < n; i++) {
        for (int j = 0; j <= i; j++) {
            double sum = 0.0;

            if (j == i) {
                for (int k = 0; k < j; k++) {
                    sum += L[j][k] * L[j][k];
                }
                double diag = A[j][j] - sum;
                if (diag <= 0) {
                    return false; // 矩阵不是正定的
                }
                L[j][j] = sqrt(diag);
            }
            else {
                for (int k = 0; k < j; k++) {
                    sum += L[i][k] * L[j][k];
                }
                L[i][j] = (A[i][j] - sum) / L[j][j];
            }
        }
    }
    return true;
}

// 解下三角方程组 L * y = b
void solve_lower_triangular(const std::vector<std::vector<double>>& L,
    const std::vector<double>& b,
    std::vector<double>& y) {
    int n = L.size();
    y.resize(n);

    for (int i = 0; i < n; i++) {
        double sum = 0.0;
        for (int j = 0; j < i; j++) {
            sum += L[i][j] * y[j];
        }
        y[i] = (b[i] - sum) / L[i][i];
    }
}

// 解上三角方程组 L^T * x = y
void solve_upper_triangular(const std::vector<std::vector<double>>& L,
    const std::vector<double>& y,
    std::vector<double>& x) {
    int n = L.size();
    x.resize(n);

    for (int i = n - 1; i >= 0; i--) {
        double sum = 0.0;
        for (int j = i + 1; j < n; j++) {
            sum += L[j][i] * x[j]; // 注意：L^T[j][i] = L[i][j]
        }
        x[i] = (y[i] - sum) / L[i][i];
    }
}

// 使用Cholesky分解的线性回归
bool simple_linear_regression_cholesky(const std::vector<double>& X,
    const std::vector<double>& y,
    double& alpha, double& beta) {
    int m = X.size();
    if (m <= 1) return false;

    // 构建正规方程: X^T * X * beta = X^T * y
    // 其中 X 是设计矩阵，第一列为1，第二列为X值

    // 计算 X^T * X
    double sum_x = 0.0, sum_x2 = 0.0;
    double sum_y = 0.0, sum_xy = 0.0;

    for (int i = 0; i < m; i++) {
        sum_x += X[i];
        sum_x2 += X[i] * X[i];
        sum_y += y[i];
        sum_xy += X[i] * y[i];
    }

    // 正规方程矩阵 A = X^T * X
    std::vector<std::vector<double>> A = {
        {static_cast<double>(m), sum_x},
        {sum_x, sum_x2}
    };

    // 右侧向量 b = X^T * y
    std::vector<double> b = { sum_y, sum_xy };

    // Cholesky分解 A = L * L^T
    std::vector<std::vector<double>> L;
    if (!cholesky_decomposition(A, L)) {
        return false;
    }

    // 解方程组: L * L^T * beta = b
    // 第一步: L * y = b
    std::vector<double> y_temp;
    solve_lower_triangular(L, b, y_temp);

    // 第二步: L^T * beta = y
    std::vector<double> beta_vec;
    solve_upper_triangular(L, y_temp, beta_vec);

    // 提取系数
    alpha = beta_vec[0]; // 截距
    beta = beta_vec[1];  // 斜率

    return true;
}

int main() {
    // 示例数据
    std::vector<double> T = { 1.0, 2.0, 3.0, 4.0 };
    std::vector<double> log_ATMIV = { 0.1, 0.2, 0.3, 0.4 };

    double alpha, beta;
    if (simple_linear_regression_cholesky(T, log_ATMIV, alpha, beta)) {
        std::cout << "回归系数:" << std::endl;
        std::cout << "alpha (截距) = " << alpha << std::endl;
        std::cout << "beta (斜率) = " << beta << std::endl;
    }
    else {
        std::cerr << "线性回归失败" << std::endl;
    }

    return 0;
}