orthogonalize.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_TENSOR_TRAIN_ALGORITHMS_ORTHOGONALIZE_H
#define pRC_TENSOR_TRAIN_ALGORITHMS_ORTHOGONALIZE_H
 
#include <prc/algorithms/lq.hpp>
#include <prc/algorithms/qr.hpp>
#include <prc/core/basic/position.hpp>
#include <prc/core/tensor/functions/norm.hpp>
#include <prc/core/tensor/functions/reshape.hpp>
#include <prc/core/tensor/type_traits.hpp>
 
namespace pRC
{
    template<Position P, class X, class R = RemoveReference<X>,
        If<IsTensorish<R>> = 0, If<IsFloat<typename R::Value>> = 0,
        If<IsSatisfied<(P == Position::Left || P == Position::Right)>> = 0>
    static inline constexpr auto orthogonalize(X &&a)
    {
        if constexpr(P == Position::Left)
        {
            return expand(makeSeries<Index, typename R::Dimension() - 1>(),
                [&](auto const... seq)
                {
                    auto const [q, r] = qr(folding<P>(forward<X>(a)));
                    Tensor const lambda = norm(r);
 
                    return tuple(lambda,
                        eval(reshape<R::size(seq)...,
                            RemoveReference<decltype(q)>::size(1)>(q)),
                        eval(r / lambda));
                });
        }
        else
        {
            return expand(makeRange<Index, 1, typename R::Dimension{}>(),
                [&](auto const... seq)
                {
                    auto const [l, q] = lq(folding<P>(forward<X>(a)));
                    Tensor const lambda = norm(l);
 
                    return tuple(lambda, eval(l / lambda),
                        eval(reshape<RemoveReference<decltype(q)>::size(0),
                            R::size(seq)...>(q)));
                });
        }
    }
}
#endif // pRC_TENSOR_TRAIN_ALGORITHMS_ORTHOGONALIZE_H