Cmdstanr: The optimize() function doesn't return output when it fails to converge

[syclik]

Describe the bug

When running optimize using CmdStan, it can return a value even if a global optimum isn't reached.

1. CmdStan does write to file, even if algorithm diagnostics state that it may not have found the right optimum.
2. RStan currently returns the value.
3. CmdStanR does not. It does report the right information.

To Reproduce

I used 8 schools and ran optimization on it.

From R (after setting up model):

> fit = model$optimize(data = schools_dat)
method = optimize 
  optimize 
    algorithm = lbfgs (Default) 
      lbfgs 
        init_alpha = 0.001 (Default) 
        tol_obj = 9.9999999999999998e-13 (Default) 
        tol_rel_obj = 10000 (Default) 
        tol_grad = 1e-08 (Default) 
        tol_rel_grad = 10000000 (Default) 
        tol_param = 1e-08 (Default) 
        history_size = 5 (Default) 
    iter = 2000 (Default) 
    save_iterations = 0 (Default) 
id = 1 
data 
  file = /tmp/RtmplVvoS1/standata-14844a9babf.json 
init = 2 (Default) 
random 
  seed = 1389664808 
output 
  file = /tmp/RtmplVvoS1/eight_schools_4-202010090300-1-a51a43.csv 
  diagnostic_file =  (Default) 
  refresh = 100 (Default) 

Initial log joint probability = -47.8461 
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
      99       228.766        0.6228   8.01767e+13      0.3575      0.3575      259    
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
     118       252.377   0.000284799   9.72824e+14       1e-12       0.001      387  LS failed, Hessian reset  
Optimization terminated with error:  
  Line search failed to achieve a sufficient decrease, no more progress can be made 
Warning: Fitting finished unexpectedly!

> fit
Error: Fitting failed. Unable to print.
> 

Expected behavior

Even though it isn't able to find an optimum, it should return what's in output.csv.

Compare this to rstan:

> fit = rstan::stan("eight_schools.stan", data = schools_dat)
...
> optimum = rstan::optimizing(fit@stanmodel, data = schools_dat)
Warning message:
In .local(object, ...) : non-zero return code in optimizing
> optimum
$par
    theta[1]     theta[2]     theta[3]     theta[4]     theta[5]     theta[6]     theta[7]     theta[8]           mu 
1.288251e+01 1.288251e+01 1.288251e+01 1.288251e+01 1.288251e+01 1.288251e+01 1.288251e+01 1.288251e+01 1.288251e+01 
         tau 
3.353735e-15 

$value
[1] 240.5844

$return_code
[1] 70

$theta_tilde
     theta[1] theta[2] theta[3] theta[4] theta[5] theta[6] theta[7] theta[8]       mu          tau
[1,] 12.88251 12.88251 12.88251 12.88251 12.88251 12.88251 12.88251 12.88251 12.88251 3.353735e-15

> 

Operating system
RStudio Cloud

CmdStanR version number

> packageVersion("cmdstanr")
[1] ‘0.1.3’

[mitzimorris]

this might be caused by fact that optimize method has non-zero return when it fails to converge - see discussion here: https://github.com/stan-dev/cmdstanpy/pull/75

tldr; CmdStan non-zero return codes reflect posix return codes which aren't designed to be user-facing. also, use of return codes in both CmdStan and underlying Stan services methods are inconsistent and not well-documented.

[rok-cesnovar]

@syclik tried looking into this but cant find a model to not converge.

Mind posting the exact eight schools model and data that does not converge? Would help with fixing this a lot. Thanks!

[jgabry]

@rok-cesnovar Any version of the 8 schools model should fail with optimization, so I think you can just use the version we have via the cmdstanr_example() function:

fit <- cmdstanr_example("schools", method = "optimize")

[syclik]

@jgabry, thanks!

@rok-cesnovar, yes, the schools example will fail. MCMC will converge, but optimization will fail. Or... rather, it'll try to get tau to be 0.

[jgabry]

This should be fixed on master now that Rok's #329 is merged. Thanks again @syclik for reporting this!

> fit <- cmdstanr_example("schools", method = "optimize", quiet = FALSE)

Model executable is up to date!
Initial log joint probability = -52.4781 
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
      99       121.977      0.191486   2.56133e+09      0.3104      0.3104      167    
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
     189       255.122     0.0029351   8.58686e+16       1e-12       0.001      384  LS failed, Hessian reset  
Optimization terminated with error: 
  Line search failed to achieve a sufficient decrease, no more progress can be made
Finished in  0.1 seconds.

> fit$print()

 variable estimate
 lp__       255.12
 mu           0.48
 tau          0.00
 theta[1]     0.48
 theta[2]     0.48
 theta[3]     0.48
 theta[4]     0.48
 theta[5]     0.48
 theta[6]     0.48
 theta[7]     0.48

 # showing 10 of 11 rows (change via 'max_rows' argument)