LBFGS Optimization to maximize function with inequality constraints on sum

Modeling
, , ,
1

I’m trying to build a stan program that optimizes a non-linear function using the $optimize function. I made a simplified example. The program works correctly and it optimizes to the max value, the problem is that it breaks when I try to add some inequality constraints on the sum, if the max value I specify for the constraint is more than the maximum possible value it works, it just fails to find out a value with a lower constraint.

data{
  // shape
  int<lower=0> N;
  
  // data
  real<lower=0> mult;
  vector[N] x;
  
  // constraints
  real<lower=0> max_value;
}

parameters{
  vector<lower=0, upper=3>[N] rates;
}

transformed parameters{
  
  vector<lower=0>[N] x_rate = x .* rates; 
  vector<lower=0>[N] x_rate_mult = x_rate .* mult;
  
  real<lower=0, upper=max_value> total = sum(x_rate_mult);
}

model{
  total ~ uniform(0, max_value);
  
  
  target += x_rate_mult;
  
}

I just generated some random data and calculated the maximum posible value.

N <- 200
data_list <- list(
  N = N,
  x = rpois(200, 1000),
  mult = 1.5,
  max_value = 10000000
)

print(str_c("Max possible value ", sum(data_list$x * data_list$mult * 3)))

If I specify the max value to the one I precalculated or something bigger it works

data_list$max_value <- 802277

fit <- optim_model$optimize(
  data = data_list
)

print(str_c("Total ", sum(fit$summary("total")$estimate)))

If max_value is smaller then it fails like this

Initial log joint probability = 30075.9 
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
       1       30075.9             0       2056.77       0.001       0.001       35    
Exception: optim_test_model_namespace::log_prob: total is 902263, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 893353, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 902277, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 902277, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 902271, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 901093, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 886001, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 844199, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 902277, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 902277, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 902276, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 901985, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 896404, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 876403, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 849046, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 826738, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 812638, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 804742, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 902277, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 902277, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 902276, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 902012, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 896818, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 877876, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 851673, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 830153, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 816497, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 808833, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Exception: optim_test_model_namespace::log_prob: total is 804778, but must be less than or equal to 802277.000000 (in '/tmp/Rtmp9znY5p/model-1192946a70e.stan', line 22, column 2 to column 58) 
Chain 1 Optimization terminated with error: 
Chain 1   Line search failed to achieve a sufficient decrease, no more progress can be made
Warning: Fitting finished unexpectedly! Use the $output() method for more information.
>

On the other hand, if I specify init = list(list(rates = rep(0, N))) it works but all parameters are set to 0, which is not useful at all.

2

Bounds aren’t calculated in the transformed parameters block, they are only enforced. To get Stan to enforce the bounds on the total you need to calculate each bound separately given the previous. You don’t need a uniform distribution prior on the total. With the code below it will automatically enforce that. Since there’s not additional data though I believe this optimization will always have the total be equal to the max.

Also, Stan will complain if the bound is included so we need to add a small amount to allow that.

This works where the bounds are calculated for each rate given the estimate of the previous rate.

functions {   
   real lb_ub_lp (real y, real lb, real ub) {
    target += log(ub - lb) + log_inv_logit(y) + log1m_inv_logit(y);
    
    return lb + (ub - lb) * inv_logit(y);
  }
}
data{
  // shape
  int<lower=0> N;
  
  // data
  real<lower=0> mult;
  vector[N] x;
  
  // constraints
  real<lower=0> max_value;
}
transformed data {
  vector[N] x_mult = x * mult;
}
parameters{
  vector[N] rates_raw;
}

transformed parameters{
  vector[N] rates;
  vector[N] x_rate_mult;
  
  rates[1] = lb_ub_lp(rates_raw[1], 0, max_value / x_mult[1]);
  x_rate_mult[1] = rates[1] * x_mult[1];
  {
    real max_ = max_value;
    for (i in 2:N) {
      max_ += -rates[i - 1] * x_mult[i - 1];
      rates[i] = lb_ub_lp(rates_raw[i], 0, max_ / x_mult[i]);
      x_rate_mult[i] = rates[i] * x_mult[i];
    }
  }
  real<lower=0, upper=max_value + 1> total = dot_product(rates, x_mult);
}
model{
  target += x_rate_mult;
  
}
1 Like
3

Hello @spinkney thanks! This works but the rates parameters are unbounded now, the sum is constrained but rates are not, they should be between 0 and 3. Seems like it is enforcing an equality constraint
and adjusting rates accordingly.

4

All of my rates are bounded. The difficulty is ensuring the upper bound of 3 just add a max(3, max_/x_mult).

1 Like
5

@spinkney amazing it works!! Had to use fmin though but it works perfectly without errors. Thanks!!

1 Like