[Numpy-discussion] Solving a NLLSQ Problem by Pieces?

Anne Archibald aarchiba@physics.mcgill...
Sat Jun 26 12:32:44 CDT 2010

The basic problem with nonlinear least squares fitting, as with other
nonlinear minimization problems, is that the standard algorithms find
only a local minimum. It's easy to miss the global minimum and instead
settle on a local minimum that is in fact a horrible fit.

To deal with this, there are global optimization techniques -
simulated annealing, genetic algorithms, and what have you (including
the simplest, explore the whole space with a "dense enough" grid then
fine-tune the best one with a local optimizer) - but they're
computationally very expensive and not always reliable. So when it's
possible to use domain-specific knowledge to make sure what you're
settling on is the real minimum, this is a better solution.

In this specific case, as in many optimization problems, the problem
can be viewed as a series of nested approximations. The crudest
approximation might even be linear in some cases; but the point is,
you solve the first approximation first because it has fewer
parameters and a solution space you understand better (e.g. maybe you
can be sure it only has one local minimum). Then because your
approximations are by assumption nested, adding more parameters
complicates the space you're solving over, but you can be reasonably
confident that you're "close" to the right solution already. (If your
approximations are "orthogonal" in the right sense, you can even fix
the parameters from previous stages and optimize only the new ones; be
careful with this, though.)

This approach is a very good way to incorporate domain-specific
knowledge into your code, but you need to parameterize your problem as
a series of nested approximations, and if it turns out the corrections
are not small you can still get yourself into trouble. (Or, for that
matter, if the initial solution space is complex enough you can still
get yourself in trouble. Or if you're not careful your solver can take
your sensible initial guess at some stage and zip off into never-never
land instead of exploring "nearby".)

If you're interested in how other people solve this particular
problem, you could take a look at the open-source panorama stitcher
"hugin", which fits for a potentially very large number of parameters,
including a camera model.

To bring this back nearer to on-topic, you will naturally not find
domain-specific knowledge built into scipy or numpy, but you will find
various local and global optimizers, some of which are specialized for
the case of least-squares. So if you wanted to implement this sort of
thing with scipy, you could write the domain-specific code yourself
and simply call into one of scipy's optimizers. You could also look at
OpenOpt, a scikit containing a number of global optimizers.

P.S. This question would be better suited to scipy-user or astropy
than numpy-discussion. -A

On 26 June 2010 13:12, Wayne Watson <sierra_mtnview@sbcglobal.net> wrote:
> The context here is astronomy and optics. The real point is in the last
> paragraph.
> I'm looking at a paper that deals with 5 NL (nonlinear) equations and 8
> unknown parameters.
> A. a=a0+arctan((y-y0)/(x-x0)
> B. z=V*r+S*e**(D*r)
>    r=sqrt((x-x0)**2+(y-y0)**2)
> and
> C.  cos(z)=cos(u)*cos(z)-sin(u)*sin(ep)*cos(b)
>     sin(a-E) = sin(b)*sin(u)/sin(z)
> He's trying to estimate parameters of a fisheye lens which has taken
> star images on a photo plate. For example, x0,y0 is the offset of the
> center of projection from the zenith (camera not pointing straight up in
> the sky.) Eq. 2 expresses some nonlinearity in the lens.
> a0, xo, y0, V, S, D, ep, and E are the parameters. It looks like he uses
> gradient descent (NLLSQ is nonlinear least squares in Subject.), and
> takes each equation in turn using the parameter values from the
> preceding one in the next, B. He provides reasonable initial estimates.
> A final step uses all eight parameters. He re-examines ep and E, and
> assigns new estimates. For all (star positions) on the photo plate, he
> minimizes SUM (Fi**2*Gi) using values from the step for A and B, except
> for x0,y0. He then does some more dithering, which I'll skip.
> What I've presented is probably a bit difficult to understand without a
> good optics understanding, but my question is something like is this
> commonly done to solve a system of NLLSQ? It looks a bit wild. I guess
> if one knows his subject well, then bringing some "extra" knowledge to
> the process helps. As I understand it, he solves parameters in A, then
> uses them in B, and so on. I guess that's a reasonable way to do it.
> --
>            Wayne Watson (Watson Adventures, Prop., Nevada City, CA)
>              (121.015 Deg. W, 39.262 Deg. N) GMT-8 hr std. time)
>               Obz Site:  39° 15' 7" N, 121° 2' 32" W, 2700 feet
>               Air travel safety Plus Three/Minus 8 rule. Eighty %
>               of crashes take place in the first 3 minutes and
>               last 8 minutes. Survival chances are good in you're
>               paying attention. No hard drinks prior to those
>               periods. No sleeping pills either. Sit within 5 rows
>               of an exit door. --The Survivors Club, Ben Sherwood
>                     Web Page:<www.speckledwithstars.net/>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

More information about the NumPy-Discussion mailing list