Spectral Decomposition

Spectral Decomposition, sometimes called Eigen Decomposition factors a matrix into a canonical form so that it is represented in terms of eigenvectors and eigenvalues. Spectral decomposition takes the form \(S = UKU^{T}\), where \(K\) is a matrix whos main diagonal contains the eigenvalues of \(S\) and the basis vectors of \(U\) are the eigenvectors of \(S\).

This decomposition is not unique. The rows of both \(U\) and \(K\) can be re-arranged in any order. The re-arranged rows can even be negated (so long as it represents a valid rotation). The formula \(S = UKU^{T}\) will remain true after doing any of these modifications to the matrix.

A matrix can be factored into eigenvectors and eigenvalues using the QR Algorithm. Like Polar Decomposition, Spectral Decompositon is also an iterative process. Also like Polar Decomposition, spectral decompositon should not need more than 20 iterations. I found that most real world matrices take about 5 to 8 iterations.

To implement the QR algorithm, set a helper matrix \(A\) to be the input matrix \(S\). Iterativeley, find the QR Factorization of \(A = QR\). Take the resulting QR factorization and set \(A = RQ\). Repeat until the lower triangular components of \(A\) are \(\approx 0\). The main diagonal of the resulting matrix \(A\) contains the eigenvectors of \(Q\).

That takes care of the eigenvalues, but what about the eigenvectors? Turns out the matrix \(U\) is an accumulation of all the orthogonal transforms taken to find \(K\). At the begenning of the algorithm set \(U\) to identity. On every iteration, multiply the \(Q\) result of the QR factorization into \(U\). The following code implements this.

struct SpectralDecompositionResult {
    Matrix U; // Each basis vector is an eigenvector
    Matrix K; // Contains eigenvalues on main diagonal
    Matrix Ut; // Transpose of U
}

bool SpectralDecompositionEarlyOut(Matrix A) {
    int indices[] = [4, 8, 9, 12, 13, 14];

    for (int i = 0; i < indices.count; ++i) {
        if (!EpsilonCompare(A[i], 0)) {
            return false;
        }
    }
    return true;
}

SpectralDecompositionResult SpectralDecomposition(S) {
    SpectralDecompositionResult result;

    result.U = Matrix();
    Matrix A = S;

    for (int i = 0; i < 20; ++i) {
        QRDecompResult decomp = QRDecomposition(A);

        result.U = result.U * decomp.Q;
        A = decomp.R * decomp.Q;

        if (SpectralDecompositionEarlyOut(A)) {
            break;
        }
    }

    result.Ut = Transpose(result.U);

    result.K = Matrix();
    result.K[0] = A[0];
    result.K[1] = A[1];
    result.K[2] = A[2];

    return result;
}

The diagonal elements matrix \(K\) are the eigen values, in descending order. While this is algebraically correct, geometrically it doesn't hold any meaning. We need to adjust the result of the spectral decomposition to be geometrically intuitive. The next section describes how to do this.