truncate.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_TENSOR_TRAIN_ALGORITHMS_TRUNCATE_H
#define pRC_TENSOR_TRAIN_ALGORITHMS_TRUNCATE_H
 
#include <prc/algorithms/svd.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<Size C, Position P, class X, class R = RemoveReference<X>,
        If<IsTensorish<R>> = 0, class V = typename R::Value, class VT = V,
        If<All<IsFloat<V>, IsFloat<VT>>> = 0,
        If<IsSatisfied<(P == Position::Left || P == Position::Right)>> = 0>
    static inline constexpr auto truncate(X &&a,
        VT const &tolerance = NumericLimits<VT>::epsilon())
    {
        using DF = RemoveReference<decltype(folding<!P>(forward<X>(a)))>;
 
        if constexpr(C < min(DF::size(0), DF::size(1)))
        {
            if constexpr(P == Position::Right)
            {
                return expand(makeSeries<Index, typename R::Dimension() - 1>(),
                    [&](auto const... seq)
                    {
                        auto const [u, s, v] =
                            svd<C>(folding<!P>(forward<X>(a)), tolerance);
                        Tensor const lambda = norm(s);
 
                        return tuple(lambda,
                            eval(reshape<R::size(seq)...,
                                RemoveReference<decltype(u)>::size(1)>(u)),
                            eval(fromDiagonal(s / lambda) * adjoint(v)));
                    });
            }
            else
            {
                return expand(makeRange<Index, 1, typename R::Dimension{}>(),
                    [&](auto const... seq)
                    {
                        auto const [u, s, v] =
                            svd<C>(folding<!P>(forward<X>(a)), tolerance);
                        Tensor const lambda = norm(s);
 
                        return tuple(lambda, eval(u * fromDiagonal(s / lambda)),
                            eval(reshape<RemoveReference<decltype(v)>::size(1),
                                R::size(seq)...>(adjoint(v))));
                    });
            }
        }
        else
        {
            return orthogonalize<!P>(forward<X>(a));
        }
    }
}
#endif // pRC_TENSOR_TRAIN_ALGORITHMS_TRUNCATE_H