如何通过PP与Q矩阵关联三角矩阵与原矩阵的eigenvalue及eigenvector？

代码实现 #

int main()
{
    int n = 4;
    double** A = createMatrix(n, n);
    double** PP = createMatrix(n, n);
    double** Q = createMatrix(n, n);
    double** X = createMatrix(n, n);
    double* a = createVector(n);

    A[0][0] = 1;
    A[0][1] = A[1][0] = 2;
    A[0][2] = A[2][0] = -3;
    A[0][3] = A[3][0] = 4;
    A[1][1] = 1;
    A[1][2] = A[2][1] = -1;
    A[1][3] = A[3][1] = 9;
    A[2][2] = 2;
    A[2][3] = A[3][2] = 10;
    A[3][3] = 0;

    HouseholderMain(A, PP, n);
    QRtridiag(A, Q, n);
    AutovTriang(A, X, a, n);

        return 0;
}

void HouseholderMain(double** A, double** PP, int n)
{
    int i, j, k, l;
    double sum;
    double** P = createMatrix(n, n);
    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            PP[i][j] = (i == j ? 1 : 0);
        }
    }
    for (int l = 0; l < n - 2; l++) {
        Householder2Cycle(A, P, n, l);
        double** PA = multMatrices(P, A, n);
        for (i = 0; i < n; i++) {
            for (j = 0; j < n; j++) {
                sum = 0;
                for (k = 0; k < n; k++) {
                    sum += PA[i][k] * P[k][j];
                }
                A[i][j] = sum;
            }
        }
        double** a = multMatrices(PP, P, n);
        for (i = 0; i < n; i++) {
            for (j = 0; j < n; j++) {
                PP[i][j] = a[i][j];
            }
        }
    }
}


void Householder2Cycle(double** A, double** P, int n, int k)
{
    int i, j;
    double sum, mod, mod2x, h;
    double* x = createVector(n);
    for (i = 0; i < k + 1; i++) {
        x[i] = 0;
    }
    for (i = k + 2; i < n; i++) {
        x[i] = A[i][k];
    }
    sum = 0;
    for (i = k + 1; i < n; i++) {
        sum += A[i][k] * A[i][k];
    }
    mod = (A[k+1][k] > 0 ? sqrt(sum) : -sqrt(sum));
    x[k + 1] = A[k + 1][k] + mod;
    sum = 0;
    for (i = 0; i < n; i++) {
        sum += x[i] * x[i];
    }
    mod2x = sum;
    h = mod2x / 2;
    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            P[i][j] = (i == j ? 1 : 0) - x[i] * x[j] / h;
        }
    }
    freeVector(x);
}


void QRtridiag(double** A, double** Q, int n)
{
    int i, j;
    double theta, c, s, t, p, pp;
    for (j = 0; j <= n - 2; j++) {
        for (i = j + 1; i < n; i++) {
            Q[i][j] = Q[j][i] = 0;
        }
    }
    for (i = 0; i < n; i++) {
        Q[i][i] = 1;
    }
    for (i = 0; i < n - 1; i++) {
        theta = A[i][i] / A[i+1][i];
        t = 1 / (fabs(theta) + sqrt(theta * theta + 1));
        if (theta < 0) t = -t;
        c = (1 - t * t) / (1 + t * t);
        s = 2 * t / (1 + t * t);
        p = A[i][i + 1];
        A[i][i] = c * A[i][i] + s * A[i + 1][i];
        A[i + 1][i] = 0;
        if (i < n - 2) {
            A[i][i + 2] = s * A[i + 1][i + 2];
            A[i + 1][i + 2] = c * A[i + 1][i + 2];
        }
        A[i][i + 1] = c * p + s * A[i + 1][i + 1];
        A[i + 1][i + 1] = -s * p + c * A[i + 1][i + 1];
        for (j = 0; j <= i; j++) {
            pp = Q[j][i];
            Q[j][i] = c * pp;
            Q[j][i+1] = -s * pp;
        }
        Q[i+1][i] = s;
        Q[i + 1][i + 1] = c;
    }
}


void AutovTriang(double** A, double** X, double* a, int n)
{
    int i, j, k;
    for (i = 0; i < n; i++) {
        a[i] = A[i][i];
    }
    for (k = n - 1; k >= 2; k--) {
        X[k][k] = 1;
        X[k - 1][k] = -A[k - 1][k] / (a[k - 1] - a[k]);
        for (i = k - 2; i >= 0; i--) {
            X[i][k] = (-A[i][i + 1] * X[i + 1][k] - A[i][i + 2] * X[i + 2][k]) / (a[i] - a[k]);
        }
    }
    X[1][1] = X[0][0] = 1;
    X[0][1] = -A[0][1] / (a[0] - a[1]);
    for (j = 0; j <= n - 2; j++) {
        for (i = j + 1; i < n; i++) {
            X[i][j] = 0;
        }
    }
}

问题背景 #

代码执行流程如下：

• 调用HouseholderMain将对称矩阵A三对角化，覆盖原A矩阵，同时累积正交变换矩阵PP；
• 调用QRtridiag对三对角化后的A执行QR分解，将A覆盖为上三角矩阵，同时得到正交变换矩阵Q；
• 调用AutovTriang计算上三角矩阵A的特征值向量a和特征向量矩阵X。

需要解决的问题：如何将最终得到的上三角矩阵的特征系统关联到原矩阵的特征系统？

解决方案 #

1. 特征值的关联 #

整个过程中所有变换都是正交相似变换（Householder矩阵、Givens矩阵都是正交矩阵，相似变换保持特征值不变），因此AutovTriang得到的特征值向量a直接就是原矩阵的特征值，无需额外转换。

2. 特征向量的关联 #

特征向量需要通过累积的正交变换矩阵逆推回原矩阵空间，步骤如下：

步骤1：从三角矩阵特征向量得到三对角矩阵的特征向量

QRtridiag中的Q是正交矩阵，它将三对角矩阵T相似变换为上三角矩阵R，即：
$$R = Q^T T Q$$
已知X是R的特征向量矩阵（满足$R X = X \text{diag}(a)$），则三对角矩阵T的特征向量矩阵为：
$$X_T = Q X$$
验证：$T X_T = T Q X = Q R X = Q X \text{diag}(a) = X_T \text{diag}(a)$，符合特征向量定义。

步骤2：从三对角矩阵特征向量得到原矩阵的特征向量

HouseholderMain中的PP是所有Householder变换矩阵的乘积，它将原矩阵A₀相似变换为三对角矩阵T，即：
$$T = PP A₀ PP^T$$
由于PP是正交矩阵，其逆矩阵等于转置矩阵PP^T，因此原矩阵A₀的特征向量矩阵为：
$$X_{\text{original}} = PP^T X_T = PP^T Q X$$
验证：$A₀ X_{\text{original}} = A₀ PP^T Q X = PP^T T Q X = PP^T Q X \text{diag}(a) = X_{\text{original}} \text{diag}(a)$，符合特征向量定义。

3. 代码实现建议 #

在现有代码基础上，添加矩阵转置和乘法逻辑即可得到原矩阵的特征向量：

• 先计算PP的转置矩阵PP_T；
• 计算Q_X = multMatrices(Q, X, n)；
• 最终原矩阵特征向量矩阵X_original = multMatrices(PP_T, Q_X, n)；
• 注意：如果需要特征向量是单位向量，可能需要对X_original的每一列做归一化处理（因为AutovTriang得到的X不一定是单位向量）。

内容的提问来源于stack exchange，提问作者D. Alfano