norm.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_TENSOR_TRAIN_COMMON_FUNCTIONS_NORM_H
#define pRC_TENSOR_TRAIN_COMMON_FUNCTIONS_NORM_H
 
#include <prc/core/basic/position.hpp>
#include <prc/core/tensor/declarations.hpp>
#include <prc/tensor_train/operator/declarations.hpp>
#include <prc/tensor_train/tensor/declarations.hpp>
 
namespace pRC
{
    template<Index P = 2, Index Q = P, class X, class R = RemoveReference<X>>
        requires(TensorTrain::IsTensorish<R> ||
                    TensorTrain::IsOperatorish<R>) &&
        (P == 2 && Q == 2)
    static inline constexpr auto norm(X &&a)
    {
        using T = NonComplex<Value<R>>;
        T lambda = zero();
 
        return RecursiveLambda(
            [&a, &lambda]<Index I = 0>(auto const &self, auto const &p) {
                auto const f = [&p, &a, &lambda]()
                {
                    Tensor const tmp = contract<1, 0>(p, a.template core<I>());
                    auto [q, r] = orthogonalize<Position::Left>(tmp);
 
                    auto const rNorm = norm(r)();
                    if(rNorm > NumericLimits<T>::min())
                    {
                        lambda += log(rNorm);
                        r /= rNorm;
                    }
 
                    return r;
                };
 
                if constexpr(I == R::Dimension - 1)
                {
                    return exp(lambda) * reshape<>(f());
                }
                else
                {
                    return self.template operator()<I + 1>(f());
                }
            })(identity<Tensor<T, 1, 1>>());
    }
}
#endif // pRC_TENSOR_TRAIN_COMMON_FUNCTIONS_NORM_H