This notebook demonstrates the solution I posted to a regression problem posted on OR Stack Exchange: “Linearizing a program with multinomial logit in the objective”. Given \(K\) observations \(y_1,\dots,y_K\) of the dependent variable, the mission is to choose nonnegative values for \(x_{i,k}\) (\(i=1,\dots,N\), \(k=1,\dots,K\)) so as to minimize \[\sum_{k=1}^K \left|y_k - \sum_{i=1}^N \frac{e^{x_{i,k}}}{\sum_{j=1}^K e^{x_{i,j}}}\right|.\]
We will generate a test problem and solve it, using OMPR + ROI + CPLEX to solve the modified model (a linear program).
library(ROI) # general solver interfaceROI: R Optimization Infrastructure
Registered solver plugins: nlminb, cplex.
Default solver: auto.library(ROI.plugin.cplex) # connects the ROI package to CPLEX
library(ompr) # MIP modeling packageRegistered S3 method overwritten by 'data.table':
method from
print.data.table library(ompr.roi) # connects OMPR to ROI
library(magrittr) # for the pipe operator
library(GA) # to compare a genetic algorithm with the LPLoading required package: foreach
Loading required package: iterators
____ _
/ ___| / \ Genetic
| | _ / _ \ Algorithms
| |_| |/ ___ \
\____/_/ \_\ version 3.2
Type 'citation("GA")' for citing this R package in publications.
Attaching package: ‘GA’
The following object is masked from ‘package:utils’:
deWe will pick some arbitrary dimensions and generate “dependent variable observations” using a uniform distribution. (The purpose here is to illustrate the computations, not to fit actual data.)
N <- 10 # number of "agents"
K <- 100 # number of sample points
set.seed(123) # for reproducibility
y <- runif(K) # sample dataWe will also code the objective function as specified in the original question, so that we can check any proposed solution.
obj <- function(x) {
# The input x will be an N x K matrix.
# We convert each x[i, j] to e^x[i, j].
xx <- exp(x)
# Next, we scale by row sums.
xx <- xx / apply(xx, 1, sum)
# Finally, we subtract column sums from y, take absolute values, and sum.
(y - apply(xx, 2, sum)) %>% abs() %>% sum()
}The proposed model involves a change of variables to \[z_{i,k}=\frac{e^{x_{i,k}}}{\sum_{j=1}^K e^{x_{i,j}}}.\] The \(z\) variables are nonnegative and must satisfy one constraint: \[\sum_{k=1}^{K}z_{i, k}=1\quad\forall i=1,\dots,N.\]The objective becomes \[\min \sum_{k=1}^K \left|y_k - \sum_{i=1}^N z_{i,k}\right|.\] We linearize this by adding new variables \(w_k\ge 0\), whose sum we minimize, together with the constraints \[-w_k \le y_k - \sum_{i=1}^N z_{i,k} \le w_k \quad \forall k=1,\dots,K.\]
Since the matrix of \(z\) variables has to be flattened into a vector for OMPR, we need indexing functions to convert between one- and two-dimensional representations.
ix1 <- function(i, k) {
# Convert i in 1...N and k in 1...K to a single index in 1...N*K.
(i - 1) * K + k
}
ix2 <- function(t) {
# Convert t in 1...N*K to vector (i,k) with i in 1...N and k in 1...K
c((t - 1) %/% K + 1, (t - 1) %% K + 1)
}Now we build the LP model. It requires a strictly positive lower bound on \(z\), which I will arbitrarily set.
zlo <- 1e-7model <-
# Create an empty model.
MIPModel() %>%
# Add the z and w variables. z.var[t] is z[ix2(t)].
add_variable(z.var[j], j = 1:(N * K), type = "continuous", lb = zlo) %>%
add_variable(w.var[k], k = 1:K, type = "continuous", lb = 0) %>%
# Set the objective to minimize the sum of the w variables.
set_objective(sum_expr(w.var[k], k = 1:K), "min") %>%
# Constrain each row of the z matrix to sum to 1.
add_constraint(sum_expr(z.var[j], j = ix1(i, 1):ix1(i, K)) == 1, i = 1:N) %>%
# Constrain w[k] to be at least the absolute error in observation k.
add_constraint(sum_expr(z.var[j], j = 1:(N * K), (j - 1) %% K + 1 == k) - w.var[k] <= y[k], k = 1:K) %>%
add_constraint(sum_expr(z.var[j], j = 1:(N * K), (j - 1) %% K + 1 == k) + w.var[k] >= y[k], k = 1:K)
adding constraints [================================>--------] 80% eta 0s
adding constraints [====================================>----] 90% eta 0s
adding constraints [=========================================] 100% eta 0s
adding constraints [============>----------------------------] 31% eta 0s
adding constraints [============>----------------------------] 32% eta 0s
adding constraints [=============>---------------------------] 33% eta 0s
adding constraints [=============>---------------------------] 34% eta 0s
adding constraints [=============>---------------------------] 35% eta 0s
adding constraints [==============>--------------------------] 36% eta 0s
adding constraints [==============>--------------------------] 37% eta 0s
adding constraints [===============>-------------------------] 38% eta 0s
adding constraints [===============>-------------------------] 39% eta 0s
adding constraints [===============>-------------------------] 40% eta 0s
adding constraints [================>------------------------] 41% eta 0s
adding constraints [================>------------------------] 42% eta 0s
adding constraints [=================>-----------------------] 43% eta 0s
adding constraints [=================>-----------------------] 44% eta 0s
adding constraints [=================>-----------------------] 45% eta 0s
adding constraints [==================>----------------------] 46% eta 0s
adding constraints [==================>----------------------] 47% eta 0s
adding constraints [===================>---------------------] 48% eta 0s
adding constraints [===================>---------------------] 49% eta 0s
adding constraints [===================>---------------------] 50% eta 0s
adding constraints [====================>--------------------] 51% eta 0s
adding constraints [====================>--------------------] 52% eta 0s
adding constraints [=====================>-------------------] 53% eta 0s
adding constraints [=====================>-------------------] 54% eta 0s
adding constraints [======================>------------------] 55% eta 0s
adding constraints [======================>------------------] 56% eta 0s
adding constraints [======================>------------------] 57% eta 0s
adding constraints [=======================>-----------------] 58% eta 0s
adding constraints [=======================>-----------------] 59% eta 0s
adding constraints [========================>----------------] 60% eta 0s
adding constraints [========================>----------------] 61% eta 0s
adding constraints [========================>----------------] 62% eta 0s
adding constraints [=========================>---------------] 63% eta 0s
adding constraints [=========================>---------------] 64% eta 0s
adding constraints [==========================>--------------] 65% eta 0s
adding constraints [==========================>--------------] 66% eta 0s
adding constraints [==========================>--------------] 67% eta 0s
adding constraints [===========================>-------------] 68% eta 0s
adding constraints [===========================>-------------] 69% eta 0s
adding constraints [============================>------------] 70% eta 0s
adding constraints [============================>------------] 71% eta 0s
adding constraints [=============================>-----------] 72% eta 0s
adding constraints [=============================>-----------] 73% eta 0s
adding constraints [=============================>-----------] 74% eta 0s
adding constraints [==============================>----------] 75% eta 0s
adding constraints [==============================>----------] 76% eta 0s
adding constraints [===============================>---------] 77% eta 0s
adding constraints [===============================>---------] 78% eta 0s
adding constraints [===============================>---------] 79% eta 0s
adding constraints [================================>--------] 80% eta 0s
adding constraints [================================>--------] 81% eta 0s
adding constraints [=================================>-------] 82% eta 0s
adding constraints [=================================>-------] 83% eta 0s
adding constraints [=================================>-------] 84% eta 0s
adding constraints [==================================>------] 85% eta 0s
adding constraints [==================================>------] 86% eta 0s
adding constraints [===================================>-----] 87% eta 0s
adding constraints [===================================>-----] 88% eta 0s
adding constraints [===================================>-----] 89% eta 0s
adding constraints [====================================>----] 90% eta 0s
adding constraints [====================================>----] 91% eta 0s
adding constraints [=====================================>---] 92% eta 0s
adding constraints [=====================================>---] 93% eta 0s
adding constraints [======================================>--] 94% eta 0s
adding constraints [======================================>--] 95% eta 0s
adding constraints [======================================>--] 96% eta 0s
adding constraints [=======================================>-] 97% eta 0s
adding constraints [=======================================>-] 98% eta 0s
adding constraints [========================================>] 99% eta 0s
adding constraints [=========================================] 100% eta 0s
adding constraints [=============>---------------------------] 33% eta 0s
adding constraints [=============>---------------------------] 34% eta 0s
adding constraints [=============>---------------------------] 35% eta 0s
adding constraints [==============>--------------------------] 36% eta 0s
adding constraints [==============>--------------------------] 37% eta 0s
adding constraints [===============>-------------------------] 38% eta 0s
adding constraints [===============>-------------------------] 39% eta 0s
adding constraints [===============>-------------------------] 40% eta 0s
adding constraints [================>------------------------] 41% eta 0s
adding constraints [================>------------------------] 42% eta 0s
adding constraints [=================>-----------------------] 43% eta 0s
adding constraints [=================>-----------------------] 44% eta 0s
adding constraints [=================>-----------------------] 45% eta 0s
adding constraints [==================>----------------------] 46% eta 0s
adding constraints [==================>----------------------] 47% eta 0s
adding constraints [===================>---------------------] 48% eta 0s
adding constraints [===================>---------------------] 49% eta 0s
adding constraints [===================>---------------------] 50% eta 0s
adding constraints [====================>--------------------] 51% eta 0s
adding constraints [====================>--------------------] 52% eta 0s
adding constraints [=====================>-------------------] 53% eta 0s
adding constraints [=====================>-------------------] 54% eta 0s
adding constraints [======================>------------------] 55% eta 0s
adding constraints [======================>------------------] 56% eta 0s
adding constraints [======================>------------------] 57% eta 0s
adding constraints [=======================>-----------------] 58% eta 0s
adding constraints [=======================>-----------------] 59% eta 0s
adding constraints [========================>----------------] 60% eta 0s
adding constraints [========================>----------------] 61% eta 0s
adding constraints [========================>----------------] 62% eta 0s
adding constraints [=========================>---------------] 63% eta 0s
adding constraints [=========================>---------------] 64% eta 0s
adding constraints [==========================>--------------] 65% eta 0s
adding constraints [==========================>--------------] 66% eta 0s
adding constraints [==========================>--------------] 67% eta 0s
adding constraints [===========================>-------------] 68% eta 0s
adding constraints [===========================>-------------] 69% eta 0s
adding constraints [============================>------------] 70% eta 0s
adding constraints [============================>------------] 71% eta 0s
adding constraints [=============================>-----------] 72% eta 0s
adding constraints [=============================>-----------] 73% eta 0s
adding constraints [=============================>-----------] 74% eta 0s
adding constraints [==============================>----------] 75% eta 0s
adding constraints [==============================>----------] 76% eta 0s
adding constraints [===============================>---------] 77% eta 0s
adding constraints [===============================>---------] 78% eta 0s
adding constraints [===============================>---------] 79% eta 0s
adding constraints [================================>--------] 80% eta 0s
adding constraints [================================>--------] 81% eta 0s
adding constraints [=================================>-------] 82% eta 0s
adding constraints [=================================>-------] 83% eta 0s
adding constraints [=================================>-------] 84% eta 0s
adding constraints [==================================>------] 85% eta 0s
adding constraints [==================================>------] 86% eta 0s
adding constraints [===================================>-----] 87% eta 0s
adding constraints [===================================>-----] 88% eta 0s
adding constraints [===================================>-----] 89% eta 0s
adding constraints [====================================>----] 90% eta 0s
adding constraints [====================================>----] 91% eta 0s
adding constraints [=====================================>---] 92% eta 0s
adding constraints [=====================================>---] 93% eta 0s
adding constraints [======================================>--] 94% eta 0s
adding constraints [======================================>--] 95% eta 0s
adding constraints [======================================>--] 96% eta 0s
adding constraints [=======================================>-] 97% eta 0s
adding constraints [=======================================>-] 98% eta 0s
adding constraints [========================================>] 99% eta 0s
adding constraints [=========================================] 100% eta 0s
We solve the LP model and extract the solution.
# Solve the LP.
solution <- model %>% solve_model(with_ROI(solver = "cplex"))CPLEX environment opened
Closed CPLEX environment# Get and print the objective value.
objval <- solution$objective_value
cat("Optimal objective value = ", objval, "\n")Optimal objective value = 39.8559 We can extract \(w\) and \(z\) from the solver output.
# The solution contains the values of w first, followed by the values of z.
w <- solution$solution[grepl("w.var*", names(solution$solution))]
z <- matrix(solution$solution[grepl("z.var*", names(solution$solution))], nrow = N, byrow = T)
The process to recover \(x\) from \(z\) is explained in the solution I posted to the original question. For each \(i=1,\dots,N\), we find the index \(j\) of the smallest entry \(z_{i,k}\) in that row, then set \(x_{i,k} = \log(z_{i,k} / z_{i,j})\) for \(k=1,\dots,K\).
# This function converts a row of the z matrix to a row of the x matrix.
z.to.x <- function(zz) {
# zz is one row of z.
log(zz / min(zz))
}
# We apply the function to each row of z to get x.
x <- t(apply(z, 1, z.to.x))As a sanity check, we confirm the objective value of the proposed solution.
obj(x)[1] 39.8559As a further sanity check, we try to minimize the original objective function using a genetic algorithm. A GA solution with a larger objective value than that of the LP will not prove that the LP solution is optimal, but a GA solution with a smaller objective value than that of the LP will definitely indicate a problem.
The “chromosome” in the GA will just be the \(x\) matrix flattened to a double precision vector, the same way that we flattened \(z\) in the LP model. The only constraint in the original formulation was the requirement that \(x\ge 0\), which we enforce as a lower bound. To compute the fitness of a chromosome, we convert the vector of values to a matrix and apply the obj() function to it. Because the GA package maximizes fitness and we are minimizing, we negate the original objective function.
fitness.fun <- function(chromosome) {
temp <-
chromosome %>%
# Reshape it to an N x K matrix.
matrix(nrow = N, byrow = T) %>%
# Apply the objective function to the matrix.
obj()
# Return the negated value.
-temp
}We build and run the genetic algorithm (using mostly default settings). The ga() function requires upper bounds for the \(x\) variables, which we set rather arbitrarily.
genalg <- ga(
type = "real-valued",
fitness = fitness.fun,
lower = rep.int(0, N * K),
upper = rep.int(100, N * K),
monitor = FALSE,
maxiter = 1000,
run = 100,
seed = 123
)
cat("The best fit from the GA has error ", -genalg@fitnessValue, " compared to the LP result ", obj(x), ".\n")The best fit from the GA has error 39.8559 compared to the LP result 39.8559 .It appears that the GA was able to match the LP result but not improve on it.
Let’s define a function to calculate the density (fraction of elements that are nonzero) in a matrix or vector.
fractionDense <- function(x) {
temp <-
x %>%
# Flatten the argument to a vector.
as.vector() %>%
# Truncate small entries to zero.
zapsmall()
sum(temp != 0) / length(temp)
}We can compare the density of the LP and GA solutions.
cat("The LP solution has density ", fractionDense(x), ".\n")The LP solution has density 0.034 .cat("The GA solution has density ", fractionDense(genalg@solution[1,]), ".\n")The GA solution has density 1 .If we want to use the LP model but desire a denser solution, then (using the invariance property mentioned in my solution online) we can add an arbitrary positive amount \(\lambda_i\) to each \(x_{i,k}\). We do that next, confirm the density, and check the objective value to confirm that it is unchanged.
lambda <- runif(N)
x2 <- apply(x, 2, function(t) t + lambda)
cat("The adjusted solution x has density ", fractionDense(x2), " and objective value ", obj(x2), ".\n")The adjusted solution x has density 1 and objective value 39.8559 .Out of curiosity, let’s see whether a similar transformation can make the GA solution very sparse.
# Get the GA solution and convert it to a matrix.
x3 <- genalg@solution[1,] %>% matrix(nrow = N, ncol = K, byrow = T)
# Find the smallest $x$ value for each $i$ (which we can subtract and still preserve nonnegativity).
x4 <- apply(x3, 1, min)
# Subtract each row min from that row.
x3 <- apply(x3, 2, function(t) t - x4)
# Confirm that the revised solution is nonnegative.
cat("The minimum of the modified solution is ", min(x3), ".\n")The minimum of the modified solution is 0 .# Confirm that the objective value is unchanged.
cat("The modified solution has objective value ", obj(x3), ".\n")The modified solution has objective value 39.8559 .# Check its density.
cat("The density of the modified solution is ", fractionDense(x3), ".\n")The density of the modified solution is 0.99 .So we reduced the density, but only slightly.