R/174_reductions_dcp2cone_canonicalizers_pnorm_canon.R

In CVXR: Disciplined Convex Optimization

#####
## DO NOT EDIT THIS FILE!! EDIT THE SOURCE INSTEAD: rsrc_tree/reductions/dcp2cone/canonicalizers/pnorm_canon.R
#####

## CVXPY SOURCE: reductions/dcp2cone/canonicalizers/pnorm_canon.py
## p=2: SOC constraint; p!=2: abs + PowCone3D (exact)
## approx: Uses gm_constrs (SOC constraints via rational approximation)


# -- pnorm_exact_canon ---------------------------------------------
pnorm_exact_canon <- function(expr, args, solver_context = NULL) {
  p <- expr@p

  ## p=2: direct SOC
  if (p == 2) {
    return(.pnorm_p2_canon(expr, args, solver_context = solver_context))
  }

  x <- args[[1L]]
  shape <- .shape(expr)
  bounds <- get_expr_bounds_if_supported(expr, solver_context)
  t <- Variable(shape = shape, bounds = bounds)
  constraints <- list()

  ## For p > 1, take absolute value first
  if (p > 1) {
    abs_expr <- Abs(x)
    abs_result <- abs_canon(abs_expr, .args(abs_expr))
    x <- abs_result[[1L]]
    constraints <- c(constraints, abs_result[[2L]])
  }

  r <- Variable(shape = .shape(x))
  ## sum(r) == t
  constraints <- c(constraints, list(sum(r) == t))

  ## promoted_t = ones * t (broadcast scalar t to match x shape)
  promoted_t <- Constant(matrix(1, nrow = .shape(x)[1L], ncol = .shape(x)[2L])) * t

  if (p < 0) {
    alpha <- as.numeric(-p / (1 - p))
    constraints <- c(constraints, list(
      PowCone3D(.cvxr_vec(x), .cvxr_vec(r), .cvxr_vec(promoted_t), alpha)
    ))
  } else if (p > 0 && p < 1) {
    alpha <- as.numeric(p)
    constraints <- c(constraints, list(
      PowCone3D(.cvxr_vec(x), .cvxr_vec(promoted_t), .cvxr_vec(r), alpha)
    ))
  } else if (p > 1) {
    alpha <- as.numeric(1 / p)
    constraints <- c(constraints, list(
      PowCone3D(.cvxr_vec(r), .cvxr_vec(promoted_t), .cvxr_vec(x), alpha)
    ))
  }

  list(t, constraints)
}

## Internal: p=2 pnorm via SOC
## CVXPY SOURCE: pnorm_canon.py lines 35-45
.pnorm_p2_canon <- function(expr, args, solver_context = NULL) {
  x <- args[[1L]]
  axis <- expr@axis
  shape <- .shape(expr)
  bounds <- get_expr_bounds_if_supported(expr, solver_context)
  t <- Variable(shape = shape, bounds = bounds)

  if (is.null(axis)) {
    ## Scalar output: SOC(t, vec(x))
    list(t, list(SOC(t, .cvxr_vec(x))))
  } else {
    ## axis-aware: SOC(vec(t), x, axis)
    list(t, list(SOC(.cvxr_vec(t), x, axis)))
  }
}

# -- pnorm_approx_canon --------------------------------------------
## CVXPY SOURCE: dcp2cone/canonicalizers/pnorm_canon.py lines 90-139
## Uses gm_constrs (SOC constraints via rational approximation)
pnorm_approx_canon <- function(expr, args, solver_context = NULL) {
  p <- expr@p

  ## p=2: direct SOC (shared with exact)
  if (p == 2) {
    return(.pnorm_p2_canon(expr, args, solver_context = solver_context))
  }

  x <- args[[1L]]
  ## Convert to exact rational for weight computation
  ## NOTE: as.bigq(float) gives huge denominators (e.g., as.bigq(1.6) -> 3.6e15/2.3e15).
  ## Use .limit_denominator() to recover clean rational, matching pow_high/mid/neg pattern.
  p_bigq <- .limit_denominator(as.bigq(p), as.bigz(expr@max_denom))
  shape <- .shape(expr)
  bounds <- get_expr_bounds_if_supported(expr, solver_context)
  t <- Variable(shape = shape, bounds = bounds)

  constraints <- list()
  if (p > 1) {
    abs_expr <- Abs(x)
    abs_result <- abs_canon(abs_expr, .args(abs_expr))
    x <- abs_result[[1L]]
    constraints <- c(constraints, abs_result[[2L]])
  }

  r <- Variable(shape = .shape(x))
  ## sum(r) == t
  constraints <- c(constraints, list(sum(r) == t))

  ## promoted_t = ones * t (broadcast scalar t to match x shape)
  promoted_t <- Constant(matrix(1, nrow = .shape(x)[1L], ncol = .shape(x)[2L])) * t

  if (p < 0) {
    w <- c(-p_bigq / (as.bigq(1L) - p_bigq), as.bigq(1L) / (as.bigq(1L) - p_bigq))
    constraints <- c(constraints, gm_constrs(promoted_t, list(x, r), w))
  } else if (p > 0 && p < 1) {
    w <- c(p_bigq, as.bigq(1L) - p_bigq)
    constraints <- c(constraints, gm_constrs(r, list(x, promoted_t), w))
  } else if (p > 1) {
    w <- c(as.bigq(1L) / p_bigq, as.bigq(1L) - as.bigq(1L) / p_bigq)
    constraints <- c(constraints, gm_constrs(x, list(r, promoted_t), w))
  }

  ## CVXPY SOURCE: pnorm_canon.py:121-137 -- filters constraints to SOCs
  ## only (the tower also contains an equality `sum(r) == t` and the abs
  ## decomposition, which shouldn't count toward num_soc).
  num_soc <- sum(vapply(constraints, function(c) .s7_is(c, SOC), logical(1L)))
  .maybe_warn_soc_approx(
    solver_context, target_cone = PowCone3D,
    approx_error = expr@approx_error, num_soc = num_soc,
    msg_head = sprintf("pnorm with p=%g is being approximated",
                       as.numeric(expr@original_p))
  )

  list(t, constraints)
}

method(dcp_canonicalize, Pnorm) <- pnorm_exact_canon
method(has_dcp_canon, Pnorm) <- function(expr) TRUE
method(dcp_canonicalize, PnormApprox) <- pnorm_approx_canon
method(has_dcp_canon, PnormApprox) <- function(expr) TRUE