Встроенная C++ реализация QR-декомпозиции с помощью преобразования Householder

Question

(Результат: пара функций orhtonormalize а также forward_transformation отсюда получается проверенная реализация метода домохозяина).

Я пытаюсь реализовать алгоритм Householder для QR-разложения прямоугольной матрицы? На месте означает, что входные данные изменяются во время вычислений и диагонали верхней треугольной матрицы R, предоставляемой дополнительно, как это было представлено в статье на стр. 12 (или что-то похожее, но экономичное). Реализация функции tred2 (широко представленной в Интернете) не совсем ясна. Мое намерение состоит в том, чтобы реализовать (или получить готовую к использованию реализацию) преобразование Householder в чистом C++.

Это результат моих собственных усилий. Это наивное прямое переопределение кода MATLAB из упомянутой статьи, и (следовательно) оно дает в основном неправильный вывод:

#include <iostream>
#include <ostream>
#include <valarray>
#include <vector>
#include <deque>
#include <algorithm>
#include <numeric>
#include <utility>
#include <functional>
#include <limits>

#include <cstdlib>
#include <cassert>

using value_type = long double;
using vector = std::valarray< value_type >;
using matrix = std::valarray< vector >;

std::ostream &
operator << (std::ostream & out, matrix const & matrix)
{
    for (vector const & point : matrix) {
        for (value_type const & x : point) {
            out << x << ' ';
        }
        out << std::endl;
    }
    return out;
}

value_type const eps = std::numeric_limits< value_type >::epsilon();
value_type const zero = value_type(0);
value_type const one = value_type(1);

int
main()
{
    matrix qr{{12, 6, -4},
              {-51, 167, 24},
              {4, -68, -41}};
    std::cout << qr << std::endl;

    std::size_t const m = qr[0].size();
    std::size_t const n = qr.size();
    vector rtrace(zero, n);
    for (std::size_t j = 0; j < n; ++j) { // Householder itself
        value_type norm = zero;
        vector & qr_j = qr[j];
        for (std::size_t i = j; i < m; ++i) {
            value_type const & qrij = qr_j[i];
            norm += qrij * qrij;
        }
        using std::sqrt;
        norm = sqrt(norm);
        value_type & qrjj = qr_j[j];
        value_type & dj = rtrace[j];
        dj = (zero < qrjj) ? -norm : norm;
        using std::abs;
        value_type f = norm * (norm + abs(qrjj));
        assert(eps < f);
        f = one / sqrt(std::move(f));
        qrjj -= dj;
        for (std::size_t k = j; k < m; ++k) {
            qr_j[k] *= f;
        }
        for (std::size_t i = j + 1; i < n; ++i) {
            vector & qr_i = qr[i];
            value_type dot_product = zero;
            for (std::size_t k = j; k < m; ++k) {
                dot_product += qr_j[k] * qr_i[k];
            }
            for (std::size_t k = j; k < m; ++k) {
                qr_i[k] -= qr_j[k] * dot_product;
            }
        }
    }
    std::cout << "output:\n" << qr << std::endl;
    std::cout << "diagonal of R: " << std::endl;
    for (value_type const & rd : rtrace) {
        std::cout << rd << ' ';
    }
    std::cout << std::endl << std::endl;
    matrix q{{1, 0, 0},
             {0, 1, 0},
             {0, 0, 1}};
    for (std::size_t i = 0; i < n; ++i) {
        vector & qi = q[i];
        std::size_t j = n;
        while (0 < j) {
            --j;
            vector & qr_j = qr[j];
            value_type s_ = zero;
            for (std::size_t k = j; k < m; ++k) {
                s_ += qr_j[k] * qi[k];
            }
            for (std::size_t k = j; k < m; ++k) {
                qi[k] += qr_j[k] * s_;
            }
        }
    }
    std::cout << "Q:\n" << q << std::endl << std::endl;
    matrix r{{0, 0, 0},
             {0, 0, 0},
             {0, 0, 0}};
    for (std::size_t i = 0; i < n; ++i) {
        r[i][i] = rtrace[i];
        for (std::size_t j = i + 1; j < n; ++j) {
            r[j][i] = qr[j][i];
        }
    }
    std::cout << "R:\n" << r << std::endl << std::endl;
    matrix a{{0, 0, 0},
             {0, 0, 0},
             {0, 0, 0}};
    for (std::size_t i = 0; i < n; ++i) {
        for (std::size_t j = 0; j < n; ++j) {
            a[j][i] += q[i][j] * r[i][j];
        }
    }
    std::cout << "regenerated input:\n" << a << std::endl << std::endl;
    return EXIT_SUCCESS;
}