JuliaDiff / SparseDiffTools.jl

Fast jacobian computation through sparsity exploitation and matrix coloring
MIT License
237 stars 41 forks source link

Different numerical behaviour after updating SparseDiffTools.jl from v1.16.1 to v1.16.3 #151

Closed ferrolho closed 2 years ago

ferrolho commented 2 years ago

Hi! I am using ForwardColorJacCache, matrix_colors, and forwarddiff_color_jacobian! in my project. I am solving an optimisation problem and I use forwarddiff_color_jacobian! to evaluate jacobians using forward-mode automatic differentiation.

The changes to SparseDiffTools.jl introduced in v1.16.2 and v1.16.3 have changed the optimiser's behaviour significantly. Below you can find the iteration logs of the solver before and after the update.

Using v1.16.1:

  Iter      Objective      FeasError   OptError    ||Step||    CGits 
--------  --------------  ----------  ----------  ----------  -------
       0    0.000000e+00   3.891e+02
      10    0.000000e+00   1.453e-03   0.000e+00   3.027e+00        0
      11    0.000000e+00   2.536e-07   0.000e+00   1.626e-02        0

EXIT: Locally optimal solution found.

Using v1.16.3:

  Iter      Objective      FeasError   OptError    ||Step||    CGits 
--------  --------------  ----------  ----------  ----------  -------
       0    0.000000e+00   3.891e+02
      10    0.000000e+00   2.616e+02   0.000e+00   1.198e+02        0
      20    0.000000e+00   2.193e+02   0.000e+00   7.919e+01        0
      30    0.000000e+00   1.973e+02   0.000e+00   1.527e+02        0
      40    0.000000e+00   1.678e+02   0.000e+00   2.000e+02        0
      50    0.000000e+00   1.567e+02   0.000e+00   1.404e+02        0
      60    0.000000e+00   1.486e+02   0.000e+00   5.208e+01        0
      70    0.000000e+00   1.400e+02   0.000e+00   1.450e+02        0
      80    0.000000e+00   1.304e+02   0.000e+00   8.402e+01        0
      90    0.000000e+00   1.248e+02   0.000e+00   7.268e+01        0
     100    0.000000e+00   1.201e+02   0.000e+00   1.946e+02        0
     110    0.000000e+00   1.158e+02   0.000e+00   1.009e+02        0
     120    0.000000e+00   1.128e+02   0.000e+00   1.455e+02        0
     130    0.000000e+00   1.103e+02   0.000e+00   3.402e+01        0
     140    0.000000e+00   1.075e+02   0.000e+00   6.832e+01        0
     150    0.000000e+00   1.052e+02   0.000e+00   9.935e+01        0
     160    0.000000e+00   1.034e+02   0.000e+00   6.049e+01        0
     170    0.000000e+00   1.011e+02   0.000e+00   6.967e+01        0
     180    0.000000e+00   9.893e+01   0.000e+00   1.303e+02        0
     190    0.000000e+00   9.733e+01   0.000e+00   7.816e+01        0
     193    0.000000e+00   9.690e+01   0.000e+00   3.735e+01        0

EXIT: Time limit reached. Current point is infeasible.

I am not sure why this is happening...

ChrisRackauckas commented 2 years ago

Can you isolate this to a Jacobian which is different? I cannot find a case so it's a bit odd, but it could just be a missing case.

ferrolho commented 2 years ago

Yes, sorry. I do recognise that what I posted above is very little info. I am trying to do that now.

I just tested the behaviour in other more simple problems and I didn't see any numerical differences there, so it must be quite specific...

ChrisRackauckas commented 2 years ago

Yeah. I know v1.16.2 had an issue, so maybe you were accidentally testing that? The OrdinaryDiffEq tests picked it up though, and v1.16.3 was the hotfix (and that was rather specific). But v1.16.3 got through all of the DiffEq tests just fine, and that tests dense and sparse, so I would be surprised if there was something still missing. That said, there definitely could such a case.

ferrolho commented 2 years ago

This is something on v1.16.3. When I was rolling back the releases one by one, I did go through v1.16.2, which for my same problem was crashing the solver; but the hotfix changed that.

ferrolho commented 2 years ago

I've found the specific lines causing this.

If I revert https://github.com/JuliaDiff/SparseDiffTools.jl/blob/47a8a9f4725e1ca8b07d5d17a8e674dd2718292d/src/differentiation/compute_jacobian_ad.jl#L320-L326 to https://github.com/JuliaDiff/SparseDiffTools.jl/blob/fb09091c195311ab622ced4889694e889631c980/src/differentiation/compute_jacobian_ad.jl#L287 the behaviour remains the same.

ChrisRackauckas commented 2 years ago

hmm, that looks correct. What if you remove the ivdep?

ferrolho commented 2 years ago

Interesting. Just removing ivdep fixes the issue. I'm not familiar with what that keyword is used for, though.

ChrisRackauckas commented 2 years ago

It should be safe if not aliasing, but apparently there's still a memory connection in here? Do a PR to remove it.

ferrolho commented 2 years ago

If not aliasing what? I've submitted the PR.

ChrisRackauckas commented 2 years ago

two different arrays in the loop

ferrolho commented 2 years ago

Thank you for your help, Chris. Prompt as always! :rackauckas: