Re: [eigen] Levenberg-marquardt and dogleg minimizer

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

On Sat, Feb 21, 2009 at 12:21 AM, Keir Mierle <mierle@xxxxxxxxx> wrote:
> [+libmv-devel, +eigen]
>
> On Fri, Feb 20, 2009 at 2:53 PM, Timothy Hunter <tjhunter@xxxxxxxxxxxx>
> wrote:
>>
>> Hello Keir,
>> thank you this post, it was a very interesting reading. I was trying to
>> find a good way to represent the gradient and the hessian of a function for
>> a convex solver, and I like your way of separating the jacobian from the
>> function object itself.
>
> I have a failed attempt which is not posted going down the other path :)

I have to say that I don't like the separation of f and the jacobian.
The current design works well with a numerical differentiation, but
when you compute the jacobian yourself it is very likely the case that
the intermediate results obtained during the computation of f can be
reused for the Jacobian.

Therefore I would suggest to compute both f and its Jacobian in a
single function call.

Moreover, I'd also suggest to pass the returned value using pointer
arguments (to avoid useless memory allocation and copies). So, for
instance, the called function could be:

void operator() (const XVectorType& x, FVectorType* f, JacobianType* j = 0)
{
  if (f) {
   // compute f
  }

  if (j) {
    // compute J
  }
}

Now the question is how to integrate a numerical differentiation or
auto-diff tool.

One possibility is to let the user call NumericalJacobian himself, in
the above function:

if(j) {
  m_numerical_jacobian(x, j);
}

where m_numerical_jacobian would be declared in the functor like this:

NumericalJacobian<MyFunctor>  m_numerical_jacobian(*this);

and with:

template<Func>
class NumericalJacobian {
const Func& m_f;
typename Func::FVectorType m_res1, m_res2;

void operator() (x,j) {
  // calls to:
  m_f(x +/- delta, m_res1);
  m_f(x +/- delta, m_res2);
  // update j with (m_res2 - m_res1)/delta....
}

};

Of course the above code is not legal C++, it is just to give the idea !

The main question is how to easily interact with an auto-diff tool.
Typically it would be nice to able to replace NumericalJacobian by a
AutoDiffJacobian class without any changes. I think that would
requires to make: Functor::operator() a template function such that
AutoDiffJacobian would be able to call f() with custom scalar types
(auto-diff tools basically used custom scalar types with overloaded
operators tracking the derivatives). I'm not sure that would work
though, and so I'd strongly suggest to try to build such an
AutoDiffJacobian class before establishing a final design (check the
AdolcForwardSupport module).

>>
>> I have a few suggestions:
>> - in the Jacobian class, instead of having a boolean as a second template
>> type, you may want a integer, so that different methods can be implemented
>> in the future
>
> This is a good idea, and as a bonus, will make calling code clearer, i.e.
>
> NumericJacobian<MyFunc, CENTRAL_DIFFERENCE> instead of
> NumericJacobian<MyFunc, false>
>
> I will make this change.

I agree too.

>>
>> - in a function, the typenames XMatrixType and FMatrixType are not so
>> obvious. Maybe InputType and ReturnType?
>
> I agree that FMatrixType and XMatrixType are not great, but I am not
> convinced InputType and ReturnType are better. Arguably the return type is
> return type of the minimization, which is the same type as the parameters (x
> in this case) and not the type of f(x). Note that inside the
> LevenbergMarquardt class, I typedef things to Parameters and FVec. Again,
> not optimal. I will think about this more.

perhaps ArgumentType and ValueType ?

Another thing is the ability to control the iterations. For instance
in the constrained gradient implementation I borrowed from GMM++ there
is an additional template argument taking an IterationController
object. Perhaps the same class could be used for both algorithm ? (it
would also replace your current Parameters class). Actually I did not
check at all if that makes sense or not, and I'm not saying the GMM++
way is better, it is just an idea.  What matters is to have a similar
mechanism for all algorithms.

cheers,
Gael.

>> Timothy Hunter
>>
>> On Thursday 19 February 2009 09:12:10 Keir Mierle wrote:
>> > The API for the LM implementation in libmv is more or less ready; since
>> > I
>> > plan to include it in Eigen at some point, I thought I'd send it along:
>> >
>> >
>> > http://code.google.com/p/libmv/source/browse/trunk/src/libmv/numeric/levenb
>> >erg_marquardt.h
>> >
>> > http://code.google.com/p/libmv/source/browse/trunk/src/libmv/numeric/levenb
>> >erg_marquardt_test.cc
>> >
>> > Mostly, I am looking for a review from an API perspective; also, I need
>> > to
>> > support sparse matrices at some point but don't have any experience with
>> > them; perhaps Gael could comment here as he is the sparse expert. In
>> > particular, I will need to support sparse Jacobians and solving sparse
>> > normal equations.
>> >
>> > I also have an implementation of Powell's dogleg method, which can be
>> > much
>> > faster than LM:
>> >
>> >
>> > http://code.google.com/p/libmv/source/browse/trunk/src/libmv/numeric/dogleg
>> >.h
>> >
>> > http://code.google.com/p/libmv/source/browse/trunk/src/libmv/numeric/dogleg
>> >_test.cc
>> >
>> > Keir
>>
>> --
>> Timothy Hunter
>> Student (Stanford University)
>>
>> T. 404 421 3075
>