QR decomposition of a matrix and eigenvectors

Question

perhaps this is more a question of theory than programming, but I can't find an acceptable answer. I wrote this code in C++, which tridiagonalizes the symmetric matrix A with the function HouseholderMain (which overwrites A and updates PP with the matrix that tridiagonalizes A), then finds the QR decomposition of A with the function QRtridiag (which overwrites A with its upper triangular form and gives the matrix Q of the decomposition) and finally calculates the vector of eigenvalues a and the matrix of eigenvectors X of the upper triangular matrix A with the function AutovTriang.

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

My question is simply the following: how can I relate eigenvalues and eigenvectors of the final upper triangular matrix to those of the original matrix? I tried to understand the QR algorithm but I am not sure it offers what I need. I have the eigensystem of the triangular matrix yet and I think there is an easy way to calculate the eigensystem of the original matrix using the matrices PP and Q, but I cannot figure out how to do it. Can someone help me? Thanks you so much.

Answer 1

A direct relation is that the product of the diagonal elements of R (which are the eignvalues of R) times the determinant of Q (which is +1 or -1), from the QR decomposition, is equal to the product of the eigenvalues of your original A.

There are iterative procedures to get the eigenvalues of QR with repeated applications of the QR decomposition. Those procedures are sped up by the tri-diagonalization you have already done. The eigenvalues of your original A are equal to the eigenvalues of its tridiagonalization, and you can back-convert the eigenvectors using the transformation matrix you accumulated in PP. If we call A' the tridiagonal form of A, then A' = PP^T A PP, so A = PP A' PP^T, where that transformation would convert the eigenvectors of A' to A.

The simplest procedure to get the eigenvalues of A' is to get its QR decomposition, and compute RQ. Do a QR decomposition on that for a new Q and R, compute RQ again, until the diagonals of successive RQ's have sufficiently converged. The sequence of RQ's will approach diagonal, with the eigenvalues of A or A' on the diagonal.

This can take many iterations. A way to significantly speed up the convergence is by subtracting a constant from the diagonal elements of QR, and then adding that constant back to RQ. Choosing the minimum absolute value of the diagonal elements of QR as the constant works well. In a test I tried, that reduced the iterations from 547 to 32.

Answer 2

I write the full code with little improvements in the following. Any suggestion would be appreciated. Thank you.

int main() 
{
    int n = 4;
    double** A = createMatrix(n, n);
    A[0][0] = 3.69583037;
    A[0][1] = A[1][0] = 0.57095943;
    A[0][2] = A[2][0] = 1.25982728;
    A[0][3] = A[3][0] = -0.91138766;
    A[1][1] = -5.72767862;
    A[1][2] = A[2][1] = 1.40997728;
    A[1][3] = A[3][1] = 3.09590731;
    A[2][2] = 0.84352098;
    A[2][3] = A[3][2] = -0.00810195;
    A[3][3] = 2.71813063;

    double** eigenvec = createMatrix(n, n);
    double* eigenval = createVector(n);

    HouseholderAlgorithm(A, eigenvec, eigenval, n);

    return 0;
}


void Tridiag(double** A, double** P, int n, int k)
/*  Tridiagonalizes a real symmetric nxn matrix A using the Householder algorithm if called repeatedly
    for k = 0, ... , n - 3.  P (passed as any nxn matrix) is the matrix acting the transformation.  */
{
    int i, j, m;
    double sum, mod, h;
    double* x = createVector(n);
    double* y = 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];
    }
    h = sum / 2;
    for (i = k + 1; i < n; i++) {
        for (j = k + 1; j < n; j++) {
            y[j] = A[j][i];
        }
        for (j = k + 1; j < n; j++) {
            sum = 0;
            for (m = k + 1; m < n; m++) {
                sum += ((j == m ? 1 : 0) - x[j] * x[m] / h) * y[m];
            }
            A[j][i] = sum;
        }
    }
    for (i = k + 1; i < n; i++) {
        for (j = k + 1; j < n; j++) {
            y[j] = A[i][j];
        }
        for (j = k + 1; j < n; j++) {
            sum = 0;
            for (m = k + 1; m < n; m++) {
                sum += y[m] * ((j == m ? 1 : 0) - x[m] * x[j] / h);
            }
            A[i][j] = sum;
        }
    }
    A[k + 1][k] = A[k][k + 1] = -mod;
    for (i = k + 2; i < n; i++) {
        A[i][k] = A[k][i] = 0;
    }
    for (i = k; i < n; i++) {
        for (j = k + 1; j < n; j++) {
            y[j] = P[i][j];
        }
        for (j = k + 1; j < n; j++) {
            sum = 0;
            for (m = k + 1; m < n; m++) {
                sum += y[m] * ((j == m ? 1 : 0) - x[m] * x[j] / h);
            }
            P[i][j] = sum;
        }
    }
    freeVector(x);
    freeVector(y);
}



void QRdecom(double** A, double** Q, double** P, int n)
/*  Gives the QR decomposition of a real symmetric matrix A. On output overwrites
    A with the upper triangular matrix R. Q (passed as the unit nxn matrix) is overwritten with the
    unitary matrix of the decomposition. On each call multiply the matrices P and Q and overwrites
    P with the result.  */
{
    int i, j;
    double theta, c, s, t, p, pp;
    for (i = 0; i <= n - 2; 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;
        for (j = 0; j <= n - 1; j++) {
            pp = P[j][i];
            P[j][i] = c * pp + s * P[j][i + 1];
            P[j][i + 1] = -s * pp + c * P[j][i + 1];
        }
    }
}


void QRinv(double** A, double** Q, int n)
/*  Takes the matrices Q and R (which is A on input) of the QR decomposition and 
    gives the inverse decomposition RQ overwriting A with the result.  */
{
    int j, k;
    double q, qq, qqq;
    for (j = 1; j <= n - 3; j++) {
        q = A[j][j];
        qq = A[j][j + 1];
        qqq = A[j][j + 2];
        A[j][j - 1] = q * Q[j][j - 1];
        A[j][j] = q * Q[j][j] + qq * Q[j + 1][j];
        A[j][j + 1] = q * Q[j][j + 1] + qq * Q[j + 1][j + 1] + qqq * Q[j + 2][j + 1];
    }
    q = A[0][0];
    qq = A[0][1];
    qqq = A[0][2];
    A[0][0] = q * Q[0][0] + qq * Q[1][0];
    A[0][1] = q * Q[0][1] + qq * Q[1][1] + qqq * Q[2][1];
    q = A[n - 2][n - 2];
    qq = A[n - 2][n - 1];
    A[n - 2][n - 3] = q * Q[n - 2][n - 3];
    A[n - 2][n - 2] = q * Q[n - 2][n - 2] + qq * Q[n - 1][n - 2];
    A[n - 2][n - 1] = q * Q[n - 2][n - 1] + qq * Q[n - 1][n - 1];
    q = A[n - 1][n - 1];
    A[n - 1][n - 2] = q * Q[n - 1][n - 2];
    A[n - 1][n - 1] = q * Q[n - 1][n - 1];
    for (j = 2; j < n; j++) {
        for (k = 0; k < j - 1; k++) {
            A[j][k] = A[k][j] = 0;
        }
    }
}


void QRdiag(double** A, double** P, int n)
/*  Diagonalizes a real symmetric tridiagonal matrix A (overwritten in output) with the QR algorithm.
    P (any nxn matrix can be passed) is the matrix performing the transformation,
    its column being the eigenvectors of the original matrix */
{
    int i, j, k;
    double** Q = createMatrix(n, n);
    for (k = 0; k < 30; k++) {
        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;
        }
        QRdecom(A, Q, P, n);
        QRinv(A, Q, n);
    }
    freeMatrix(Q, n);
}



void HouseholderAlgorithm(double** A, double** P, double* a, int n)
/*  Diagonalizes a real symmetric matrix A (overwritten on output) with the QR algorithm.
    P (passed as any nxn matrix) is the matrix of eigenvectors and a is the vector of eigenvalues. */
{
    int i, j;
    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            P[i][j] = (i == j ? 1 : 0);
        }
    }
    for (int i = 0; i < n - 2; i++) {
        Tridiag(A, P, n, i);
    }
    QRdiag(A, P, n);
    for (i = 0; i < n; i++) {
        a[i] = A[i][i];
    }
}

double** createMatrix(int r, int c)
{
    double** matrix = (double**)malloc(r * sizeof(double*));
    if (matrix) {
        for (int i = 0; i < r; i++) {
            matrix[i] = (double*)malloc(c * sizeof(double));
            if (matrix[i]) {}
            else {
                printf("Memory not allocated.\n");
                free(matrix[i]);
                exit(0);
            }
        }
    }
    else {
        printf("Memory not allocated.\n");
        free(matrix);
        exit(0);
    }
    return matrix;
}

void freeMatrix(double** a, int r)
{
    for (int i = 0; i < r; i++) {
            free(a[i]);
        }
        free(a);
}