optim-julep: define AbstractProblem and make use of it in optimizers

[timholy]

Today I encountered an optimization problem with special structure: the Hessian (which I believe I can calculate analytically) is, say, typically around 20k-by-20k (ranging from 2k-by-2k to 200k-by-200k), and essentially every entry is nonzero---so it is not sparse. However, it can be parametrized as a rank-4 modification of a block-diagonal matrix (with 3x3 blocks), and in this representation it is amenable to very efficient solution by the Woodbury matrix identity.

Assuming one can't coax Ipopt or similar packages to make use of such special structure, this reminds me yet again of the advantages of being able to perform optimization using purely julia code.

However, it looks like Optim may need some work to support this kind of thing. In particular, statements like this one will wreak havoc. I'm guessing the best approach would be to have problems specified as subtypes of an AbstractProblem type, and the user calls optimize(problem, x, options...). This could contain theTwiceDifferentiableFunction, but might also allow specialization of a method like H = allocate_hessian(problem). Does this seem reasonable?

I could just add this to #50 (which is the branch I use for all my own work anyway), but obviously I should eventually get around to putting that in shape for merging, so it would be best to get some consensus about any other APIs I introduce. This looks like a relatively independent refactoring from the very nice #149, and we might want to get that finished off before undertaking this.

[mlubin]

OOQP was designed to handle cases like this. Definitely a useful idea to abstract over the linear algebra to the extent possible.

[Evizero]

Personally, I would also prefer the route of optimizing types instead of functions (which is what you mean, I think?). This would allow for neat tricks when it comes to Machine Learning.

But as far as I was able to tell there was a strong preference towards optimizing Julia functions directly (#128 (comment)), so I abandoned that branch. Has this changed? I remember there was a lot of code to go through to abstract the f.f(...) / f.g!(...) calls into grad!(problem, ....) etc. but I still have that branch around locally

[timholy]

I like #128, thanks for pointing me to it. Very much in line with this proposal.

[pkofod]

Hm, I'm wondering if I can fit this into my GSoC schedule. If not, then maybe after. Could I convince @timholy and @Evizero to give brief examples of what you have in mind? For example, what would you wish you had had @timholy , and what kind of tricks @Evizero ?

[Evizero]

Hard to say now because AFAIK 0.5 brings new goodies to functions (i.e. the call replacement) that could make this proposal obsolete (?) Do not take my word for it, though, as I haven't spend any time investigating the change in detail.

Anyhow the idea was to replace the concrete call sites like https://github.com/JuliaOpt/Optim.jl/blob/master/src/bfgs.jl#L71 , with custom functions like grad!(gr, problem, x) (similar to https://github.com/Evizero/LearnBase.jl/blob/master/src/optim/optim.jl#L6-L26 , which was me testing that idea a while ago), where problem can be any custom data type. that would allow endusers to make use of dispatch to their liking.

In machine learning this would be especially nice since the optimization problems there usually depend on constant parameters (i.e. the data). Right now solving such problems require boilerplate and closures like this https://github.com/JuliaOpt/Optim.jl#optimizing-functions-that-depend-on-constant-parameters . Allowing custom types could hide this complexity, enable compile time dispatch, and allow for inlining.

So was my thinking at least.

[pkofod]

i think this is pretty much possible now. Just create your differentiable using the full constructor. Please reopen or open another issue if you have a specific problem (or you can share the original one)