Non Linear Optimization in R with nloptr vs Excel

Question

I can't seem to find an answer for my particular problem. I'm trying to minimize the mean absolute error between a real value vector and a linear combination of models as follows:

library(nloptr)

df <- data.frame(
  real = c(24.2418, 21.7374, 7.203, 115.0233, 16.632, 5.4644, 27.8917, 0.5904, 0.633, 105.3656, 110.0564, 122.9399, 23.0418, 4.2186, 109.5453), 
  m_1 = c(17.7790979854662, 32.93639375, 6.94294375000066, 98.909065625, 80.1068562499999, 11.656556250002, 39.868921875, 0.859157480988586, 0.612625482376216, 112.580383416359, 151.896765176957, 155.81521460987, 7.3257, 4.1268, 41.5711703205879), 
  m_2 = c(25.5062900474607, 32.7709877317137, 6.94317649489279, 98.898593448344, 80.1114267521597, 11.6563450007055, 39.8691722409057, 0.907260912997819, 0.602795676399197, 114.183526809465, 139.51052151724, 149.993624420536, 6.85002142907728, 4.66668862149305, 70.7527906311631), 
  m_3 = c(27.1495912199999, 40.2339353399999, 7.10339542, 87.1444967133334, 58.4563384933332, 11.1378198366666, 37.6030141333333, 0.852288459999999, 0.681724560000002, 100.101136186666, 118.536525109999, 136.923264319999, 5.64763034333333, 4.8659515, 70.12675835), 
  m_4 = c(25.511590625, 32.9363937499999, 7.00050769504031, 98.3660974929738, 80.10685625, 11.65655625, 39.868921875, 0.665580984791989, 0.498756215272925, 85.791042265746, 135.619508469251, 140.946144588455, 5.05824305930683, 3.25333636845094, 22.0908752674237), 
  m_5 = c(25.6118152340655, 34.5646719234769, 6.82416840383483, 91.5582383465651, 84.4603847826215, 11.3405620277701, 40.7906062254877, 0.908706704665592, 0.602817399156822, 114.326905157898, 139.595783699511, 150.046375909198, 6.8511793011574, 4.6622942290559, 56.2753212961812), 
  m_6 = c(21.9868574376585, 44.3147731773466, 6.38315486686481, 100.303757097094, 9.13921010739697, 7.83817900918309, 31.5458855316741, 1.09960505333834, 0.817751834425678, 101.110814846224, 145.55847538105, 142.82362305075, 7.61732986965459, 4.6774198307473, 67.5821464371521)
)

best_dist <- function(x) {
  output <- df$m_1 * x[1] + df$m_2 * x[2] + df$m_3 * x[3] + 
    df$m_4 * x[4] + df$m_5 * x[5] + df$m_6 * x[6]
  mean(abs(output - df$real))
}

restriction <- function(x) sum(x) - 1

nloptr(
  x0 = rep(1 / 6, 6), 
  eval_f = best_dist, 
  lb = rep(0, 6), 
  ub = rep(1, 6), 
  eval_g_eq = restriction, 
  opts = list(algorithm = "NLOPT_GN_ISRES", xtol_rel = 1e-16, maxeval = 1e4)
)

As you could read I'm using the nloptr package. The above code yields the not optimal result of 14.85 for the objective function and the parameters are all the inital parameters. You may change the initial parameters to some other vector and still won't get the optimal solution.

However, using the excel solver one can easily get a result of 10.77 for the objective function and (0, 0, .15, 0, 0, .85) for the parameters.

I've tried using an algorithm with gradient, however I can't seem to get the syntax right. Here's my other attempt.

gradient <- function(x) {
  output <- df$m_1 * x[1] + df$m_2 * x[2] + df$m_3 * x[3] + 
    df$m_4 * x[4] + df$m_5 * x[5] + df$m_6 * x[6]
  err <- output - df$real

  c(
    - sum(sign(err) * df$m_1), 
    - sum(sign(err) * df$m_2), 
    - sum(sign(err) * df$m_3), 
    - sum(sign(err) * df$m_4), 
    - sum(sign(err) * df$m_5), 
    - sum(sign(err) * df$m_6)
  )
}

nloptr(
  x0 = runif(6), 
  eval_f = best_dist, 
  eval_grad_f = gradient, 
  lb = rep(0, 6), 
  ub = rep(1, 6), 
  eval_g_eq = restriction, 
  opts = list(algorithm = "NLOPT_GN_ISRES", xtol_rel = 1e-16, maxeval = 1e4)
)

Answer 1

This is like LAD Regression with a side constraint. LP (linear programming) formulations for LAD regression are shown here. Using your naming of the data (real[i], m[i,j]), we have:

min sum(i, | real[i] - sum(j,m[i,j]*x[j]) |  ) 
subject to
   sum(j, x[j]) = 1
   0 <= x[j] <= 1

This can be linearized (using variable splitting) as follows:

min sum(i, r1[i]+r2[i])
subject to
   r1[i]-r2[i] = real[i] - sum(j,m[i,j]*x[j])
   sum(j, x[j]) = 1
bounds:
   r1[i], r2[i] >= 0       (positive variables)
   0 <= x[j] <= 1

This is a pure LP and can be solved with any LP solver.