矩阵和特征向量的 QR 分解

Question

也许这更多是一个理论问题而不是编程问题，但我找不到可接受的答案。我用 C++ 编写了这段代码，它使用 HouseholderMain 函数对对称矩阵 A 进行三对角化（该函数会覆盖 A 并使用对 A 进行三对角化的矩阵更新 PP），然后使用函数 QRtridiag 求出 A 的 QR 分解（该函数会使用其上层覆盖 A ）三角形式并给出分解的矩阵 Q），最后用函数 AutovTriang 计算上三角矩阵 A 的特征值向量 a 和特征向量 X 的矩阵。

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;
        }
    }
}

我的问题很简单：如何将最终上三角矩阵的特征值和特征向量与原始矩阵的特征值和特征向量相关联？我试图理解 QR 算法，但我不确定它是否能满足我的需要。我已经有了三角矩阵的特征系统，并且我认为有一种简单的方法可以使用矩阵 PP 和 Q 来计算原始矩阵的特征系统，但我不知道该怎么做。 有人能帮我吗？非常感谢你。

Answer 1

你不断改变

A

，所以我们把原来的

A

称为A，Householder A'之后得到的

A

，以及QR分解后得到的

A

A''。

然后我们有：

A' = PP^T A PP

和：

A'' = Q^T A'

所以：

A''= Q^T PP^T A PP

或者走另一条路，这就是您正在寻找的转变：

A = PP Q A'' PP^T