Hi! The GSoC is coming to an end so this will be the last technical post. I am a bit tired of all this, so it will be short.
The last feature I was working on is robust loss functions support. The results are again available as IPython Notebook, look here https://gist.github.com/nmayorov/dac97f3ed9d638043191 (I’m struggling to get “&” work correctly in LaTeX blocks, so the formatting is a bit off at the moment). The plan is to provide this example as tutorial for scipy.
As a demonstration of large-scale capabilities of the new least-squares algorithms we decided to provide an example of solving a real industry problem called bundle adjustment. The most convenient form is IPython Notebook, later it will serve as an example/tutorial for scipy (perhaps hosted on some server). Here I just give a link to the static version http://nbviewer.ipython.org/gist/nmayorov/6098a514cc3277d72dd7
I finally made my code available as PR to scipy https://github.com/scipy/scipy/pull/5044 This PR contains all code, but was branched from the previous one and focuses on sparse Jacobian support. In this post I’ll explain the approach I chose to handle large and sparse Jacobian matrices.
Trust-region type optimization algorithms solve the following quadratic minimization problem at each iteration:
If such problem is too big to solve, the following popular approach can be used. Continue reading 2D Subspace Trust-Region Method
This post was updated to improve its clarity and to incorporate new information about “leastsqbound”.
Before I present the results I want to make a few notes.
- Initially I wanted to find a very accurate reference optimal value for each problem and measure the accuracy of an optimization process by comparison with it. I abandoned this idea for several reasons. a) In local optimization there isn’t a single correct minimum, all local minima are equivalently good. So ideally we should find all local minima which can be hard and comparison logic with several minima becomes awkward. b) Sources with problem descriptions often provide inaccurate (or plain incorrect) reference values or provide them with single precision, or doesn’t provide them at all. Finding optimal values with MATLAB (for example) is cumbersome and still we can’t 100% assure the required accuracy.
- It is desirable to compare algorithms with identical termination conditions. But this requirement is never satisfied in practice as we work with already implemented algorithms. Also there is no one correct way to specify termination condition. So the termination conditions for all algorithms will be somewhat different, but nothing we can do about it.
The idea for this simple algorithm is taken from this paper. We use a rectangular trust region, so intersection of a trust region and a rectangular feasible region is again some rectangle. Thus at each iteration we need to solve the following constrained quadratic problem
The most relevant description of this algorithm can be found in the paper “A subspace, interior and conjugate gradient method for large-scale bound-constrained minimization problems” by Coleman and Li, some insights on its implementation can be found in MATLAB documentation here and here. The difficulty was that the algorithm incorporates several ideas, but it was not very clear how to combine them all together in the actual code. I will describe each idea separately and then outline the algorithm in general. I will consider the algorithm applied to a problem with known Hessian, in least squares we replace it by . I won’t give any explanation or motivation for some things, if you are really interesting try digging into the original papers.
Hi! The last two weeks in GSoC I spent time developing algorithm drafts, testing and tuning them. It is an addicting process and I spent perhaps already too much time on it. The results are pretty good, but now I seriously need to stop it and start writing production quality code and so on.
The first thing we decided I should do is to provide a plain English description of the methods I studied and implemented. To make it more understandable I decided to make this post as an introduction to nonlinear least-squares optimization. So here I start.
As nothing exciting is done yet I will just briefly review the current situation and perspectives.