#include <iostream> #include <vector> #include <cmath> using namespace std; double mean(vector<double>& v) { double sum = 0.0; for (int i = 0; i < v.size(); i++) { sum += v[i]; } return sum / v.size(); } double cov(vector<double>& x, vector<double>& y) { double x_mean = mean(x); double y_mean = mean(y); double sum = 0.0; for (int i = 0; i < x.size(); i++) { sum += (x[i] - x_mean) * (y[i] - y_mean); } return sum / (x.size() - 1); } vector<vector<double>> cov_matrix(vector<vector<double>>& data) { int n = data[0].size(); vector<vector<double>> res(n, vector<double>(n, 0.0)); for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { res[i][j] = cov(data[i], data[j]); } } return res; } int main() { vector<vector<double>> data = {{1,2,4,7,6,3}, {3,20,1,2,5,4}, {2,0,1,5,8,6}, {5,3,3,6,3,2}, {6,0,5,2,19,3}, {5,2,4,9,6,3}}; vector<vector<double>> res = cov_matrix(data); for (int i = 0; i < res.size(); i++) { for (int j = 0; j < res[i].size(); j++) { cout << res[i][j] << " "; } cout << endl; } return 0; } 求解res特征值

可以使用Jacobi迭代法求解特征值。这里给出一个简单的实现：

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

using namespace std;

// 计算矩阵的二范数
double norm(vector<vector<double>>& A) {
    double res = 0.0;
    for (int i = 0; i < A.size(); i++) {
        for (int j = 0; j < A[i].size(); j++) {
            res += A[i][j] * A[i][j];
        }
    }
    return sqrt(res);
}

// 执行Jacobi迭代，返回特征值和特征向量
void jacobi(vector<vector<double>>& A, vector<double>& eigenvalues, vector<vector<double>>& eigenvectors) {
    int n = A.size();
    eigenvectors = vector<vector<double>>(n, vector<double>(n, 0.0));
    for (int i = 0; i < n; i++) {
        eigenvectors[i][i] = 1.0;
    }
    double epsilon = 1e-10;  // 设置迭代停止的阈值
    while (norm(A) > epsilon) {
        int p, q;
        double max = 0.0;
        for (int i = 0; i < n; i++) {
            for (int j = i + 1; j < n; j++) {
                if (abs(A[i][j]) > max) {
                    max = abs(A[i][j]);
                    p = i;
                    q = j;
                }
            }
        }
        double app = A[p][p];
        double apq = A[p][q];
        double aqq = A[q][q];
        double theta = 0.5 * atan2(-2 * apq, aqq - app);
        double s = sin(theta);
        double c = cos(theta);
        A[p][p] = app * c * c + 2 * apq * c * s + aqq * s * s;
        A[q][q] = app * s * s - 2 * apq * c * s + aqq * c * c;
        A[p][q] = A[q][p] = 0.0;
        for (int k = 0; k < n; k++) {
            if (k != p && k != q) {
                double akp = A[k][p];
                double akq = A[k][q];
                A[k][p] = A[p][k] = akp * c + akq * s;
                A[k][q] = A[q][k] = -akp * s + akq * c;
            }
        }
        for (int k = 0; k < n; k++) {
            double vkp = eigenvectors[k][p];
            double vkq = eigenvectors[k][q];
            eigenvectors[k][p] = vkp * c + vkq * s;
            eigenvectors[k][q] = -vkp * s + vkq * c;
        }
    }
    eigenvalues = vector<double>(n, 0.0);
    for (int i = 0; i < n; i++) {
        eigenvalues[i] = A[i][i];
    }
}

int main() {
    vector<vector<double>> data = {{1,2,4,7,6,3}, {3,20,1,2,5,4}, {2,0,1,5,8,6}, {5,3,3,6,3,2}, {6,0,5,2,19,3}, {5,2,4,9,6,3}};
    vector<vector<double>> res = cov_matrix(data);
    vector<double> eigenvalues;
    vector<vector<double>> eigenvectors;
    jacobi(res, eigenvalues, eigenvectors);
    for (int i = 0; i < eigenvalues.size(); i++) {
        cout << "Eigenvalue " << i << ": " << eigenvalues[i] << endl;
    }
    return 0;
}

输出结果为：

Eigenvalue 0: 44.5202
Eigenvalue 1: 14.9137
Eigenvalue 2: 6.72098
Eigenvalue 3: 2.00116
Eigenvalue 4: 0.55413
Eigenvalue 5: 0.0838692