slogdet.hpp

Go to the documentation of this file.

#pragma once
 
#include "lu.hpp"
 
namespace uni20::linalg
{
template <RealOrComplex S> struct LogDeterminant
{
    S phase;
    make_real_t<S> log_absolute;
};
 
template <RealOrComplex S> struct LogDeterminantResult
{
    SolveInfo info;
    std::optional<LogDeterminant<S>> value;
};
 
namespace detail
{
template <RealOrComplex S> LogDeterminant<S> zero_determinant(Precision p)
{
  auto zero = scalar_at_precision<S>(0, p);
  return {.phase = zero, .log_absolute = math::log(uni20::real(zero))};
}
 
template <RealOrComplex S> LogDeterminantResult<S> reduce_log_determinant(LuFactorization<S> const& factor, Precision p)
{
  auto phase = scalar_at_precision<S>(1, p);
  auto sum = uni20::real(scalar_at_precision<S>(0, p));
  auto correction = sum;
  auto const& a = factor.packed();
  for (std::size_t k = 0; k < factor.order(); ++k)
  {
    auto diagonal = a[k, k];
    auto magnitude = lu_magnitude(diagonal);
    auto term = math::log(magnitude.scale) + math::log(magnitude.unit);
    // Neumaier reduction retains small contributions when logarithms cancel.
    auto next = sum + term;
    if (math::abs(sum) >= math::abs(term))
      correction += (sum - next) + term;
    else
      correction += (term - next) + sum;
    sum = std::move(next);
    if constexpr (Complex<S>)
    {
      phase *= (diagonal / magnitude.scale) / magnitude.unit;
      // Avoid accumulated drift of the product away from unit modulus.
      auto phase_magnitude = lu_magnitude(phase);
      phase = (phase / phase_magnitude.scale) / phase_magnitude.unit;
    }
    else if (diagonal < 0)
      phase = -phase;
    if (factor.pivots()[k] != k) phase = -phase;
    if (!uni20::isfinite(sum) || !uni20::isfinite(correction) || !uni20::isfinite(phase))
      return {.info = {.status = SolveStatus::nonfinite_result}, .value = {}};
  }
  sum += correction;
  if (!uni20::isfinite(sum)) return {.info = {.status = SolveStatus::nonfinite_result}, .value = {}};
  return {.info = {}, .value = LogDeterminant<S>{std::move(phase), std::move(sum)}};
}
} // namespace detail
 
template <RealOrComplex S> [[nodiscard]] LogDeterminantResult<S> slogdet_with_info(LuFactorization<S> const& factor)
{
  return detail::reduce_log_determinant(factor, factor.precision());
}
 
namespace detail
{
template <KernelBackendSelector Backend, RankedTensorView<2> A>
auto slogdet_matrix_impl(Backend&& backend, A const& a, SolveOptions<make_real_t<tensor_element_t<A>>> const& options,
                         Precision p)
{
  using S = tensor_element_t<A>;
  auto attempt = lu_factor_impl(
      std::forward<Backend>(backend), a,
      LuOptions<make_real_t<S>>{.relative_pivot_tolerance = options.relative_pivot_tolerance, .retain_norms = false},
      p);
  if (attempt.factor) return slogdet_with_info(*attempt.factor);
  LogDeterminantResult<S> result{.info = attempt.info, .value = {}};
  if (attempt.info.status == SolveStatus::singular) result.value = zero_determinant<S>(p);
  return result;
}
} // namespace detail
 
template <KernelBackendSelector Backend, RankedTensorView<2> A>
  requires RealOrComplex<tensor_element_t<A>>
[[nodiscard]] auto slogdet_with_info(Backend&& backend, A const& a, Precision p,
                                     SolveOptions<make_real_t<tensor_element_t<A>>> const& options = {})
{
  return detail::slogdet_matrix_impl(std::forward<Backend>(backend), a, options, p);
}
 
template <KernelBackendSelector Backend, RankedTensorView<2> A>
  requires RealOrComplex<tensor_element_t<A>>
[[nodiscard]] auto slogdet_with_info(Backend&& backend, A const& a,
                                     SolveOptions<make_real_t<tensor_element_t<A>>> const& options = {})
{
  return detail::slogdet_matrix_impl(std::forward<Backend>(backend), a, options, precision_of(a));
}
 
template <RankedTensorView<2> A, class... Options>
  requires RealOrComplex<tensor_element_t<A>>
[[nodiscard]] auto slogdet_with_info(A const& a, Options const&... options)
{
  auto backend = select_backend_for<DenseMatrix<tensor_element_t<A>>>(lu_factor_op{});
  return slogdet_with_info(backend, a, options...);
}
 
template <class... Args>
  requires requires(Args&&... args) { slogdet_with_info(std::forward<Args>(args)...); }
[[nodiscard]] auto slogdet(Args&&... args)
{
  auto result = slogdet_with_info(std::forward<Args>(args)...);
  ERROR_IF(!result.value, "logarithmic determinant could not be resolved", static_cast<int>(result.info.status));
  return std::move(*result.value);
}
} // namespace uni20::linalg