householder.hpp

Go to the documentation of this file.

/*
 * Copyright © 2025-2026 Akreon LLC
 * SPDX-License-Identifier: Apache-2.0
 */
 
#pragma once
 
#include <algorithm>
 
#include "akreon/linalg/assert.hpp"
#include "akreon/linalg/linalg_ops.hpp"
#include "akreon/linalg/linalg_util.hpp"
#include "akreon/linalg/norms/vector_norms.hpp"
 
namespace linalg {
template <typename T, std::size_t N>
  requires(Field<T> && N > 0)
struct HouseholderVector {
  constexpr HouseholderVector(const Vector<T, N>& x) : beta(real_type_t<T>(0)) {
    make_householder(x, std::integral_constant<std::size_t, N>{});
  }
 
  Vector<T, N> v;      
  real_type_t<T> beta; 
 
private:
  constexpr void make_householder(const Vector<T, N>& x, auto) {
    // Scale x before computing its norm and reflector direction.
    real_type_t<T> scale = 0;
    for (const auto& value : x) scale = std::max(scale, constexpr_abs(value));
    v(0) = T(1);
    if (scale == real_type_t<T>(0)) return;
    for (std::size_t i = 0; i < N; ++i) v(i) = x(i) / scale;
    const auto norm = L2<Vector<T, N>>()(v);
    if constexpr (Complex<T>) {
      // Choose the leading phase so x + phase*||x||*e_0 avoids cancellation.
      // Normalize the phase at unit scale. The constexpr magnitude's
      // absolute square-root tolerance otherwise distorts tiny leading
      // entries and leaves a residual that packed QR silently discards.
      const auto phase_scale =
          std::max(constexpr_abs(v(0).real()), constexpr_abs(v(0).imag()));
      const auto leading =
          phase_scale == real_type_t<T>(0) ? T(1) : v(0) / phase_scale;
      const auto phase = leading / constexpr_abs(leading);
      v(0) += phase * norm;
      // QR and Hessenberg packing reconstruct the leading entry as one.
      v = (T(1) / v(0)) * v;
      const auto vv = inner_product(v, v);
      beta = real_type_t<T>(2) / vv.real();
    } else {
      const auto tail = get_elements<1, N>(v);
      if (inner_product(tail, tail) == T(0)) {
        beta = x(0) < T(0) ? T(2) : T(0);
        v(0) = T(1);
        return;
      }
      // Use the opposite-sign reflection and normalize its leading entry to
      // one.
      const auto leading = v(0) + (v(0) >= T(0) ? norm : -norm);
      v = (T(1) / leading) * v;
      v(0) = T(1);
      beta = T(2) / inner_product(v, v);
    }
  }
 
  constexpr void make_householder(const Vector<T, N>&,
                                  std::integral_constant<std::size_t, 1>) {
    beta = real_type_t<T>(0);
    v(0) = T(1);
  }
};
 
} // namespace linalg