# Linear Least Squares with Bounds

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.

# Robust nonlinear regression in scipy

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

# Large-Scale Least Squares

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.

# 2D Subspace Trust-Region Method

Trust-region type optimization algorithms solve the following quadratic minimization problem at each iteration: $\displaystyle \min m(p) = \frac{1}{2} p^T B p + g^T p, \text { s. t. } \lVert p \rVert \leq \Delta$

If such problem is too big to solve, the following popular approach can be used. Continue reading 2D Subspace Trust-Region Method

# Algorithm Benchmarks

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.

1. 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.
2. 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.

# Dogbox Algorithm

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 $\displaystyle \min_p m(p) = \frac{1}{2} p^T B p + g^T p \text{, s. t. } \tilde{l} \leq p \leq \tilde{u}$