Newbie trying to use EvaluationCallback with AutoDiff

[delorytheape]

Hi,
Thanks for creating ceres, it looks like a good/useful library. I am a newbie trying to understand how to structure code to solve my problem, and was wondering if someone could point to an example etc.. that may be readily adapted to my use-case. My problem is as follows;

I wish to fit a non-linear, parameterised curve f(x; a_1, a_2, a_3, a_4, a_5) to an observed dataset y(x), where a_1...a_5 are parameters to be estimated.
I am able to implement f(x; a_1, a_2, a_3, a_4, a_5) with templates.
The 'curve' is mathematically complicated, and difficult to determine closed-form expressions for the derivatives, so am hoping to use auto-derivatives.
There are several thousand observations (eg. values of y(x) for x in [-2000,...2000] ) with which to fit f(x; a_1, a_2, a_3, a_4, a_5) .
It is much more efficient to evaluate the function for all values of x at the same time (rather than evaluating individual points in isolation).

As I currently understand (after first reading of the documentation), the best way forward would be similar to ./examples/evaluation_callback_example.cc. ie;

• Implement a class derived from EvaluationCallback (to efficiently evaluate the function at all values of 'x' at the same time)
• Implement a 'CostAndJacobianCopyingCostFunction' to return residuals and Jacobians evaluated at individual values of 'x'.

The part I am not quite sure of is;
Q: Can I use ceres::AutoDiffCostFunction in the EvaluationCallback function to avoid having to calculate derivatives by hand? (The simple example has a closed form expression for the derivatives, but my function is too complicated for that...)

[sandwichmaker]

Have you tried doing something much simpler, which is having a single residual block per observation and then using multithreading so that ceres can use multiple threads to do the cost and Jacobian computation?

[delorytheape]

Thanks for your reply sandwichmaker. I have tried various approaches (unsuccessfully), including a single residual block per observation (eg. as per the ./example/evaluation_callback_example.cc).

I have done a deep-dive into the implementation of auto-derivatives, and am currently pursuing a solution that (hopefully) will allow me to use auto-derivatives (with my own custom extensions).

[sandwichmaker]

What do you mean by these approaches being unsuccessful. Can you share the output of Summary::FullReport when you had one residualblock per observation?

[delorytheape]

By 'unsuccessful' I mean that I could not see how to structure my code to be able to use autodifferentiation within an evaluation_callback function. The only example of a callback (evaluation_callback_example.cc) assumes a closed form of the derivative is available - and as a newbie I am trying to work out if this is a limitation of Ceres, or of my understanding...

As an aside, I am also trying to determine if Ceres autodifferentiation can handle complex numbers in the calculation of residuals (wanting to use a templated FFT/IFFT library)?

[delorytheape]

Hi @sandwichmaker, thanks again for your reply, and apologies for my slow response.

Further context: I am ultimately trying to target an embedded processor, so multithreading is not really a viable option - I need to structure the code to minimize calculations, memory footprint, etc as much as possible - hence my motivation for calculating all residuals once (and with one thread).

I have attempted to put together a "MWE" to show what I am trying to do. Note: This MWE is not the actual problem I am trying to solve, but is structured in a similar way. It is the software structure that I am interested in (rather than a closed form solution to this MWE).

Consider an unknown pulse in the form f(t) = aexp(-bt)+c; whose parameters (a,b,c) we wish to estimate. We do not have data samples of f(t) itself, but we do have access to noisy values of the sampled autocorrelation function R(\tau) = IFFT( || FFT(f(n)) ||^2 ), where FFT(f(n)) represents the fast fourier transform of a sampled version of f(t), and IFFT(f) is the inverse FFT. We have several thousand samples of the autocorrelation function R(\tau)

Option 1: (Following examples/curve_fitting.cc paradigm); Create a vector residual functor something like the following;

struct vector_residual_functor
{ 
    vector_residual_functor( const double*& sample_autocorrelation_function, const int & acf_index)
    {
        // Pseudo-code
        local_acf = sample_autocorrelation_function; // Store the autocorrelation function locally for residual computing later.
        tau_value = acf_index;
    }
    
    template <typename T>
    bool operator()(const T* const a, const T* const b, const T* const c, T* residuals) const
    { 
        // Pseudocode to;
        //     i) Using parameters (a,b,c),  reconstruct estimate of f(t) at 4096 points, then 
        //     ii) calculate estimate of observed autocorrelation using FFT & IFFT
        std::vector<T> f_hat(4096, 0.0);                  // Vector for storing estimate of f(t)
        std::vector<double>  t = LinearSpace(0, 4096);  // Vector containing (0,1,2,3,.....,4095)
        
        // Estimate pulse shape
        f_hat = a*exp(-b*t)+c;    // Assume this works and we get f(t) evaluated at 4096 points..
        
        // Pseudo code to estimate autocorrelation function.
        // Assume that conj(z) calculates the complex conjugate of  vector 'z'
        // And that I have created a class for manipulating complex numbers that contain Jets....
        std::vector<T> acf_hat =  IFFT(FFT(f_hat) * conj(FFT(f_hat)) );  
        
        // Calculate and return residual for this point
        residuals[0] = local_acf[tau_value] - acf_hat[tau_value];
        return true;
    }
    
    std::vector<double> local_acf;  // Local version of autocorrelation function.
    int tau_value;                               // Indicates  autocorrelation lag for this struct instance.. 
}

Then create lots of residual blocks in the 'Problem' (one for each observation)

....etc...

double a,b,c = 0.0;

Problem problem;
for (int ii = 0; ii < 4096; ii++) {
  problem.AddResidualBlock(  
          new AutoDiffCostFunction<vector_residual_functor, 1, 1, 1, 1>(
                 new vector_residual_functor(autocorrelation_data, ii));
          , nullptr, 
          &a,
          &b, 
          &c);
}

In my understanding, the problem with this structure is the very inefficient way the residual is calculated. The auto-correlation function estimate O(N logN) is computed 'N' (=4096) times (once for each residual block), making the whole approach O(N^2 log N). In reality, the auto-correlation only needs to be calculated once per residual update cycle, and the appropriate value (depending on the index tau_value) used to update the appropriate residual block. This motivates a structure similar to the original post, which is based on examples/evaluation_callback_example.cc ie;

Option 2:

• Implement a class derived from EvaluationCallback to efficiently estimate autocorrelation for all residual blocks with one calculation (eg. PrepareForEvaluation(...) )
• Implement a 'CostAndJacobianCopyingCostFunction' to return residuals and Jacobians evaluated at individual values of 'tau'.

The question I have for the second approach is; Q: Can I use ceres::AutoDiffCostFunction in the EvaluationCallback function to avoid having to calculate derivatives by hand?

[sandwichmaker]

Thanks for sharing the details of your problem. A couple of thoughts:

1. I agree computing the autocorrelation function in each residual is not a great idea. If you are going to do autodiff for all this anyways, then you should just put all the residuals in a single functor that computes a multidimensional residual. That will give you the same or better performance than doing the EvaluationCallback.
2. If you really want peak performance, then you should do analytic derivatives. You will likely be able to use analytic properties of the FFT to make mathematical shortcuts which autodiff will never be able to. This sort of thing comes up for example when you are trying to differentiate the solution to a linear system Ax = b, w.r.t say the coefficient matrix A. Trying to differentiate the linear solver is very expensive, it is simpler to analytically derive the derivative of x w.r.t A and use it, which turns out to be quite cheap once you have done to factor A to solve for x anyways.
3. To answer your question about using AutoDiffCostFunction inside your EvaluationCallback, sure you can create whatever objects you want to create and call them in the EvaluationCallback. You are responsible for managing their lifetime etc. But I do not see how this would be any better than the approach I have mentioned in (1).

[delorytheape]

Thanks again for your perspective. Much appreciated.

[sandwichmaker]

My suggestion to you is (especially since you are new to ceres) is to start by not using the evaluationcallback, and just do one residualblock per observation and enable threading and see what kind of performance you get. I suspect you will get pretty good performance with that and not have to use evaluationcallback.