Re: [eigen] Levenberg-marquardt and dogleg minimizer

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

On Mon, Feb 23, 2009 at 8:32 PM,  <tjhunter@xxxxxxxxxxxx> wrote:
> We are starting to have a lot of ideas to put together :)
> I agree with Gael, I remember some cases where you can compute the value and
> the jacobian (or the gradient) at the same time. That being said, it will
> complicate the automation of differentiation. The main design issue I see
> here is: do we assume that the writer of the functor object will call
> himself provided differentiation tools? From a user perspective, it would be
> nice if we could require the most basic information only to be present, such
> as:
>
> class MyFunctor
> {
> // some enum variables
>
> // some typedefs
>
> void operator()(const SomeInputType &, SomeOutputType *)
> {
> //main code here
> }
> };
>
> and at the same time be able to find out if the user provides a computation
> of the gradient/jacobian/hessian like:
>
>
> class MyMoreComplicatedFunctor
> {
> // some enum variables
>
> // some typedefs
>
> void operator()(const SomeInputType &, SomeOutputType *=NULL,
> SomeGradientType *=NULL, SomeHessianType *=NULL)
> {
> //main code here
> }
> };
>
>
> My knowledge of c++ is still limited so I am not sure to see how to do it. I
> guess there is no way to detect if the operator() is overloaded, so maybe
> some additional enum flags or inheriting from a base class? Furthermore, I
> am not sure to see how adolc would fit here.
> Given your experience, do you think such minimal requirements are possible?

this would be easy to do via optional base classes for the functors,
to give an idea:

template<typename Derived> class NumericalJacobianBase
{
  void operator() (x, f, j)
  {
      static_cast<Derived&>(*this).operator()(x,f);
      NumeriaclJacobian<Derived>()(static_cast<Derived&>(*this),x,j);
  }
};

class myfunc : public NumericalJacobianBase<myfunc> {
  void operator() (x, f)
  {
    ....
  }
};

as soon as I have little time I'll try to implement such an autodiff
Jacobian class using adolc to see what are the exact requirements
(according to my experience with adolc this should be easy, but we
never know).

> Best regards
>
> Timothy Hunter
>
>
> On Saturday 21 February 2009 02:40:33 Gael Guennebaud wrote:
>> 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/l
>> >> >evenb erg_marquardt.h
>> >> >
>> >> >
>> >> > http://code.google.com/p/libmv/source/browse/trunk/src/libmv/numeric/l
>> >> >evenb 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/d
>> >> >ogleg .h
>> >> >
>> >> >
>> >> > http://code.google.com/p/libmv/source/browse/trunk/src/libmv/numeric/d
>> >> >ogleg _test.cc
>> >> >
>> >> > Keir
>> >>
>> >> --
>> >> Timothy Hunter
>> >> Student (Stanford University)
>> >>
>> >> T. 404 421 3075
>
> --
> Timothy Hunter
> Master Student
> Electrical Engineering
> Stanford University
> http://ai.stanford.edu/~tjhunter/
> T. 404 421 3075