least_squares.h

Go to the documentation of this file.

#ifndef FUSED_ALGORITHMS_LEAST_SQUARES_H
#define FUSED_ALGORITHMS_LEAST_SQUARES_H
 
#include "config.h"
#include "utilities/expected.h"
#include "matrix_traits.h"
#include "fused/fused_matrix.h"
#include "fused/fused_vector.h"
#include "algorithms/decomposition/qr.h"
#include "algorithms/solvers/triangular_solve.h"
 
namespace matrix_algorithms
{
 
    using matrix_traits::MatrixStatus;
 
    template <typename T, my_size_t M, my_size_t N>
    Expected<FusedVector<T, N>, MatrixStatus> least_squares(
        const FusedMatrix<T, M, N> &A,
        const FusedVector<T, M> &b)
    {
        static_assert(is_floating_point_v<T>,
                      "least_squares requires a floating-point scalar type");
        static_assert(M >= N, "least_squares requires M >= N (overdetermined system)");
 
        // 1. QR decompose A
        auto qr = qr_householder(A);
 
        // 2. Compute c = Qᵀ·b
        auto c = qr.apply_Qt(b);
 
        // 3. Extract R₁ (top N×N of R) and c₁ (top N of c)
        FusedMatrix<T, N, N> R1(T(0));
 
        for (my_size_t i = 0; i < N; ++i)
        {
            for (my_size_t j = i; j < N; ++j)
            {
                R1(i, j) = qr.QR(i, j);
            }
        }
 
        FusedVector<T, N> c1(T(0));
 
        for (my_size_t i = 0; i < N; ++i)
        {
            c1(i) = c(i);
        }
 
        // 4. Back-substitute R₁·x = c₁
        return back_substitute(R1, c1);
    }
 
} // namespace matrix_algorithms
 
#endif // FUSED_ALGORITHMS_LEAST_SQUARES_H