scipt.opt.approx_fprime not good enough

[samcoveney]

I am checking through some of the tests for reconst/dti, and I am finding (after some changes made for an up and coming PR) that the checks on the non-linear-least-squares error function derivate using scipy.opt.approx_fprime are failing.

The tests assert that the analytical Jacobian and the numerical Jacobian are not equal.

This is not the first time I have encountered this issue (I came across it while doing this #2730 as well, and have found the same problem in other projects as well).

I decided to quickly check if approximation of the derivatives were the problem, so I installed some autograd code and got to work in the tests in dipy\reconst\tests\test_dti.py (look for # === SEE HERE ===):

    # === SEE HERE: import autograd libraries ===
    import autograd.numpy as npa  # Thinly-wrapped numpy
    from autograd import grad as grd    # The only autograd function you may ever need

    nlls = dti._NllsHelper()
    sigma = 1.4826 * np.median(np.abs(error - np.median(error)))
    for D in [D_orig, np.zeros_like(D_orig)]:

        # Test Jacobian at D
        args = [D, X, Y, 1/sigma**2]
        # NOTE: 1. call 'err_func', to set internal stuff in the class
        nlls.err_func(*args)
        # NOTE: 2. call 'jabobian_func', corresponds to last err_func call
        analytical = nlls.jacobian_func(*args)

        for i in range(len(X)):

            args = [X[i], Y[i], 1/sigma**2]

            # === SEE HERE: commented out ===
            #approx = opt.approx_fprime(D, nlls.err_func, 1e-8, *args)

            # === SEE HERE: use autodiff to approximate the derivative ===
            grad_err_func = grd(nlls.err_func) 
            approx = grad_err_func(D, X[i], Y[i], 1/sigma**2)

            print(approx - analytical[i])
            assert np.allclose(approx, analytical[i])

The tests now pass because, of course, the function nlls.jacobian_func really is the correct jacobian of nlls.err_func... I am printing the difference just to check, and the autograd perfectly performs the differentiation.

I suspect the problem is with scipy.opt.approx_fprime using only forward differences.

Suggestions welcome please, I am happy to attempt a fix as part of the PR I am finishing up. We could switch these tests off, but given that my next PR has involved making changes to this function, it would be good to have a working test.

(edit: I had to do import autograd.numpy as np in dipy\reconst\dti.py for this autograd to work)

[samcoveney]

Any thoughts @skoudoro @arokem ?

[skoudoro]

Any thoughts @skoudoro @arokem ?

not yet 😅 , I need to take the time to get deeper because there is no easy answer.

[samcoveney]

I've turned off the test in #3170 since I've already verified the gradient. Can change later if required.

One issue is that the location (in parameter space) at which the gradient is evaluated is the minimum of the cost function (sum of squares), which means that a forward difference approximation at this point is going to be fairly poor. If the location of testing the gradient was not the minimum of the cost function, perhaps this might work better... I will have a look

[tanishka321]

Hi, @samcoveney Can I look at this issue?
I am not sure but I think considering the epsilon=1e-6 parameter in the approx_fprime function may help as reducing the step size (epsilon) in finite difference approximation decreases the spacing between the points used to calculate the gradient, making it behave more like central differences.

Line 584: approx = opt.approx_fprime(D, nlls.err_func, 1e-8, *args, epsilon=1e-6)

I am stuck with very silly doubt on how to compare and test it 😅, Could you please guide me on this and how to proceed further on this approach?

[samcoveney]

I have tested different epsilon - I am now quite sure that the problem is testing using first-order finite difference approximation at the minimum, where it is not a good approx (I used to do many CUDA simulations using finite differences).

The fact that the auto-grad gave an exact match seems to prove to me that first order finite difference is the issue, since a different numerical calculation of the gradient gave an exact match - remember, we are not just dealing with 1 variable here, but with 7.

If you have time, what you could try it changing D to not be the same as the "ground truth D", so that we are not approximating the derivative at the minimum? Perhaps paired with another epsilon, if needed. Other than that, I am a bit stumped as to what else to do.

To compare and test, just change epsilon and D and run the test, something like:

coverage run -m pytest -s --doctest-modules --verbose ~/Software/dipy_fork/dipy/reconst/tests/test_dti.py -k 'test_nnls_jacobian_func'

I would also suggest to try implementing the autograd option as per my first comment (note the last line of that comment), just to convince you that this is what is going on

[samcoveney]

@tanishka321 any progress?