[Insight-users] Gradient computation bug in SingleValuedVnlCostFunctionAdaptor when scales are used (patch included)

Thu Apr 24 13:25:49 EDT 2008

Hi Luis, Tom, others,

Related to this:

I have written a optimizer base class that handles the scaling 
functionality of optimizers. By letting all itk optimizers inherit from 
this class, the scales are dealt with in exactly the same manner. In fact, 
the code is similar to the itkSingleValuedVnlCostFunctionAdaptor.cxx 
mechanism, including the proposed modification of Tom.

You may have a look at the doxygen documentation:

http://elastix.isi.uu.nl/doxygen/a00554.html
which describes the ScaledCostFunction class

http://elastix.isi.uu.nl/doxygen/a00555.html
which described the ScaledOptimizer class

I've programmed quite some itk::optimizers that inherit from this base 
class. On the site mentioned above you can find links to the inheriting 
optimizers, including a quasi-Newton LBFGS optimizer (which does exactly 
the same as the current vnl-based-version, but is written entirely in 
ITK-style) and a conjugate gradient optimizer, which implements various 
types of conjugate gradient (polakribier, fletcher reeves, dai-yuan etc). 
Also a derivative-free evolution-strategy-based optimizer (CMA-ES) is 
implemented.

It is part of a larger piece of software for image registration. You could 
inspect the source code at http://elastix.isi.uu.nl/download.php .

Kind regards,
Stefan.

At 12:32 24-4-2008, you wrote:

>Hi Tom,
>
>Thanks for pointing this out.
>
>You are right, the current code doesn't perform the correct
>mathematical operation that corresponds to scaling the argument
>of a cost function.
>
>
>However, that was not the intent of introducing the scaling
>in these classes.
>
>
>The motivation for this scaling operation was to extend the
>functionality of VNL optimizers to be usable with domains outside
>of the cannonical [-1:+1] domain.
>
>
>Since these functions are used for the purpose of optimization,
>we intended simply to shrink (or expand) the domain of the original
>problem in order to map to the [-1:+1] that the algorithmic
>implementation of vnl optimizers expect.
>
>To make sense of it, you may find interesting to look at both the
>ITK wrappers for VNL optimizer, and the ITK wrapper for the cost
>function.
>
>
>Take for example,
>the itkAmoebaOptimizer.cxx:
>
>   it its call to GetValue() in line 84, it uses the scales
>   to multiply the parameters *before* they are passes to
>   the cost function.
>
>      84:      parameters[i] *= scales[i];
>
>   In its call to StartOptimization(), before invoking the
>   VNL optimizer, we scale the parameters of the transform
>
>     218:       parameters[i] *= scales[i];
>
>   then we invoke the VNL optimizer
>
>     227:      m_VnlOptimizer->minimize( parameters );
>
>    and then we restore the scale of the transform parameters
>
>    242:       parameters[i] /= scales[i];
>
>
>The goal, is to create a situation where the VNL optimizer
>"thinks" that is optimizing a function whose domain is close
>to [-1:+1], as VNL algorithms expect or assume.
>
>
>However, the VNL optimizer is calling the cost function
>in order to get values from it. That cost function will
>happen to be an ITK ImageMetric (most of the time), and
>therefore the VNL cost function adaptor must "restore"
>the scale of the parametric space to what the metric
>expects as parameters. That's why the cost function adaptor
>
>     itkSingleValuedVnlCostFunctionAdaptor.cxx
>
>in its call to "f()" rescales the parameters *back* by using
>in line 61:
>
>    61:   parameters[i] = inparameters[i]/m_Scales[i];
>
>and then calling the actual ITK cost function in line
>
>    69:   m_CostFunction->GetValue( parameters ));
>
>
>In this way, ITK metrics will still work in a space where
>the domain of translations may be [-500mm, +500mm] and
>[-0.1degrees,+0.1degrees], while the VNL optimizer works
>in a space with a domain close to [-1:+1].
>
>
>The point is well taken in that we could have applied the
>scaling to the resulting values of the gradients. Something
>that could conveniently be done in the method:
>
>        ConvertExternalToInternalGradient()
>
>in line 176 of itkSingleValuedVnlCostFunctionAdaptor.cxx
>
>
>This additional (and mathematically correct scaling) will
>indicate the following:
>
>
>    A) if we have an image of 100x100 pixels with spacing
>       2.0mm x 2.0mm and we are optimizing the Mean square
>       metric against another image of similar dimensions,
>       we see this problem happening in a domain of
>       [200x200]mm.  If the Metric changes by 35 intensity
>       values in average when you translate the images by
>       5 millimeters, that corresponds to a derivative of
>
>                        7 iv / mm
>
>
>    B) The scaling conversion will be equivalent to replace
>       the units of mm with units of "image width".  That is
>       somehow normalizing the space to the size of the image.
>
>       In this context, the conversion
>
>                F(x) = f(x/s)
>
>                grad(F(x)) = grad(f(x/s))/s
>
>       will give us the value of the gradient on the units
>       of "image width", as:
>
>                     (7/200)  iv / iw
>
>
>
>      At this point, we could make this correction,
>      but this may require to recalibrate the parameters
>      of most of image registration examples and tests
>      that use VNL optimizers...
>
>
>Have you run an experimental build with this change
>to evaluate the impact of the modification ?
>
>
>----
>
>
>The operation is not done in the same way either across optimizers.
>
>For example,
>if you look at the GradientDescent and RegularStepGradientDescent
>optimizers, the scales are applied only to the *resulting* gradient
>vector returned from the cost function, and then, the resulting vector
>is used during the computation of the step to be taken. In these two
>optimizers, the scale factor is not used for feeding the values
>inside the cost function.
>
>
>
>     Regards,
>
>
>         Luis
>
>
>-----------------------
>Tom Vercauteren wrote:
>>Hi all,
>>I think I found a bug in
>>SingleValuedVnlCostFunctionAdaptor::ConvertExternalToInternalGradient.
>>Here is a brief description:
>>SingleValuedVnlCostFunctionAdaptor can be used to wrap an itk cost
>>function into a vnl cost function.
>>Additionally the user can specify a set of scales such that the
>>function to be optimized is F(x) = f(x./s) instead of f(x).
>>The gradient of F, grad(F)(x) is thus given by grad(F)(x) = grad(f)(x./s)./s
>>Instead of that, the cost function adpator returns grad(f)(x./s)
>>
>>I have filed a bug report in the bug tracker:
>>http://www.itk.org/Bug/view.php?id=6887
>>Let me know if I am missing something.
>>Best regards,
>>Tom Vercauteren
>>_______________________________________________
>>Insight-users mailing list
>>Insight-users at itk.org
>>http://www.itk.org/mailman/listinfo/insight-users
>_______________________________________________
>Insight-users mailing list
>Insight-users at itk.org
>http://www.itk.org/mailman/listinfo/insight-users

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20080424/8841801a/attachment-0001.htm>