ops/svd_least_squares.hpp

Go to the documentation of this file.

#pragma once
 
#include <uni20/linalg/backends/lapack/svd_least_squares.hpp>
#include <uni20/linalg/precision.hpp>
#include <uni20/tensor/precision.hpp>
 
namespace uni20::linalg
{
template <RealOrComplex S> struct SvdLeastSquaresResult
{
    LeastSquaresInfo info = {};
    std::optional<DenseMatrix<S>> solution = {};
    std::vector<make_real_t<S>> residual_norms = {};
    std::optional<SvdLeastSquaresDiagnostics<make_real_t<S>>> diagnostics = {};
};
 
template <KernelBackendSelector Backend, RankedTensorView<2> A, RankedTensorView<2> B>
  requires RealOrComplex<tensor_element_t<A>> && std::same_as<tensor_element_t<A>, tensor_element_t<B>>
[[nodiscard]] auto
svd_least_squares_with_info(Backend&& backend, A const& a, B const& b, Precision p,
                            SvdLeastSquaresOptions<make_real_t<tensor_element_t<A>>> const& options = {})
{
  using S = tensor_element_t<A>;
  require_precision<S>(p);
  (void)p.bit_count();
  ERROR_IF(a.extent(0) != b.extent(0), "SVD least squares requires matching RHS rows");
  least_squares_detail::require_svd_options(options);
  SvdLeastSquaresResult<S> result;
  result.solution.emplace(uninitialized, a.extent(1), b.extent(1), p);
  result.residual_norms.resize(b.extent(1));
  result.diagnostics.emplace();
  auto ad = mdspec_of(a);
  auto bd = mdspec_of(b);
  auto xd = mdspec_of(*result.solution);
  detail::dispatch_at_precision<S>(std::forward<Backend>(backend), svd_least_squares_op{}, p, xd, ad, bd,
                                   std::span(result.residual_norms), result.info, *result.diagnostics, options);
  if (!result.info.succeeded())
  {
    result.solution.reset();
    result.residual_norms.clear();
    result.diagnostics.reset();
  }
  return result;
}
 
template <KernelBackendSelector Backend, RankedTensorView<2> A, RankedTensorView<2> B>
  requires RealOrComplex<tensor_element_t<A>> && std::same_as<tensor_element_t<A>, tensor_element_t<B>>
[[nodiscard]] auto
svd_least_squares_with_info(Backend&& backend, A const& a, B const& b,
                            SvdLeastSquaresOptions<make_real_t<tensor_element_t<A>>> const& options = {})
{
  return svd_least_squares_with_info(std::forward<Backend>(backend), a, b, common_default_precision(a, b), options);
}
 
template <RankedTensorView<2> A, RankedTensorView<2> B>
  requires RealOrComplex<tensor_element_t<A>> && std::same_as<tensor_element_t<A>, tensor_element_t<B>>
[[nodiscard]] auto
svd_least_squares_with_info(A const& a, B const& b, Precision p,
                            SvdLeastSquaresOptions<make_real_t<tensor_element_t<A>>> const& options = {})
{
  return svd_least_squares_with_info(select_backend_for<DenseMatrix<tensor_element_t<A>>>(svd_least_squares_op{}), a, b,
                                     p, options);
}
 
template <RankedTensorView<2> A, RankedTensorView<2> B>
  requires RealOrComplex<tensor_element_t<A>> && std::same_as<tensor_element_t<A>, tensor_element_t<B>>
[[nodiscard]] auto
svd_least_squares_with_info(A const& a, B const& b,
                            SvdLeastSquaresOptions<make_real_t<tensor_element_t<A>>> const& options = {})
{
  return svd_least_squares_with_info(a, b, common_default_precision(a, b), options);
}
 
template <class... Args>
  requires requires(Args&&... args) { svd_least_squares_with_info(std::forward<Args>(args)...); }
[[nodiscard]] auto svd_least_squares(Args&&... args)
{
  auto result = svd_least_squares_with_info(std::forward<Args>(args)...);
  ERROR_IF(!result.info.succeeded(), "SVD least squares failed", static_cast<int>(result.info.status));
  return std::move(*result.solution);
}
} // namespace uni20::linalg