Preconditioning

[cortner]

I was happy to find the CG has preconditioning implemented, but why stop there? Virtually all of the optimisation algorithms in this package can benefit from preconditioning. Personally, I'd be most interested in seeing it for L-BFGS. I can try to put together an implementation, but thought I should check first what others think about this.

I would propose to

• move the preconditioning functionality cg_precondfwd, cg_precondfwddot, cg_precondinvdot out of cg.jl into a new file precon.jl
• implement preconditioning for L-BFGS; this would require changing just the line 40, copy!(s, q) which effectively uses the identity matrix as preconditioner
• sketch out for other algorithms what would need to be done to give them preconditioning capabilities.

[pkofod]

Improvements are always welcomed, and it's great to see that you're willing to do a PR yourself.

It may be worth doing some benchmarks or examples showing how well it works, and, if it comes with any, the cost of adding preconditioning. I am aware that it can speed up convergence (in iterations) for ill-conditioned problems, but what cost does it add to optimizing well-conditioned problems? Should the identity preconditioning be default, with a user-choice of other preconditioners?

[cortner]

That is fair if course. But a proper implementation of preconditioning will have zero overhead in fact the only potential overhead is creating some extra copies which the current implementation does I think, and which L-BFGS does anyways in this line. As an initial implementation, I might keep this behaviour to be on the safe side.

[cortner]

Is there a standard set of benchmarks I can run?

[pkofod]

Okay, I may have slightly misunderstood your intent then. We have some simple unconstrained problems in https://github.com/JuliaOpt/Optim.jl/tree/master/src/problems, but it is on the drawing board to extend this part of the package over the coming months.

If you have something but no benchmarks feel free to submit a PR without them.

[cortner]

@pkofod (also cc @timholy - as I thought you implemented preconditioning in CG?)

Currently, the precondprep function is passed to the optimiser as parameter. To me it seems more natural to use dispatch here.

• ADVANTAGE: cleaner (in my view) design; further it would shield the user from having to supply the update functionality. E.g., one could then invoke mo = Optimiser(P=SomePrecon()) and the dispatch would automatically determine how to update it.
• DISADVANTAGE: more boiler-plate code. E.g., in fminbox.jl, one would need to define a new preconditioner type to create a dispatch that calls the correct precondprep implementation.

A related issue with the current implementation, as I see it, is that it is unclear what P is. Is it the primal metric or the dual metric? E.g., in Tim's implementation it is the dual metric, which I think is common in statistics, but in PDEs e.g. it is more common to work with the primal metric. Creating a layer of abstraction and fixing the terminology a bit could fix some of this.

Any views would be appreciated.

P.S.: For now I will make a pull request without making that change, since it will otherwise lie formant for another week or so, while the rest is actually done.

[cortner]

After this is all done: do you have any interest in Issue #188, or do you prefer I let it go?

[pkofod]

After this is all done: do you have any interest in Issue #188, or do you prefer I let it go?

I'm not sure I know what you mean. This is issue #188. Or do you mean the last bullet on your list?

[cortner]

sorry - I had meant to post this in the pull-request. I'd just like to know whether you have an interest in re-design of the preconditioning functionality, or whether you prefer to keep status quo.

It might be worth discussing this together with with the Linear Algebra packages.

[pkofod]

If you feel like there's more to discuss (for example adding preconditioning to other algorithms) I don't think you need to close this issue right now.

When you say Linear Algebra packages do you mean the linear algebra in the base library, or packages in general, that implement linear algebra routines? If it's the latter, I'm not sure such a coordination is feasible. Maybe I am misunderstanding you.

[cortner]

It is easy to add preconditioning to other algorithms. This is not so much of an issue.

What I am proposing here is to change the preconditioning infrastructure. I now notice a related discussion in IterativeSolvers.jl.

Specifically, I propose to define

Initialise: P = MyPrecon()
Update after state x has been updated: update!(P, x)
Apply preconditioner and take inner products: essentially as it is done now.

The main difference would be that precondprep!/ update! is no longer passed to the optimiser, but is chosen via dispatch. This seems to me is consistent with the redesign of the Optim.jl infrastructure anyhow.

[cortner]

I close this and start a new discussion.