# Exposing quadratic optimizer

**URL:** <https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684>\
**Category:** Developers\
**Created:** [June 29, 2018, 1:02am UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684 "2018-06-29T01:02:41Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![charlesm93](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charlesm93/32/4634_2.png) [@charlesm93](https://discourse.mc-stan.org/u/charlesm93)\
**Post date:** [June 29, 2018, 1:02am UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684/1 "2018-06-29T01:02:41Z")

</div>

I’m trying to expose a higher-order function (the quadratic optimizer with linear constraints, see [https://github.com/stan-dev/stan/issues/2556](https://github.com/stan-dev/stan/issues/2556)). I followed closely what was done for `algebra_solver` and then ran an already existing unit test, using:

`./runTests.py src/test/unit/lang/parser/matrix_functions_test.cpp`

The goal was to see if the code would compile. It doesn’t. The error message is in the attached file was return. [error\_message.txt](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/6/69423d71939154e30d027043c6ac31e31768b5f5.txt) (103.5 KB)

Ok. Both errors point to the same line of code, which I wrote: `term_grammar_def.hpp:287:9`: As far as I can tell, this line is fine and the error lies somewhere else. I’m not quiet sure where and how to look. After going through the files I made a few corrections, but none that fixed the above error. Note there is a total of 21 edited files.

Does anyone more familiar with the language parser have an idea about what may be causing the error? Thank you!

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [July 20, 2018, 4:34pm UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684/2 "2018-07-20T16:34:51Z")

</div>

Sorry I missed this, @charlesm93. If you haven’t sorted it out yet, could you push a branch and email me? I can check it out.

---

<div class="post-metadata">

**Author:** ![charlesm93](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charlesm93/32/4634_2.png) [@charlesm93](https://discourse.mc-stan.org/u/charlesm93)\
**Post date:** [July 20, 2018, 11:41pm UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684/3 "2018-07-20T23:41:28Z")

</div>

Hey Bob,  
branch is pushed and linked in my original post. I haven’t worked it out yet. I might give it a stab on the airplane.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [July 21, 2018, 10:41am UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684/4 "2018-07-21T10:41:46Z")

</div>

> [@charlesm93](#):
>
> branch is pushed and linked in my original post.

You linked the issue. The branch is [feature/issue-2556-quadratic\_optimizer](https://github.com/stan-dev/stan/tree/feature/issue-2556-quadratic_optimizer).

I’d recommend simplifying things. Remove the semantic actions to make sure the FUSION declarations and grammar rule types are all correct. Then compile the function and make sure that works. Then try to tie it into the grammar. It’s painful, but it’s the only way we were ever able to debug these Boost Spirit Qi grammars (in @syclik’s words, you want to take “baby steps”).

So I checked out your branch and removed the quadratic optimizer semantic action declaration, definition, and call from the grammar and everything compiles other than `generate_expression` because the case isn’t defined for the quadratic optimzer. Fix that (with a stub), then start with a stub for the semantic action to get it to compile. Then start filling it in piece by piece.

---

<div class="post-metadata">

**Author:** ![charlesm93](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charlesm93/32/4634_2.png) [@charlesm93](https://discourse.mc-stan.org/u/charlesm93)\
**Post date:** [July 31, 2018, 12:16am UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684/5 "2018-07-31T00:16:52Z")

</div>

Following Bob’s advice to break things down, I was able to fix a myriad of small errors which ultimately caused the bug. The changes were committed and the code now compiles.

It however seems I’m still missing a piece to expose the function. See the following unit test `src/test/unit/lang/parser/quadratic_optimizer_test.cpp`, according to which `quadratic_optimizer` is not a known function. Which file would the parser use to check if a higher-order function exists? Do I perhaps need to recompile certain files?

Tagging @mitzimorris.

---

<div class="post-metadata">

**Author:** ![mitzimorris](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mitzimorris/32/29_2.png) [@mitzimorris](https://discourse.mc-stan.org/u/mitzimorris)\
**Post date:** [July 31, 2018, 8:07am UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684/7 "2018-07-31T08:07:16Z")

</div>

checked out branch feature/issue-2556\_etc  
but looks like you haven’t added (or committed or pushed) file  
`src/test/unit/lang/parser/quadratic_optimizer_test.cpp`

cheers,  
Mitzi

---

<div class="post-metadata">

**Author:** ![mitzimorris](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mitzimorris/32/29_2.png) [@mitzimorris](https://discourse.mc-stan.org/u/mitzimorris)\
**Post date:** [August 1, 2018, 5:54pm UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684/8 "2018-08-01T17:54:18Z")

</div>

@charlesm93 - there were a few remaining small errors -

1. inconsistant name for semantic actions function `validate_quadratic_optimizer_control_f`

2. test-model was missing arg #8

also tweaked the grammar rule so that it will backtrack if it hits functions named something like `quadratic_optimizer_with_fries(` etc.

committed and pushed my changes - hope this helps.  
cheers,  
Mitzi

---

<div class="post-metadata">

**Author:** ![mitzimorris](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mitzimorris/32/29_2.png) [@mitzimorris](https://discourse.mc-stan.org/u/mitzimorris)\
**Post date:** [August 2, 2018, 11:52am UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684/9 "2018-08-02T11:52:36Z")

</div>

here is what I did to uncover these problems:

1. started looking at the unit test - cross-referenced the line in the test-model file with the quadratic optimizer code and noticed that there was a missing argument.

2. uncommented the grammar rule in `src/stan/lang/grammars/term_grammar_def.hpp`

at this point, the code wasn’t compiling. captured the c++ boost/spirit/qi compiler error novel to an output file. in the error novel, looked for the very first line of spew that referenced files in `src/stan/lang/grammars` which showed that there was a problem with the term grammar rule. careful cross-referencing of the semantic action and its corresponding declarations and definitions in the `semantic_actions.hpp`, `semantic_actions_def.cpp` files showed that the problem was with the function name.

cheers,  
Mitzi

---

<div class="post-metadata">

**Author:** ![mitzimorris](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mitzimorris/32/29_2.png) [@mitzimorris](https://discourse.mc-stan.org/u/mitzimorris)\
**Post date:** [August 3, 2018, 2:06pm UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684/10 "2018-08-03T14:06:47Z")

</div>

added tests for all 11 or so ways that one could mess up on the arguments to the quadratic\_optimizer function call.

@Bob_Carpenter are there other tests that should be done w/r/t adding this to the language?

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [August 11, 2018, 3:15pm UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684/11 "2018-08-11T15:15:00Z")

</div>

Not that I’m aware of.

---

<div class="post-metadata">

**Author:** ![charlesm93](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charlesm93/32/4634_2.png) [@charlesm93](https://discourse.mc-stan.org/u/charlesm93)\
**Post date:** [August 19, 2018, 3:56pm UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684/12 "2018-08-19T15:56:52Z")

</div>

@mitzimorris I finally got a chance to check out your commits (have been a little mia this month). Thank you so much for resolving the issue.

---

<div class="post-metadata">

**Author:** ![andre.pfeuffer](https://avatars.discourse-cdn.com/v4/letter/a/e480ec/32.png) [@andre.pfeuffer](https://discourse.mc-stan.org/u/andre.pfeuffer)\
**Post date:** [June 16, 2019, 4:00am UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684/13 "2019-06-16T04:00:43Z")

</div>

Is there a chance to get the solve in the next Stan release?  
I would be beneficial for function approximation as asked here:  
[https://discourse.mc-stan.org/t/curve-fitting-including-convexity-constraint/9077/15](https://discourse.mc-stan.org/t/curve-fitting-including-convexity-constraint/9077/15)

---

<div class="post-metadata">

**Author:** ![charlesm93](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charlesm93/32/4634_2.png) [@charlesm93](https://discourse.mc-stan.org/u/charlesm93)\
**Post date:** [June 17, 2019, 6:33am UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684/14 "2019-06-17T06:33:39Z")

</div>

There are technical issues for the general implementation of a quadratic optimizer with linear constraints. Yes, we can solve the problem, and yes we can compute derivatives where they exist. But the resulting posterior is non-smooth at the boundaries, which is a problem for HMC. See this [discussion here](https://discourse.mc-stan.org/t/non-smooth-posterior-and-kkt-problem/6281).

---

<div class="post-metadata">

**Author:** ![andre.pfeuffer](https://avatars.discourse-cdn.com/v4/letter/a/e480ec/32.png) [@andre.pfeuffer](https://discourse.mc-stan.org/u/andre.pfeuffer)\
**Post date:** [June 17, 2019, 7:42am UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684/15 "2019-06-17T07:42:35Z")

</div>

But then I would expect this is also true for the current algebra\_solver. Depending on the  
tolerance \epsilon \> 0 this will create a non-smooth posterior. Is that true?

---

<div class="post-metadata">

**Author:** ![charlesm93](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charlesm93/32/4634_2.png) [@charlesm93](https://discourse.mc-stan.org/u/charlesm93)\
**Post date:** [June 17, 2019, 9:49am UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684/16 "2019-06-17T09:49:42Z")

</div>

No. \epsilon is not a model parameter, you can only pass it as a `double`.

Maybe I misunderstood your query?

---

<div class="post-metadata">

**Author:** ![betanalpha](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/betanalpha/32/4_2.png) [@betanalpha](https://discourse.mc-stan.org/u/betanalpha)\
**Post date:** [June 17, 2019, 2:32pm UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684/17 "2019-06-17T14:32:31Z")

</div>

The algebraic solver requires full rank algebraic systems, otherwise the Jacobian is ill-posed and the function errs out. The quadratic solver is full rank in the interior but becomes rank deficient on the boundaries. In other words there are no inconsistencies.

---

<div class="post-metadata">

**Author:** ![andre.pfeuffer](https://avatars.discourse-cdn.com/v4/letter/a/e480ec/32.png) [@andre.pfeuffer](https://discourse.mc-stan.org/u/andre.pfeuffer)\
**Post date:** [June 18, 2019, 3:22am UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684/18 "2019-06-18T03:22:09Z")

</div>

Thank you. I was thinking about the solver as a numerical approximation. The solver may stop in any interval [-\epsilon, +\epsilon]. What happen to the derivatives? They may be far off.  
This all would be no problem, if not my (or others) Stan programs gets stuck from time to time using  
the algebra solver.

---

<div class="post-metadata">

**Author:** ![charlesm93](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charlesm93/32/4634_2.png) [@charlesm93](https://discourse.mc-stan.org/u/charlesm93)\
**Post date:** [June 18, 2019, 6:40am UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684/19 "2019-06-18T06:40:23Z")

</div>

It’s true that the solution and furthermore all the derivatives we compute with autodiff are numerical approximations. What matters is that the error is not too large.

Ideally, you would set \epsilon to be small enough to get “sufficient” precision (i.e. your inference doesn’t get better if you make \epsilon smaller). The error in the solution propagates to the derivatives. So it’s not that the posterior is not smooth, rather that we’re doing a poor job evaluating the log density and the gradient.

> my (or others) Stan programs gets stuck from time to time using the algebra solver.

Could you describe the exact symptoms? I’ve observed shortcomings with the solver and there may be tricks to get improved behavior. I’m also working on a better solver.

---

<div class="post-metadata">

**Author:** ![andre.pfeuffer](https://avatars.discourse-cdn.com/v4/letter/a/e480ec/32.png) [@andre.pfeuffer](https://discourse.mc-stan.org/u/andre.pfeuffer)\
**Post date:** [June 18, 2019, 10:15am UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684/20 "2019-06-18T10:15:12Z")

</div>

The Stan program with solver runs forever in the sampling period. I’m calculating the inverse of some log\_sum\_exp expected value problem in a planning system. It’s Markov chain systems (M/M systems).  
The equations should be well behaved.

Currently I’m solving.

```
  vector algebra_system(vector y, // unknowns
                        vector theta, // parameters
                        real[] dat,
                        int[] dat_int) {
    vector[2] EV = E_value(y[1], y[2], theta[26:28], theta[1:25]);
    vector[2] f_x;
    f_x[1] = dat[1] - EV[1];
    f_x[2] = dat[2] - EV[2];
    return f_x;
  }

```

(I’d rather not like to share the function E\_value, not because I’m not willing, but to avoid intellectual problems.)  
dat is my target, but I could also minimize:

(EV[1] - dat[1])^2 + (EV[2] - dat[2])^2

I naively changed Stan’s algebra\_solver.hpp to try use Eigen’s LevenbergMarquardt Algorithm:

line 156  
Eigen::LevenbergMarquardt solver(fx);  
line 167  
solver.minimize(theta\_dbl);

Since it also “solves” minimizes to f(x) = 0 by minimizing the squared error at 0.

See references here: [https://stackoverflow.com/questions/18509228/how-to-use-the-eigen-unsupported-levenberg-marquardt-implementation](https://stackoverflow.com/questions/18509228/how-to-use-the-eigen-unsupported-levenberg-marquardt-implementation)

But it gives me some resize() exception. Although I admit I didn’t understand what’s going on in  
the C++ file in all details.

---

<div class="post-metadata">

**Author:** ![charlesm93](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charlesm93/32/4634_2.png) [@charlesm93](https://discourse.mc-stan.org/u/charlesm93)\
**Post date:** [June 18, 2019, 11:10am UTC](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684/21 "2019-06-18T11:10:36Z")

</div>

You could try writing a vanilla Newton solver in Stan’s `functions` block and use this instead of the algebraic solver. If the problem is quasi-convex and doesn’t require too many iterations to solve, you might get some speed up. I have a branch with the Newton solver in `math`, but I haven’t exposed it to the language yet.

You should also make sure the calculations in `E_value` are efficient. Bottlenecks can arise with matrix calculations. Some general guidelines: avoid direct matrix inversion, use functions that take advantage of symmetry, etc. Of course, without seeing `E_value` it’s hard to prescribe anything specific.

[Next page](https://discourse.mc-stan.org/t/exposing-quadratic-optimizer/4684.md?page=2)