Newton solver in Stan

I am using MacPorts and never had trouble installing cmake with it. Maybe use that?

Try installing Sphinx 1.3.1 with Anaconda?

I took an odd path: I “lost” homebrew, as in brew was no longer a valid command on my terminal. So I re-installed homebrew and ran brew install cmake. This time, I got the following message: cmake 3.13.4 is already installed and up-to-date. Ok, good enough. Not sure what happened.

Updates

I have some working code in C that uses KINSOL to solve a simple system, see laplace_code.txt (4.2 KB). Note the extension should be .c.

The problem is the code is quite different than what we would do in C++ with Stan. Stan’s CVODES integrators pave the way on how to do this, but we’ll still need to discuss the specs. Also are there functions in cvodes_ode_data that we can use for the KINSOL solver?

I have tests in Stan, with which we can evaluate KINSOL’s performance (and furthermore, try out its different features). So now it’s a matter of linking to the KINSOL library, and converting Stan input (written in C++, with Eigen matrices, etc.) into C code KINSOL can handle.

I think you can wrap the N_Vector things from sundials around Eigen matrices. At least I have setup an N_Vector to use the data stored in a std::vector as I recall in the cvode_ode_data object, I think.

I am curious as to how good KINSOL performs… is the examples directory from KINSOL helpful?

is the examples directory from KINSOL helpful?

Yes, it got me started fairly quickly (a few hours) and the accompanying documentation is good. It also walks you through a bunch of features you can try out with KINSOL.

But I haven’t worked through all the examples.

Here’s a compiling error that I could use help with. I’m re-writing the above C-code inside a Stan unit test. The goal is to expose kinsol to Stan. I simply added the kinsol directories under the Sundial library. When compiling my unit test, I get the following error:

Undefined symbols for architecture x86_64:
  "_N_VNew_Serial", referenced from:
      laplace_kinsol_Test::TestBody() in i_test.o
ld: symbol(s) not found for architecture x86_64
clang: error: linker command failed with exit code 1 (use -v to see invocation)

This error extends to other “symbols”, such as KINCreate(), KINInit(), etc. Now, some of these should work, since they are used in other test files. For example, N_VNew_Serial() appears in idas_system_test.cpp. The header files between this test and mine are comparable. I actually got the compiler to recognize some of the symbols (but not all of them) when renaming (!) my test idas_test.cpp. Conversely, I got idas_system_test.cpp to fail by changing the file’s name. Either manipulation shouldn’t have zero effect – unless I’m missing something.

In make/tests there’s

IDAS_TESTS := $(subst .cpp,$(EXE),$(shell find test -name *idas*_test.cpp) $(shell find test -name *_dae*_test.cpp))
$(IDAS_TESTS) : $(LIBSUNDIALS)

that indicates $(LIBSUNDIALS) dependency. You’ll need something similar to include sundials for your KINSOL tests.

1 Like

@yizhang thanks for the pointer. I ended up editing make/libraries and make/tests to expose KINSOL and link the unit test.

@wds15 I closely followed what you did with cvodes_data and created a prototype kinsol_data. Thank you for the properly written and well documented code. kinsol_daa is a bit awkward in certain places, because of the need to handle various vector types (std::vector, Eigen::VectorXd, N_Vector), so this will need cleaning up.

Here’s a test with the prototype. Not perfect, but should be enough to run performance tests.

For the sensitivities, I’ll simply use the vari class we already have in place for algebra_solver().

Glad it helped. @betanalpha was the first to start this - I did refine that code quite a bit. Anyway… finally someone saying good things about this code. Others were very confused by it (it’s not the easiest for sure).

Results??? Results?? I am curious how KINSOL does - looking forward to numbers.

The issue with cvode_data is that to understand it, you probably need a good grasp of how sundials works and how the sensitivities get computed. Neither is trivial. But, as far documenting a hard problem goes, I think this was done quite well.

Regarding performance: right now, I’m getting discrepancies in the output when benchmarking against Powell and a custom Newton solver. The latter two are in agreement. Things like 0.084199049 vs 0.084199123, and it gets worse in high dimensions. So I need to fix this before looking at performance.

Note related to reported performance but you may have a memory leak as the destructor


doesn’t clean up N_Vector nv_x_.

EDIT: same problem happens to kinsol_memory.

Are you sure? Quickly looking over it, I think that the memory is not managed by the NVector, but rather by the x_ object. The N_VMake_serial only wraps around it, but the data is not owned by it, I think.

Do you mean in line 199 nv_x_(N_VMake_Serial(N_, &to_array_1d(x_)[0])), x_ would manage the data of std::vector allocated from to_array_1d? It doesn’t make sense to me.

Ah… yeah. The issue is the temporary - that’s right, but it does not help to clean up the N_Vector.

The better thing is probably to wrap N_Vector around the Eigen x_ input argument directly. Then things should be fine.

@wds15 yep, that’s better.

@charlesm93 Since KINSetMaxSetupCalls is set to 1, it’s exact Newton where Jacobian gets updated every step. How’s that going to affect the performance of a large system? Will modified newton achieve the same, say, for Laplace approximation, with better efficiency?

Something simple like nv_x_(N_VMake_Serial(N_, &x_(0))) does not compile, so I’ll need to use Eigen::Map to wrap x_ directly. Added a call to free kinsol_memory . Also moved code to header files, and outside the unit tests, and expanded kinsol_test (note: the test inla system is not written optimally, but it works as a test).

Since KINSetMaxSetupCalls is set to 1, it’s exact Newton where Jacobian gets updated every step. How’s that going to affect the performance of a large system? Will modified newton achieve the same, say, for Laplace approximation, with better efficiency?

There are several features to try out, including inexact Newton and linesearch methods. This is the next step in the performance tests. Eventually we’ll need to figure out which tuning parameters we expose, what default values we set, etc.

I’ve dug a bit deeper into why the kinsol solver disagrees with the Newton solver, and it seems to be a matter of adjusting the tuning parameters, so that both solvers have the same “precision”. In low dimension, that’s not hard to tweak.

When solving a 100-d system however, I get a flag -7 error: “Five consecutive steps have been taken that satisfy a scaled step length test.” or in more details “Five consecutive steps have been taken that satisfy the inequality ∥D_u p∥_{L2} > 0.99 mxnewtstep, where p denotes the current step and mxnewtstep is a scalar upper bound on the scaled step length. Such a failure may mean that ∥D_F F (u)∥_{L2} asymptotes from above to a positive value, or the real scalar mxnewtstep is too small”. Haven’t figured out what this means, but I’m digging through the manual.

A bit more on mxnewstep, the maximum allowable scaled length of the Newton step. According to the manual, the default value is 1000 ||u_0||_D, where u_0 is the initial guess. The problem is I pass a vector of 0 as my initial guess (a reasonable guess), hence mxnewstep = 0. If I pass in 1000||u||_D, where u is the solution I obtained from another solver, the kinsol solver produces the right value. So the default needs to be adjusted when the initial guess has norm 0, though I’m not sure what an optimal default might be.

Some preliminary performance results. Kinsol is challenging because it has a lot of tuning parameters. Not tuning them well leads to error or suboptimal performance. Note kinsol always gives an error message when it fails. Coming up with the right implementation will require some work.

For these tests, the linesearch method is kinsol’s defaut, KIN_LINESEARCH. The quasi Newton method is implemented by setting mbset = 5, which yields stable results. The problem at hand is a Laplace approximation – not written in the most optimal manner – with dimensions varying from 10 to 500. I got an error message (flag -7) with linesearch and quasi methods at dimension 500.

performance01
performance02 performance03

Note: these are evaluation times. The sensitivities are computed separately and in the same way for all solvers, using the same vari class.

2 Likes

@yizhang, @wds15, and I discussed an API for the kinsol solver. A lot of features seem overkill for users. We agreed on a basic implementation: exact Newton (i.e. Jacobian computed at each step) and the usual tuning parameters (relative and absolute tolerance, and max num steps).

Staring at the graphs I posted above, I now think offering linesearch and quasi-newton would be nice features too. If others corroborate I’ll add these in, else I’ll roll with the basic implementation.