qr.cpp

Go to the documentation of this file.

#include "linalg/factorization/qr.hpp"
#include <cmath>
#include <algorithm>
#include <vector>
 
namespace num {
 
// qr()   --  Householder reflections
//
// A Householder reflector is H = I - 2*v*v^T where v is a unit vector.
// H is symmetric and orthogonal: H^T = H,  H^2 = I,  det(H) = -1.
//
// For a vector x, choose v so that H*x = -sign(x[0])*||x||*e_1.
// This zeros all entries of x below the first.
//
// SIGN CONVENTION
// v = x + sign(x[0]) * ||x|| * e_1
//
// Using the same sign as x[0] avoids catastrophic cancellation.
// If x[0] > 0: v[0] = x[0] + ||x||  (sum of positives  -- well conditioned)
// If x[0] < 0: v[0] = x[0] - ||x||  (sum of negatives  -- still fine)
// Contrast with v[0] = x[0] - ||x|| when x[0] > 0: cancellation when x[0] ~= ||x||.
//
// APPLYING H TO A SUBMATRIX (RANK-1 UPDATE)
// H * A[k:m, k:n] = A[k:m, k:n] - 2*v*(v^T * A[k:m, k:n])
//
// This is an outer product update: one inner product per column, one axpy.
// Cost: O(m*n) per step.  Total cost: O(m*n^2 - n^3/3).
//
// BUILDING Q EXPLICITLY
// Q = H_0 * H_1 * ... * H_{r-1}   where r = min(m-1, n)
//
// Build by applying reflectors in reverse order to the identity:
//   Q = I
//   for k = r-1 downto 0:
//     Q[k:m, k:m] = H_k * Q[k:m, k:m]   (only columns >= k are non-zero)
//
// Correctness: After applying H_{r-1},...,H_{k+1}, Q[:,j<k] = e_j (untouched
// by all previous steps since H_i only acts on rows >= i > j).
 
QRResult qr(const Matrix& A) {
    constexpr real householder_tol = 1e-14;
    const idx m = A.rows();
    const idx n = A.cols();
    const idx r = (m > n) ? n : m - 1;    // number of Householder steps
 
    Matrix R = A;   // working copy -> becomes R at the end
 
    // Store the normalised Householder vectors for later Q accumulation.
    // vs[k] has length m-k; it spans the rows k..m-1.
    std::vector<std::vector<real>> vs(r);
 
    for (idx k = 0; k < r; ++k) {
        const idx len = m - k;
 
        std::vector<real> x(len);
        for (idx i = 0; i < len; ++i) x[i] = R(k + i, k);
 
        real norm_x = real(0);
        for (idx i = 0; i < len; ++i) norm_x += x[i] * x[i];
        norm_x = std::sqrt(norm_x);
 
        std::vector<real> v = x;
        v[0] += (x[0] >= real(0)) ? norm_x : -norm_x;   // sign trick
 
        real norm_v = real(0);
        for (idx i = 0; i < len; ++i) norm_v += v[i] * v[i];
        norm_v = std::sqrt(norm_v);
 
        if (norm_v < householder_tol) {
            // Column already zero below diagonal; Householder is identity.
            vs[k].assign(len, real(0));
            continue;
        }
        for (idx i = 0; i < len; ++i) v[i] /= norm_v;
        vs[k] = v;
 
        // For each column j: R[k:m, j] -= 2 * v * (v^T * R[k:m, j])
        for (idx j = k; j < n; ++j) {
            real vTr = real(0);
            for (idx i = 0; i < len; ++i) vTr += v[i] * R(k + i, j);
            const real two_vTr = real(2) * vTr;
            for (idx i = 0; i < len; ++i) R(k + i, j) -= two_vTr * v[i];
        }
    }
 
    Matrix Q(m, m, real(0));
    for (idx i = 0; i < m; ++i) Q(i, i) = real(1);   // start with identity
 
    for (idx k = r; k-- > 0; ) {
        const std::vector<real>& v = vs[k];
        const idx len = static_cast<idx>(v.size());
 
        // Apply H_k to Q[k:m, k:m]:  Q[k:m, j] -= 2*v*(v^T * Q[k:m, j])
        // Only columns j >= k can be non-zero at this point (see note above).
        for (idx j = k; j < m; ++j) {
            real vTq = real(0);
            for (idx i = 0; i < len; ++i) vTq += v[i] * Q(k + i, j);
            const real two_vTq = real(2) * vTq;
            for (idx i = 0; i < len; ++i) Q(k + i, j) -= two_vTq * v[i];
        }
    }
 
    for (idx i = 1; i < m; ++i)
        for (idx j = 0; j < std::min(i, n); ++j)
            R(i, j) = real(0);
 
    return {std::move(Q), std::move(R)};
}
 
// qr_solve()   --  least-squares via Q^T multiply + back-substitution
//
// min ||A*x - b||  <=>  R[:n,:n] * x = (Q^T * b)[:n]
//
// Proof: ||A*x - b||^2 = ||Q*R*x - b||^2 = ||R*x - Q^T*b||^2
//                      = ||R[:n,:n]*x - c[:n]||^2 + ||c[n:]||^2
// where c = Q^T*b.  The second term is fixed (residual); minimise the first.
 
void qr_solve(const QRResult& f, const Vector& b, Vector& x) {
    const idx m = f.Q.rows();
    const idx n = f.R.cols();
 
    // Q^T[i,j] = Q[j,i] so y[i] = sum_j Q[j,i] * b[j]
    Vector y(m, real(0));
    for (idx i = 0; i < m; ++i)
        for (idx j = 0; j < m; ++j)
            y[i] += f.Q(j, i) * b[j];
 
    Vector xv(n, real(0));
    for (idx i = n; i-- > 0; ) {
        xv[i] = y[i];
        for (idx j = i + 1; j < n; ++j)
            xv[i] -= f.R(i, j) * xv[j];
        xv[i] /= f.R(i, i);
    }
 
    x = std::move(xv);
}
 
} // namespace num