QR Decomposition

QR Decomposition, also called QR Factorization is the process of breaking down a matrix into the product of an orthogonal matrix and an upper triangular matrix. The formula for this is expressed as \(A=QR\) where \(A\) is the input matrix, \(Q\) is an orthogonal matrix and \(R\) is an upper triangular matrix.

QR Decomposition is important because it is the basis of the QR Algorithm, the method we will use to find the eigenvectors and eignvalues of a matrix. The variables \(Q\), \(A\) and \(R\) presented in this page are not a part of the decomposition method described before.

\(Q\) is simply the ortho-normal basis of \(A\). To get the ortho normal basis of \(A\) first orthogonalize the basis vectors matrix using the Gram Schmidt algorithm. Once the basis vectors are orthogonal, normalize them. That's it, we now know what \(Q\) is.

Remember, we're trying to solve for \(A=QR\). To find \(R\) multiply both sides of the equation by the transpose of \(Q\) like so

$$ Q^{T}A = Q^{T}QR $$

Because \(Q\) is ortho normal, \(Q^{T}\) is the same as the inverse of \(Q\). The product of multiplying a matrix by its inverse is the identity matrix. The \(Q^{T}Q\) part of the right hand side of the above equation simplifies to the identity matrix \(I\). This means the equation becomes:

$$ Q^{T}A = R $$

Now that we know that \(R\) equals \(Q^{T}A\), we can write the QRDecomposition function.

struct QRDecompResult {
    Matrix Q; // Q is an orthogonal matrix
    Matrix R; // R os an upper triangular matrix
}

QRDecompResult QRDecomposition(Matrix A) {
    QRDecompResult result;

    Vector x(A[0], A[1], A[2]);
    Vector y(A[4], A[5], A[6]);
    Vector z(A[8], A[9], A[10]);

    y = y -  ProjectionV3(x, y); // y minus projection of y onto x
    z = (z - Projection(x, z)) - Projection(y, z);

    x = Normalize(x);
    y = Normalize(y);
    z = Normalize(z);

    result.Q[0] = x[0]; resuly.Q[1] = x[1]; result.Q[2] = x[2];
    result.Q[4] = y[0]; result.Q[5] = y[1]; result.Q[6] = y[2];
    result.Q[8] = z[0]; result.Q[9] = z[1]; result.Q[10]= z[2];

    result.R = Transpose(result.Q) * A;

    return result;
}