Problem with nested optimization functions

georgios-vassos1 commented 10 months ago

Issue Description

The code encounters an issue when attempting to solve a nested optimization problem using CVXR. While the initial convex minimization problem is solved successfully, extending the optimization to include the sum of the output from the first minimization and an auxiliary function leads to code failure. Despite the convex nature of the objective function in the second optimization problem, visually confirmed through representative instances, a bug prevents the successful resolution of the nested optimization.

Code Snippet

library(CVXR)

Cost <- function(w, x, z, solver = "ECOS", ...) {
  p <- length(w)
  q <- Variable(p)

  objective   <- Minimize(t(w) %*% q)

  constraints <- list(
    q >= 0,
    q <= x,
    matrix(1.0, ncol = p) %*% q == z)

  problem <- Problem(objective, constraints)

  solve(problem, solver)$value
}

# Auxiliary component 
f <- function(z) (z-10.0)^2.0

# Manual optimization
graphical.opt <- function(w, x, solver = "ECOS", plt = TRUE, ...) {

  total.x  <- sum(x)
  span     <- seq(0.0, total.x)
  obj.vals <- sapply(span, function(z) (Cost(w, x, z, solver) + f(z)))
  opt.key  <- which.min(obj.vals)
  opt.z    <- span[opt.key]
  opt.v    <- obj.vals[opt.key]

  if (plt) {
    # Visual verification
    plot(span, obj.vals, ylab = "Objective", xlab = "Order Quantity", type = 'o', lwd = 2L)
    points(opt.z, opt.v, lwd = 5L, col = 2L)
    abline(v = opt.z, col = 2L, lwd = 2L)
    abline(h = opt.v, col = 2L, lwd = 2L)
  }

  list(
    "opt.z" = opt.z,
    "opt.v" = opt.v
  )
}

# CVXR model
cvxr.opt <- function(w, x, solver = "ECOS") {
  z <- Variable(1L)
  objective <- Minimize(Cost(w, x, z, solver) + f(z))
  constraints <- list(
    z >= 0.0,
    z <= sum(x))
  problem <- Problem(objective, constraints)

  # Solution
  result <- solve(problem, solver = solver)

  list(
    "opt.z" = result$getValue(z),
    "opt.v" = result$value
  )
}

## Problem instance
w <- matrix(c(2.0, 1.0, 3.0, 5.0), ncol = 1L)
x <- c(5.0, 5.0, 5.0, 5.0)

## Reported issue
manual <- graphical.opt(w, x)
cvx <- cvxr.opt(w, x)

# The output from manual and CVXR-based optimization modules should have been equal.
manual$opt.z == cvx$opt.z  # Should be TRUE
manual$opt.v == cvx$opt.v  # Should be TRUE

# Evaluating the objective at the optimal decision returned from solving the CVXR model should have been
# equal to the optimal value returned by solving the same model
Cost(w, x, cvx$opt.z) + f(cvx$opt.z) == cvx$opt.v # Should be TRUE

Example with Output

Running the problem instance in my code above, if you print manual$opt.z, it will return the correct value for the optimal order (manual$opt.z) and the correct minimum value of the objective (manual$opt.v).

manual$opt.z [1] 9

manual$opt.v [1] 14

You can visually verify these are the correct values with the help of the plot produced by graphical.opt when you set plt=TRUE.

To ensure the correctness of cvxr.opt, I compare its results with those of graphical.opt. However, there is a discrepancy, as the outputs from cvx.opt are:

cvx$opt.z [1] 10

cvx$opt.v [1] 4.454504e-10

The optimal ordering quantity is incorrect as it should have been 9 instead of 10, and same goes for the optimal value of objective that is 4.454504e-10 when it should have been 14. Moreover, when I evaluate the objective function Cost(w, x, cvx$opt.z) + f(cvx$opt.z) at cvx$opt.z I get the value 15 which is different from cvx$opt.v.

Version

R version 4.2.3 (2023-03-15)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS 14.3

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] CVXR_1.0-11

loaded via a namespace (and not attached):
 [1] gmp_0.7-3             Rcpp_1.0.12           compiler_4.2.3        pillar_1.9.0          nloptr_2.0.3          tools_4.2.3          
 [7] bit_4.0.5             boot_1.3-28.1         lme4_1.1-35.1         jsonlite_1.8.8        lifecycle_1.0.4       tibble_3.2.1         
[13] nlme_3.1-164          lattice_0.22-5        pkgconfig_2.0.3       rlang_1.1.3           Matrix_1.6-4          igraph_1.6.0         
[19] gurobi_10.0-2         cli_3.6.2             rstudioapi_0.15.0     microbenchmark_1.4.10 mvtnorm_1.2-4         ECOSolveR_0.5.5      
[25] Rmpfr_0.9-4           dplyr_1.1.4           systemfonts_1.0.5     generics_0.1.3        vctrs_0.6.5           bit64_4.0.5          
[31] grid_4.2.3            tidyselect_1.2.0      glue_1.7.0            R6_2.5.1              textshaping_0.3.7     fansi_1.0.6          
[37] carData_3.0-5         minqa_1.2.6           purrr_1.0.2           tidyr_1.3.0           car_3.1-2             magrittr_2.0.3       
[43] MASS_7.3-60           splines_4.2.3         abind_1.4-5           ragg_1.2.7            utf8_1.2.4            slam_0.1-50          
[49] truncnorm_1.0-9

bnaras commented 10 months ago

Thank you for the runnable example. However, I don't see a bug: you are solving a different problem in each case. For example, see what happens when you invoke:

> Cost(w, x, Variable(1))
[1] 4.807175e-11

CVXR correctly zeroes out the variables and returns the z that minimizes f(z). In your graphical.opt, you are providing the z value, which affects the last constraint in Cost. Concretely, the call to Cost in your cvxr.opt has two CVXR unknown variables, whereas your graphical solution has only one. So what problem are you trying to solve? Perhaps a mathematical formulation would clarify everything.

georgios-vassos1 commented 10 months ago

Problem Statement

Please be aware that the provided problem statement is intentionally simplified and does not fully represent the complexity of my actual problem. Its purpose is solely to illustrate the error I am encountering. While this simplified example does not encompass the entirety of my main issue, the error persists in this straightforward illustration. This simplified problem is implemented in the code I provided.

Given Parameters:

$p$: The number of capacitated sources.
$x_i$: The capacity reserved with source $i=1,\dots,p$.
$w_i$: The execution cost associated with exercising $q_i$ units from source $i=1,\dots,p$.

Objective:

Identify the optimal quantity $z\ge0$ to order from the $p$ sources.

Inner Optimization: Linear Program to Compute Cost of ordering $z$ units.

This computation involves solving a linear program defined as:

$$C(z)\doteq C(z;w,x)\doteq\min_{q_1,...,qp}\sum{i=1}^{p}w{i}q{i} $$

subject to the contraints: $$\sum{i=1}^{p}q{i}=z,\quad\text{and}\quad0 \leq q{i} \leq x{i} \ \text{for all} \ i=1,...,p$$

Outer Optimization: Optimize order quantity z.

Suppose $f$ is a known auxiliary convex function.

$$\min_{z}[C(z) + f(z)]$$

subject to the constraints:

$$0 \leq z \leq \sum{i=1}^{p}x{i}$$

In summary, the problem involves minimizing the cost associated with ordering from $p$ sources, considering the capacity constraints and execution costs. The overall objective is to find the optimal order quantity $z$. The computation involves solving a linear program to determine the cost for each feasible order quantity and evaluating the auxiliary convex function $f(z)$ with the additional constraints on $z$.

Concise Issue Description

Failure to Generate Correct Results in Nested Optimization Scenario

The code fails to generate accurate results when handling a nested optimization situation as detailed in the previous comment. The primary issue arises from utilizing the decision variable of the outer optimization program z to constrain computations within the inner optimization program, particularly in computing allocation costs.

I am new to CVXR and perhaps I am overlooking the correct approach to model this problem.

Concise Illustration of the Code

library(CVXR)

# Problem instance
p <- 4L
w <- matrix(c(2.0, 1.0, 3.0, 5.0), ncol = 1L)
x <- c(5.0, 5.0, 5.0, 5.0)

# Auxiliary component 
f <- function(z) (z-10.0)^2.0

# Decision variables
q <- Variable(p)
z <- Variable(1L)

# CVXR model
problem <- Problem(
  Minimize(
    f(z) + solve(
      Problem(
        Minimize(t(w) %*% q), 
        list(
          q >= 0.0, 
          q <= x,
          matrix(1.0, ncol = p) %*% q == z
        )
      ), 
      solver = "ECOS"
    )$value
  ), 
  list(
    z >= 0.0,
    z <= sum(x)
  )
)

# Solution inspection
result <- solve(problem, solver = "ECOS")

result$getValue(z) [1] 10 result$value [1] 4.454504e-10

I have previously clarified in my comment that these values are incorrect.

Actual Complication

The crux of the matter lies in how we integrate the decision variable from the outer optimization program into the inner optimization program within this nested structure. In such a scenario, the inner program lacks a mechanism to discern that this decision variable is not meant for optimization but rather should be treated as a fixed value, solely governed by the outer optimization.

cvxgrp / CVXR