svd.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_ALGORITHMS_SVD_H
#define pRC_ALGORITHMS_SVD_H
 
#include <prc/config.hpp>
#include <prc/algorithms/jacobi_rotation.hpp>
#include <prc/core/basic/functions/min.hpp>
#include <prc/core/functors/max.hpp>
#include <prc/core/functors/min.hpp>
#include <prc/core/tensor/functions/chip.hpp>
#include <prc/core/tensor/functions/reduce.hpp>
#include <prc/core/tensor/functions/slice.hpp>
#include <prc/core/tensor/operator/functions/extract_diagonal.hpp>
#include <prc/core/tensor/operator/functions/from_diagonal.hpp>
#include <prc/core/tensor/operator/functions/is_unitary.hpp>
#include <prc/core/tensor/type_traits.hpp>
#include <prc/core/value/functions/is_approx.hpp>
#include <prc/core/value/limits.hpp>
 
namespace pRC
{
    template<class X, class R = RemoveReference<X>, If<IsTensorish<R>> = 0,
        If<IsFloat<typename R::Value>> = 0,
        If<IsSatisfied<(typename R::Dimension() == 2)>> = 0>
    static inline constexpr auto svd(X &&input)
    {
        using T = typename R::Type;
        constexpr auto M = R::size(0);
        constexpr auto N = R::size(1);
        constexpr auto D = min(M, N);
 
        Tensor scale = reduce<Max>(abs(input));
        if(abs(scale)() <= NumericLimits<typename T::Value>::min())
        {
            scale = identity();
        }
 
        tuple<Tensor<T, M, D>, Tensor<typename T::NonComplex, D>,
            Tensor<T, N, D>>
            result;
        auto &u = get<0>(result);
        auto &s = get<1>(result);
        auto &v = get<2>(result);
        Tensor<T, D, D> a;
 
        if constexpr(M == N)
        {
            a = input / scale;
            u = identity();
            v = identity();
        }
        else if constexpr(M > N)
        {
            tie(u, a) = qr(input / scale);
            v = identity();
        }
        else if constexpr(M < N)
        {
            tie(v, a) = qr(adjoint(input / scale));
            a = eval(adjoint(a));
            u = identity();
        }
 
        auto maxDiagEntry = reduce<Max>(extractDiagonal(abs(a)))();
 
        auto const updateMaxDiagEntry =
            [](auto const &a, Index const p, Index const q, auto const &old)
        {
            return max(abs(a(p, p)), abs(a(q, q)), old);
        };
 
        auto const threshold = [](auto const &diag)
        {
            return max(NumericLimits<typename T::Value>::min(),
                NumericLimits<typename T::Value>::epsilon() *
                    identity<typename T::Value>(2) * diag);
        };
 
        Bool finished = false;
        while(!finished)
        {
            finished = true;
            range<1, D>(
                [&](auto const p)
                {
                    for(Index q = 0; q < p; ++q)
                    {
                        if(max(abs(a(p, q)), abs(a(q, p))) >
                            threshold(maxDiagEntry))
                        {
                            finished = false;
 
                            if constexpr(typename T::IsComplexified{})
                            {
                                if(auto const n = sqrt(norm<2, 1>(a(p, p)) +
                                       norm<2, 1>(a(q, p)));
                                    n == zero())
                                {
                                    a(p, p) = zero();
                                    a(q, p) = zero();
 
                                    if(abs(imag(a(p, q))) >
                                        NumericLimits<typename T::Value>::min())
                                    {
                                        auto const z = abs(a(p, q)) / a(p, q);
                                        chip<0>(a, p) *= z;
                                        chip<1>(u, p) *= conj(z);
                                    }
                                }
                                else
                                {
                                    JacobiRotation rot(conj(a(p, p)) / n,
                                        a(q, p) / n);
                                    apply(rot, a, p, q);
                                    apply(u, adjoint(rot), p, q);
 
                                    if(abs(imag(a(p, q))) >
                                        NumericLimits<typename T::Value>::min())
                                    {
                                        auto const z = abs(a(p, q)) / a(p, q);
                                        chip<1>(a, q) *= z;
                                        chip<1>(v, q) *= z;
                                    }
                                }
                                if(abs(imag(a(q, q))) >
                                    NumericLimits<typename T::Value>::min())
                                {
                                    auto const z = abs(a(q, q)) / a(q, q);
                                    chip<0>(a, q) *= z;
                                    chip<1>(u, q) *= conj(z);
                                }
 
                                maxDiagEntry =
                                    updateMaxDiagEntry(a, p, q, maxDiagEntry);
 
                                if(max(abs(a(p, q)), abs(a(q, p))) <=
                                    threshold(maxDiagEntry))
                                {
                                    continue;
                                }
                            }
 
                            Tensor<typename T::NonComplex, 2, 2> m;
                            m(0, 0) = real(a(p, p));
                            m(1, 0) = real(a(q, p));
                            m(0, 1) = real(a(p, q));
                            m(1, 1) = real(a(q, q));
 
                            auto const t = m(0, 0) + m(1, 1);
                            auto const d = m(1, 0) - m(0, 1);
                            auto const rot =
                                JacobiRotation<T>::MakeGivens(t, d);
 
                            apply(rot, m, 0, 1);
                            auto const jRight = JacobiRotation(m, 0, 1);
                            auto const jLeft = rot * transpose(jRight);
 
                            apply(jLeft, a, p, q);
                            apply(a, jRight, p, q);
                            apply(u, transpose(jLeft), p, q);
                            apply(v, jRight, p, q);
 
                            maxDiagEntry =
                                updateMaxDiagEntry(a, p, q, maxDiagEntry);
                        }
                    }
                });
        }
 
        range<D>(
            [&a, &u, &s](auto const i)
            {
                if(abs(imag(a(i, i))) > NumericLimits<typename T::Value>::min())
                {
                    s(i) = abs(a(i, i));
                    chip<1>(u, i) *= a(i, i) / s(i);
                }
                else
                {
                    auto r = real(a(i, i));
                    s(i) = abs(r);
                    if(r < zero())
                    {
                        chip<1>(u, i) = -chip<1>(u, i);
                    }
                }
            });
 
        s *= scale;
 
        for(Index i = 0; i < D; ++i)
        {
            Index iMax = i;
            for(Index p = i + 1; p < D; ++p)
            {
                if(s(p) > s(iMax))
                {
                    iMax = p;
                }
            }
 
            if(s(iMax) == zero())
            {
                break;
            }
 
            if(iMax != i)
            {
                swap(s(i), s(iMax));
                swap(chip<1>(u, i), chip<1>(u, iMax));
                swap(chip<1>(v, i), chip<1>(v, iMax));
            }
        }
 
        if constexpr(cDebugLevel >= DebugLevel::High)
        {
            if(!isApprox(input, u * fromDiagonal(s) * adjoint(v)))
            {
                Logging::error("SVD decomposition failed.");
            }
            if(!isUnitary(u))
            {
                Logging::error("SVD decomposition failed: U is not unitary.");
            }
            if(!isUnitary(v))
            {
                Logging::error("SVD decomposition failed: V is not unitary.");
            }
            if(reduce<Min>(s)() < zero())
            {
                Logging::error(
                    "SVD decomposition failed: Found negative singular "
                    "values.");
            }
        }
 
        return result;
    }
 
    template<Size C, 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<(typename R::Dimension() == 2)>> = 0,
        If<IsSatisfied<(C < min(R::size(0), R::size(1)))>> = 0>
    static inline constexpr auto svd(X &&input,
        VT const &tolerance = NumericLimits<VT>::tolerance())
    {
        auto const [fU, fS, fV] = svd(forward<X>(input));
 
        constexpr auto M = decltype(fU)::size(0);
        constexpr auto D = decltype(fS)::size();
        constexpr auto N = decltype(fV)::size(0);
 
        auto const sNorm = norm<2, 1>(fS)();
        auto sError = norm<2, 1>(slice<D - C>(fS, C))();
 
        auto result = tuple(eval(slice<M, C>(fU, 0, 0)), eval(slice<C>(fS, 0)),
            eval(slice<N, C>(fV, 0, 0)));
 
        auto &u = get<0>(result);
        auto &s = get<1>(result);
        auto &v = get<2>(result);
 
        if(sError > square(tolerance) * sNorm)
        {
            /*
            Logging::info(
                "Truncated SVD: Relative error is larger than requested "
                "tolerance due to specified cutoff:",
                sqrt(sError / sNorm));
            */
 
            return result;
        }
 
        for(Index k = C; k != 0;)
        {
            if(sError + square(s(--k)) > square(tolerance) * sNorm)
            {
                break;
            }
 
            sError += square(s(k));
 
            s(k) = zero();
            chip<1>(u, k) = zero();
            chip<1>(v, k) = zero();
        }
 
        return result;
    }
 
    template<Size C, 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<(typename R::Dimension() == 2)>> = 0,
        If<IsSatisfied<(C == min(R::size(0), R::size(1)))>> = 0>
    static inline constexpr auto svd(X &&input,
        VT const &tolerance = NumericLimits<VT>::tolerance())
    {
        auto result = svd(forward<X>(input));
        auto &u = get<0>(result);
        auto &s = get<1>(result);
        auto &v = get<2>(result);
 
        auto const sNorm = norm<2, 1>(s)();
        auto sError = zero<decltype(sNorm)>();
 
        for(Index k = C; k != 0;)
        {
            if(sError + square(s(--k)) > square(tolerance) * sNorm)
            {
                break;
            }
 
            sError += square(s(k));
 
            s(k) = zero();
            chip<1>(u, k) = zero();
            chip<1>(v, k) = zero();
        }
 
        return result;
    }
}
#endif // pRC_ALGORITHMS_SVD_H