lmi.hpp

Go to the documentation of this file.

/*
 * Copyright © 2025-2026 Akreon LLC
 * SPDX-License-Identifier: Apache-2.0
 */
 
#pragma once
#include "akreon/convex/cones/psd_cone.hpp"
 
#include "convex_set.hpp"
 
namespace convex {
 
template <linalg::RnSpace Space>
using LMIDom = ConvexDomain<Space, VarType::FREE>;
 
template <linalg::RnSpace Space>
using LMICodom = ConvexDomain<Space, VarType::FREE>;
 
template <linalg::RnSpace Space>
using LMIMeta =
    ConvexAtomMeta<LMIDom<Space>, LMICodom<Space>, CurvatureType::CONVEX>;
 
template <typename Meta, typename Child, linalg::SquareMatrix MatIn>
struct LMI : ConvexSet<Meta, LMI<Meta, Child, MatIn>, Child> {
  using ConvexBase = ConvexSet<Meta, LMI<Meta, Child, MatIn>, Child>;
  using typename ConvexBase::PointT;
  using typename ConvexBase::T;
 
  using ConvexBase::N;
 
  using Mat = linalg::remove_cvref_t<MatIn>;
  static constexpr std::size_t MatDim = linalg::mat_rows_v<Mat>;
  static_assert(linalg::mat_cols_v<Mat> == MatDim,
                "LMI expects square coefficient matrices");
 
  constexpr explicit LMI(const linalg::Vector<Mat, N>& A_in, const Mat& B_in,
                         const Child& child_in)
      : ConvexBase(child_in), A(A_in), B(B_in) {}
 
  constexpr bool contains_impl(const PointT& x, T /* tol */) const {
    auto res = B;
    for (std::size_t i = 0; i < N; ++i) res = res + (x(i) * A(i));
    linalg::LDL ldl(res);
    return ldl.is_psd();
  }
 
  constexpr bool is_well_formed() const { return true; }
 
  template <typename Other>
  constexpr bool is_equal_impl(const Other& other) const {
    if (!linalg::all_equal(B, other.B)) return false;
    for (std::size_t i = 0; i < N; ++i)
      if (!linalg::all_equal(A(i), other.A(i))) return false;
    return true;
  }
 
  template <linalg::MatrixContainer Mat>
  constexpr auto transform_impl(const affine::AffineMap<Mat>& /* map */) const {
    AK_ASSERT(false && "LMI does not support transformations");
  }
 
  template <AtomType Type, conic::BuilderLike Builder>
  constexpr auto canonicalize(const Builder& builder) const {
    const auto child_builder =
        std::get<0>(this->children)
            .template canonicalize<AtomType::SET>(builder);
    const auto child_coord = child_builder.coord_map;
 
    static constexpr std::size_t ChildCodom =
        linalg::mat_cols_v<decltype(child_coord.linear_part())>;
    static_assert(ChildCodom == N,
                  "Convex set domain does not equal the child coordinate "
                  "linear map's codomain in size");
 
    static constexpr std::size_t K = MatDim * (MatDim + 1) / 2;
    linalg::Matrix<T, K, N> F_local{};
    linalg::Vector<T, K> g_local{};
 
    auto update_F = [&F_local](const auto& i, const auto& M) {
      const auto m_svec = linalg::get_svec(M);
      for (std::size_t row_idx = 0; row_idx < K; ++row_idx)
        F_local(row_idx, i) = m_svec(row_idx);
    };
 
    auto update_g = [&g_local](const auto& M) {
      g_local = linalg::get_svec(M);
    };
 
    for (std::size_t i = 0; i < N; ++i) update_F(i, A(i));
    update_g(B);
 
    const auto F = F_local * child_coord.linear_part();
    const auto g = F_local * child_coord.translation() + g_local;
 
    const auto map = affine::AffineMap<decltype(F)>(F, g);
    const auto psd_cone =
        conic::PositiveSemidefiniteCone<linalg::Vector<T, K>, MatDim>();
    const auto cones = std::make_tuple(psd_cone);
    const auto emission =
        conic::LocalEmission(map, cones, child_builder.eq_map);
 
    return child_builder.template realize<Type>(emission);
  }
 
  linalg::Vector<Mat, N> A; 
  Mat B;                    
};
 
template <typename Mat, std::size_t N, typename Child>
LMI(const linalg::Vector<Mat, N>&, const Mat&, const Child&)
    -> LMI<LMIMeta<dom_space_t<meta_codom_t<atom_meta_t<Child>>>>,
           atom_decay_t<Child>, Mat>;
 
} // namespace convex