R/164_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))
  }

  x <- args[[1L]]
  shape <- expr@shape
  t <- Variable(shape = shape)
  constraints <- list()

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

  r <- Variable(shape = x@shape)
  ## 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 = x@shape[1L], ncol = x@shape[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) {
  x <- args[[1L]]
  axis <- expr@axis
  shape <- expr@shape
  t <- Variable(shape = shape)

  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))
  }

  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 <- expr@shape
  t <- Variable(shape = shape)

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

  r <- Variable(shape = x@shape)
  ## 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 = x@shape[1L], ncol = x@shape[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))
  }

  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