qr.hpp

#ifndef QR_DECOMP_HPP
#define QR_DECOMP_HPP
 
#include "../../math/constmath.hpp"
#include "../matrix.hpp"
#include "../operations.hpp"
 
namespace consteig
{
 
 
template <typename T, Size S> struct QRMatrix
{
    Matrix<T, S, S> _q;
    Matrix<T, S, S> _r;
};
 
// Algorithm: QR Decomposition (Hessenberg Optimized)
// Factors an upper Hessenberg matrix A = QR using a sequence of Givens
// rotations. Because A is Hessenberg, each column has at most one subdiagonal
// nonzero, so each rotation only needs to update two adjacent rows of R and
// two adjacent columns of Q.
//
// The rotation angle is chosen (Golub & Van Loan, sec. 5) so that the
// subdiagonal entry r(i+1, i) is zeroed:
//   c = x / mag,  s = y / mag,  mag = sqrt(x^2 + y^2)
// where x = r(i,i) and y = r(i+1,i).
template <typename T, Size R, Size C>
constexpr QRMatrix<T, R> qr_hessenberg(const Matrix<T, R, C> a)
{
    static_assert(R == C, "QR decomposition must be a square matrix");
 
    Matrix<T, R, R> q = eye<T, R>();
    Matrix<T, R, R> r = a;
 
    for (Size col = 0; col < R - 1; ++col)
    {
        // Compute Givens rotation to zero r(col+1, col)
        T x = r(col, col);
        T y = r(col + 1, col);
 
        if (abs(y) > static_cast<T>(1e-15))
        {
            T mag = sqrt(x * x + y * y);
            T c = x / mag;
            T s = y / mag;
 
            // R = G * R
            for (Size k = col; k < R; ++k)
            {
                T r_cur = r(col, k);
                T r_next = r(col + 1, k);
                r(col, k) = c * r_cur + s * r_next;
                r(col + 1, k) = -s * r_cur + c * r_next;
            }
 
            // Q = Q * G^T
            for (Size row = 0; row < R; ++row)
            {
                T q_cur = q(row, col);
                T q_next = q(row, col + 1);
                q(row, col) = c * q_cur + s * q_next;
                q(row, col + 1) = -s * q_cur + c * q_next;
            }
            r(col + 1, col) = 0;
        }
    }
 
    QRMatrix<T, R> res{};
    res._q = q;
    res._r = r;
    return res;
}
 
// Algorithm: QR Decomposition (General)
// Factors a square matrix A = QR using a sequence of Givens rotations, working
// column by column from bottom to top to zero all subdiagonal entries. Givens
// rotations are preferred over Gram-Schmidt here because they avoid the loss of
// orthogonality that Gram-Schmidt can suffer from floating-point cancellation.
//
// Reference: Golub & Van Loan, "Matrix Computations" (4th ed.), sec. 5
template <typename T, Size R, Size C>
constexpr QRMatrix<T, R> qr(const Matrix<T, R, C> a)
{
    static_assert(R == C, "QR decomposition must be a square matrix");
 
    Matrix<T, R, R> q = eye<T, R>();
    Matrix<T, R, R> r = a;
 
    for (Size col = 0; col < C; ++col)
    {
        for (Size row = R - 1; row > col; --row)
        {
            // Compute Givens rotation to zero r(row, col) using r(row-1, col)
            T x = r(row - 1, col);
            T y = r(row, col);
 
            if (abs(y) > static_cast<T>(1e-15))
            {
                T mag = sqrt(x * x + y * y);
                if (mag > static_cast<T>(1e-18))
                {
                    T c = x / mag;
                    T s = y / mag;
 
                    // R = G * R
                    // Only need to update columns from col to C-1
                    for (Size k = col; k < C; ++k)
                    {
                        T r_above = r(row - 1, k);
                        T r_cur = r(row, k);
                        r(row - 1, k) = c * r_above + s * r_cur;
                        r(row, k) = -s * r_above + c * r_cur;
                    }
 
                    // Q = Q * G^T
                    for (Size k = 0; k < R; ++k)
                    {
                        T q_above = q(k, row - 1);
                        T q_cur = q(k, row);
                        q(k, row - 1) = c * q_above + s * q_cur;
                        q(k, row) = -s * q_above + c * q_cur;
                    }
                    r(row, col) = 0;
                }
            }
        }
    }
 
    QRMatrix<T, R> res{};
    res._q = q;
    res._r = r;
    return res;
}
 
 
} // namespace consteig
#endif