4 ms·

Automatic differentiation can be incorrect

74 points by abetusk 1y ago

NooneAtAll3 1y ago
[video]?
ComplexSystems 1y ago
This article is saying that it can be numerically unstable in certain situations, not that it's theoretically incorrect.
- omnicognate 1y ago
  Yeah, perhaps the actual title would be better: "The Numerical Analysis of Differentiable Simulation". (Rather than the subtitle, which is itself a poor rewording of the actual subtitle in the video.)
- goosedragons 1y ago
  It can be both. A mistake in AD primitives can lead to theoretically incorrect derivatives. With the system I use I have run into a few scenarios where edge cases aren't totally covered leading to the wrong result. I have also run into numerical instability too.
  - froobius 1y ago
    > A mistake in AD primitives can lead to theoretically incorrect derivatives Ok but that's true of any program. A mistake in the implementation of the program can lead to mistakes in the result of the program...
    - goosedragons 1y ago
      That's true! But it's also true that any program dealing with floats can run into numerical instability if care isn't taken to avoid it, no? It's also not necessarily immediately obvious that the derivatives ARE wrong if the implementation is wrong.
      - srean 1y ago
        > It's also not necessarily immediately obvious that the derivatives ARE wrong if the implementation is wrong. It's neither full proof or fool proof but an absolute must is a check that the loss function is reducing. It quickly detects a common error that the sign came out wrong in my gradient call. Part of good practice one learns in grad school.
      - froobius 1y ago
        You can pretty concretely and easily check that the AD primatives are correct by comparing them to numerical differentiation.
    - godelski 1y ago
      I haven't watched the video but the text says they're getting like 60+% error on simple linear ODEs which is pretty problematic. You're right, but the scale of the problem seems to be the issue
- adgjlsfhk1 1y ago
  This is less about numerical instability and more that iterative algorithms with error control their error, but when you run AD on them you are ADing the approximation and a derivative of an approximation can be arbitrarily different from an approximation of a derivative.
  - ogogmad 1y ago
    That makes more sense. The title is flat out wrong IMO.
    - adgjlsfhk1 1y ago
      I think it is correct. lots of people view AD as a black box that you can throw algorithms to and get derivatives out, and this shows that that isn't true
      - wakawaka28 1y ago
        If you wrote code that failed to compile, you wouldn't impulsively call your compiler incorrect. This title sounds like it puts the blame in the wrong place. You can get error accumulation from even a basic calculation in a loop. We could try to solve these problems but it's not the algorithm's fault you don't know what you're doing.
        ChrisRackauckas 1y ago
        This has nothing to do with floating point error accumulation or numerical stability in the floating point sense. You can do this with arbitrary sized floating point values and you will still get the same non-convergence result.
        wakawaka28 1y ago
        Yes, it does. I admit that I have not tried to get deep into this, but it says right in the summary that this is due to error propagation and numerical instability. You can have such errors in ordinary hand-written code. The size of the float does not matter (much) for some examples, which wipe out a bunch of significant figures from the result. I'm not going to sit here and pretend I know exactly the details going on here but I studied numerical analysis and AD quite a bit back in the day. Even in the 2010s people knew that there was a chance of the resulting expressions having inherent instability. Even simple things like the quadratic formula have more stable and less stable forms. How much worse might it be for auto-generated expressions (or expression sequences, or equivalent; AD can be done a few ways) such as in AD? AD guarantees analytically correct logic (in infinite precision, for example) if you use it right but error analysis is not even attempted by most libraries.
Lerc 1y ago
With floats getting smaller and smaller in ML, It's hard to imagine anyone failing to learn this as one of their early experiences in the field. The focus should not be on the possibility of error, but managing the error to be within acceptable limits. There's a hour long video there, and it's 3am, so I'm not sure how much of this covers management. Anyone familiar with it care to say?
- Legend2440 1y ago
  Floats are getting smaller and smaller only in neural networks - which works ok because NNs are designed for nice easy backprop. Normalization and skip connections really help avoid numerical instability. The author is talking about ODEs and PDEs, which come from physical systems and don't have such nice properties.
  - storus 1y ago
    Skip connections are about preserving identity, not to make it more numerically stable.
    - D-Machine 1y ago
      I'm not sure your claim is correct, nor particularly precise. Skip connections do propagate feature maps from earlier to later/deeper in the network, and one can model this mathematically as using an identity function to do this, but I wouldn't say that this is "about preserving identity". Also, the reason the identity is useful here is because it prevents vanishing gradients (due to the identity function having an identity tensor as the gradient). Arguably, preventing vanishing gradients is exactly about numerical stability, so I would in fact argue that skip connections are there to make gradient descent of certain architectures more numerically stable. EDIT: Just to provide some research support: https://arxiv.org/html/2405.01725v1#S5.SS4.5.1.1 https://arxiv.org/html/2405.01725v1#S5.SS4.5.1.1. I.e. skip connections do a lot of things, but most of these things boil down to making gradients better conditioned. Insofar as as one problem of poor conditioning is gradient vanishing, and vanishing is really only a problem because of the limited precision of floating point, it is hard for me to not think skip connections are primarily there for numerical stability, broadly. One case supporting your point, though, where one might argue the "main" reason for skip connections is propagating information rather than for stability, is in LONG skip connections, as in U-Nets. There, since feature maps shrink with network depth, long skip connections are needed to bring high-resolution information from earlier layers to deeper ones.
deckar01 1y ago
I have been using sympy while learning electron physics to automatically integrate linear charge densities. It works great symbolically, but often fails silently when the symbols are substituted with floats before integration. https://github.com/sympy/sympy/issues/27675 https://github.com/sympy/sympy/issues/27675
ismailmaj 1y ago
The way I explained it to myself in the past why so much of the CUDA algorithms don't care much about numerical stability is that the error is a form of regularization (i.e. less overfitting over the data) in deep learning.
- Nevermark 1y ago
  I am not quite sure what that means! :) But reasons why deep learning training is very robust to moderate inaccuracy in gradients: 1. Locally, sigmoid and similar functions are the simplest smoothest possible non-linearity to propagate gradients through. 2. Globally, outside of deep recurrent networks, there is no recursion which makes the total function smooth and well behaved. 3. While the perfect gradient indicates the ideal direction to adjust parameters, for fastest improvement, all that is really needed to reduce error is to move parameters in the direction of the gradient signs, with a small enough step. That is a very low bar. It's like telling an archer they just need to shoot an arrow so it lands closer to the target than where the archer is standing, but not worry about hitting it! 4. Finally, the perfect first order gradient is only meaningful at one point of the optimization surface. Moving away from that point, i.e. updating the parameters at all, and the gradient changes quickly. So we are in gradient heuristic land even with "perfect" first order gradients. The most perfectly calculated gradient isn't actually "accurate" already. To actually get an accurate gradient over a parameter step, would take fitting the local gradient with a second or third order polynomial. I.e. not just first, but second and third order derivatives. At vastly greater computational and working memory cost. -- The only critical issue for calculating gradients, is that there is enough precision that at least directional gradient information makes it from errors back to the parameters to update. If precision is too low, then the variable magnitude rounding inherent to floating point arithmetic can completely drop directional information for smaller gradients. Without accurate gradient signs, learning stalls.
  - ismailmaj 1y ago
    Typically for matrix multiplications there is a wide range of algorithms you could use to compute it, on one extreme end you could use numerically stable summation and the other extreme you could have tiled matmul with FP8, the industry trend seems to go further away from numerical stable algorithms without much quality drop it seems. My claim is probably unfair since it ignores the scale you gain from the speed/precision tradeoff, so I assumed numerical stability is not that beneficial compared to something precision heavy like physics simulation in HPC.
ChrisRackauckas 1y ago
Interesting to see this here! This example is one of the ones mentioned in the appendix of https://arxiv.org/abs/2406.09699 https://arxiv.org/abs/2406.09699, specifically "4.1.2.4 When AD is algorithmically correct but numerically wrong". If people want a tl;dr, the main idea is you can construct an ODE where the forward pass is trivial, i.e. the ODE solution going forwards is exact, but its derivative is "hard". An easy way to do this is to make it so you have for example `x' = x - y, y' = y - x`, with initial conditions x(0)=y(0). If you start with things being the same value, then the solution to the ODE is constant since `x' = y' = 0`. But the derivative of the ODE solution with respect to its initial condition is very non-zero: a small change away from equality and boom the solution explodes. You write out the expression for dy(t)/dy(0) and what you get is a non-trivial ODE that has to be solved. What happens in this case though is that automatic differentiation "runs the same code as the primal case". I.e., automatic differentiation has the property that it walks through your code and differentiates the steps of your code, and so the control flow always matches the control flow of the non-differentiated code, slapping the derivative parts on each step. But here the primal case is trivial, so the ODE solver goes "this is easy, let's make dt as big as possible and step through this easily". But this means that `dt` is not error controlled in the derivative (the derivative, being a non-trivial ODE, needs to have a smaller dt and take small steps in order to get an accurate answer), so the derivative then gets error due to this large dt (or small number of steps). Via this construction you can make automatic differentiation give as much error as you want just by tweaking the parameters around. Thus by this construction, automatic differentiation has no bound to the error it can give, and no this isn't floating point errors this is strictly building a function that does not converge to the correct derivative. It just goes to show that automatic differentiation of a function that computes a nice error controlled answer does not necessarily give a derivative that also has any sense of error control, that is a property that has to be proved and is not true in general. This of course is a bit of a contrived example to show the point that the error is unbounded, but then it points to real issues that can show up in user code (in fact, this example was found because a user opened an issue with a related model). Then one thing that's noted in here too is that the Julia differential equation solvers hook into the AD system to explicitly "not do forward-mode AD correctly", incorporating the derivative terms into the time stepping adaptivity calculation, so that it is error controlled. The property that you get is that for these solvers you get more steps to the ODE when running it in AD mode than outside of AD mode, and that is a requirement if you want to ensure that the user's tolerance is respected. But that's explicitly "not correct" in terms of what forward-mode AD is, so calling forward-mode AD on the solver doesn't quite do forward mode AD of the solver's code "correctly" in order to give a more precise solution. That of course is a choice, you could instead choose to follow standard AD rules in such a code. The trade-off is between accuracy and complexity.
drnick1 1y ago
The authors seem to push a Julia implementation, but the issue is that no one uses Julia...
- eigenspace 1y ago
  The context here is that the author is the principal developer of one of the biggest (if not the biggest) and most performant differential equation solver suites in the world[1]. It's written in julia, and it is very widely used, and known to be state of the art. [1] https://github.com/SciML/DifferentialEquations.jl https://github.com/SciML/DifferentialEquations.jl
  - drnick1 1y ago
    Would it not make more sense to write such a tool in C? The idea here is that C is universal, and higher level languages like Python, R, Julia and others all have a foreign function interface for C.
    - eigenspace 1y ago
      Good luck being that productive in C, distributing the software adequately, and attracting collaborators.
    - adgjlsfhk1 1y ago
      No. C isn't an expressive enough language to write many solvers in. Many of the best solvers are basically impossible to write without meta-programming (where you write a generator for the n-stage version that is then able to generate the code for all the possible permuations), and there are also massive performance gains through symbolic-numeric approaches where the solver does introspection of the user function to rearrangements that aren't possible in a black-box framework.
    - ChrisRackauckas 1y ago
      This tool already exists in Python and R. See https://cran.r-project.org/web/packages/diffeqr/index.html https://cran.r-project.org/web/packages/diffeqr/index.html and https://anaconda.org/conda-forge/diffeqpy https://anaconda.org/conda-forge/diffeqpy. Julia has language bindings that make it simple enough to bind to other high level languages that these projects are maintained by the core team and supports lots of the library, including forms of automatic differentiation and GPU kernel generation. See for example the bindings that allow for usage within PyTorch https://github.com/SciML/juliatorch https://github.com/SciML/juliatorch and the Python-based Collab notebooks showing the GPU usage https://colab.research.google.com/drive/1bnQMdNvg0AL-LyPcXBiH10jBij5QUmtY?usp=sharing https://colab.research.google.com/drive/1bnQMdNvg0AL-LyPcXBi.... With Julia v1.12's small binary generation, we plan to release forms via binaries with C ABIs over the next year as well. Are those sufficient or should we consider supporting other deployments?

74 points by abetusk 1y ago

[video]?

This article is saying that it can be numerically unstable in certain situations, not that it's theoretically incorrect.

omnicognate 1y ago
Yeah, perhaps the actual title would be better: "The Numerical Analysis of Differentiable Simulation". (Rather than the subtitle, which is itself a poor rewording of the actual subtitle in the video.)
goosedragons 1y ago
It can be both. A mistake in AD primitives can lead to theoretically incorrect derivatives. With the system I use I have run into a few scenarios where edge cases aren't totally covered leading to the wrong result. I have also run into numerical instability too.
- froobius 1y ago
  > A mistake in AD primitives can lead to theoretically incorrect derivatives Ok but that's true of any program. A mistake in the implementation of the program can lead to mistakes in the result of the program...
  - goosedragons 1y ago
    That's true! But it's also true that any program dealing with floats can run into numerical instability if care isn't taken to avoid it, no? It's also not necessarily immediately obvious that the derivatives ARE wrong if the implementation is wrong.
    - srean 1y ago
      > It's also not necessarily immediately obvious that the derivatives ARE wrong if the implementation is wrong. It's neither full proof or fool proof but an absolute must is a check that the loss function is reducing. It quickly detects a common error that the sign came out wrong in my gradient call. Part of good practice one learns in grad school.
    - froobius 1y ago
      You can pretty concretely and easily check that the AD primatives are correct by comparing them to numerical differentiation.
  - godelski 1y ago
    I haven't watched the video but the text says they're getting like 60+% error on simple linear ODEs which is pretty problematic. You're right, but the scale of the problem seems to be the issue
adgjlsfhk1 1y ago
This is less about numerical instability and more that iterative algorithms with error control their error, but when you run AD on them you are ADing the approximation and a derivative of an approximation can be arbitrarily different from an approximation of a derivative.
- ogogmad 1y ago
  That makes more sense. The title is flat out wrong IMO.
  - adgjlsfhk1 1y ago
    I think it is correct. lots of people view AD as a black box that you can throw algorithms to and get derivatives out, and this shows that that isn't true
    - wakawaka28 1y ago
      If you wrote code that failed to compile, you wouldn't impulsively call your compiler incorrect. This title sounds like it puts the blame in the wrong place. You can get error accumulation from even a basic calculation in a loop. We could try to solve these problems but it's not the algorithm's fault you don't know what you're doing.
      - ChrisRackauckas 1y ago
        This has nothing to do with floating point error accumulation or numerical stability in the floating point sense. You can do this with arbitrary sized floating point values and you will still get the same non-convergence result.
        wakawaka28 1y ago
        Yes, it does. I admit that I have not tried to get deep into this, but it says right in the summary that this is due to error propagation and numerical instability. You can have such errors in ordinary hand-written code. The size of the float does not matter (much) for some examples, which wipe out a bunch of significant figures from the result. I'm not going to sit here and pretend I know exactly the details going on here but I studied numerical analysis and AD quite a bit back in the day. Even in the 2010s people knew that there was a chance of the resulting expressions having inherent instability. Even simple things like the quadratic formula have more stable and less stable forms. How much worse might it be for auto-generated expressions (or expression sequences, or equivalent; AD can be done a few ways) such as in AD? AD guarantees analytically correct logic (in infinite precision, for example) if you use it right but error analysis is not even attempted by most libraries.

Lerc 1y ago

With floats getting smaller and smaller in ML, It's hard to imagine anyone failing to learn this as one of their early experiences in the field. The focus should not be on the possibility of error, but managing the error to be within acceptable limits. There's a hour long video there, and it's 3am, so I'm not sure how much of this covers management. Anyone familiar with it care to say?

Legend2440 1y ago
Floats are getting smaller and smaller only in neural networks - which works ok because NNs are designed for nice easy backprop. Normalization and skip connections really help avoid numerical instability. The author is talking about ODEs and PDEs, which come from physical systems and don't have such nice properties.
- storus 1y ago
  Skip connections are about preserving identity, not to make it more numerically stable.
  - D-Machine 1y ago
    I'm not sure your claim is correct, nor particularly precise. Skip connections do propagate feature maps from earlier to later/deeper in the network, and one can model this mathematically as using an identity function to do this, but I wouldn't say that this is "about preserving identity". Also, the reason the identity is useful here is because it prevents vanishing gradients (due to the identity function having an identity tensor as the gradient). Arguably, preventing vanishing gradients is exactly about numerical stability, so I would in fact argue that skip connections are there to make gradient descent of certain architectures more numerically stable. EDIT: Just to provide some research support: https://arxiv.org/html/2405.01725v1#S5.SS4.5.1.1 https://arxiv.org/html/2405.01725v1#S5.SS4.5.1.1. I.e. skip connections do a lot of things, but most of these things boil down to making gradients better conditioned. Insofar as as one problem of poor conditioning is gradient vanishing, and vanishing is really only a problem because of the limited precision of floating point, it is hard for me to not think skip connections are primarily there for numerical stability, broadly. One case supporting your point, though, where one might argue the "main" reason for skip connections is propagating information rather than for stability, is in LONG skip connections, as in U-Nets. There, since feature maps shrink with network depth, long skip connections are needed to bring high-resolution information from earlier layers to deeper ones.

deckar01 1y ago

I have been using sympy while learning electron physics to automatically integrate linear charge densities. It works great symbolically, but often fails silently when the symbols are substituted with floats before integration. https://github.com/sympy/sympy/issues/27675 https://github.com/sympy/sympy/issues/27675

ismailmaj 1y ago

The way I explained it to myself in the past why so much of the CUDA algorithms don't care much about numerical stability is that the error is a form of regularization (i.e. less overfitting over the data) in deep learning.

Nevermark 1y ago
I am not quite sure what that means! :) But reasons why deep learning training is very robust to moderate inaccuracy in gradients: 1. Locally, sigmoid and similar functions are the simplest smoothest possible non-linearity to propagate gradients through. 2. Globally, outside of deep recurrent networks, there is no recursion which makes the total function smooth and well behaved. 3. While the perfect gradient indicates the ideal direction to adjust parameters, for fastest improvement, all that is really needed to reduce error is to move parameters in the direction of the gradient signs, with a small enough step. That is a very low bar. It's like telling an archer they just need to shoot an arrow so it lands closer to the target than where the archer is standing, but not worry about hitting it! 4. Finally, the perfect first order gradient is only meaningful at one point of the optimization surface. Moving away from that point, i.e. updating the parameters at all, and the gradient changes quickly. So we are in gradient heuristic land even with "perfect" first order gradients. The most perfectly calculated gradient isn't actually "accurate" already. To actually get an accurate gradient over a parameter step, would take fitting the local gradient with a second or third order polynomial. I.e. not just first, but second and third order derivatives. At vastly greater computational and working memory cost. -- The only critical issue for calculating gradients, is that there is enough precision that at least directional gradient information makes it from errors back to the parameters to update. If precision is too low, then the variable magnitude rounding inherent to floating point arithmetic can completely drop directional information for smaller gradients. Without accurate gradient signs, learning stalls.
- ismailmaj 1y ago
  Typically for matrix multiplications there is a wide range of algorithms you could use to compute it, on one extreme end you could use numerically stable summation and the other extreme you could have tiled matmul with FP8, the industry trend seems to go further away from numerical stable algorithms without much quality drop it seems. My claim is probably unfair since it ignores the scale you gain from the speed/precision tradeoff, so I assumed numerical stability is not that beneficial compared to something precision heavy like physics simulation in HPC.

ChrisRackauckas 1y ago

Interesting to see this here! This example is one of the ones mentioned in the appendix of https://arxiv.org/abs/2406.09699 https://arxiv.org/abs/2406.09699, specifically "4.1.2.4 When AD is algorithmically correct but numerically wrong". If people want a tl;dr, the main idea is you can construct an ODE where the forward pass is trivial, i.e. the ODE solution going forwards is exact, but its derivative is "hard". An easy way to do this is to make it so you have for example `x' = x - y, y' = y - x`, with initial conditions x(0)=y(0). If you start with things being the same value, then the solution to the ODE is constant since `x' = y' = 0`. But the derivative of the ODE solution with respect to its initial condition is very non-zero: a small change away from equality and boom the solution explodes. You write out the expression for dy(t)/dy(0) and what you get is a non-trivial ODE that has to be solved. What happens in this case though is that automatic differentiation "runs the same code as the primal case". I.e., automatic differentiation has the property that it walks through your code and differentiates the steps of your code, and so the control flow always matches the control flow of the non-differentiated code, slapping the derivative parts on each step. But here the primal case is trivial, so the ODE solver goes "this is easy, let's make dt as big as possible and step through this easily". But this means that `dt` is not error controlled in the derivative (the derivative, being a non-trivial ODE, needs to have a smaller dt and take small steps in order to get an accurate answer), so the derivative then gets error due to this large dt (or small number of steps). Via this construction you can make automatic differentiation give as much error as you want just by tweaking the parameters around. Thus by this construction, automatic differentiation has no bound to the error it can give, and no this isn't floating point errors this is strictly building a function that does not converge to the correct derivative. It just goes to show that automatic differentiation of a function that computes a nice error controlled answer does not necessarily give a derivative that also has any sense of error control, that is a property that has to be proved and is not true in general. This of course is a bit of a contrived example to show the point that the error is unbounded, but then it points to real issues that can show up in user code (in fact, this example was found because a user opened an issue with a related model). Then one thing that's noted in here too is that the Julia differential equation solvers hook into the AD system to explicitly "not do forward-mode AD correctly", incorporating the derivative terms into the time stepping adaptivity calculation, so that it is error controlled. The property that you get is that for these solvers you get more steps to the ODE when running it in AD mode than outside of AD mode, and that is a requirement if you want to ensure that the user's tolerance is respected. But that's explicitly "not correct" in terms of what forward-mode AD is, so calling forward-mode AD on the solver doesn't quite do forward mode AD of the solver's code "correctly" in order to give a more precise solution. That of course is a choice, you could instead choose to follow standard AD rules in such a code. The trade-off is between accuracy and complexity.

drnick1 1y ago

The authors seem to push a Julia implementation, but the issue is that no one uses Julia...

eigenspace 1y ago
The context here is that the author is the principal developer of one of the biggest (if not the biggest) and most performant differential equation solver suites in the world[1]. It's written in julia, and it is very widely used, and known to be state of the art. [1] https://github.com/SciML/DifferentialEquations.jl https://github.com/SciML/DifferentialEquations.jl
- drnick1 1y ago
  Would it not make more sense to write such a tool in C? The idea here is that C is universal, and higher level languages like Python, R, Julia and others all have a foreign function interface for C.
  - eigenspace 1y ago
    Good luck being that productive in C, distributing the software adequately, and attracting collaborators.
  - adgjlsfhk1 1y ago
    No. C isn't an expressive enough language to write many solvers in. Many of the best solvers are basically impossible to write without meta-programming (where you write a generator for the n-stage version that is then able to generate the code for all the possible permuations), and there are also massive performance gains through symbolic-numeric approaches where the solver does introspection of the user function to rearrangements that aren't possible in a black-box framework.
  - ChrisRackauckas 1y ago
    This tool already exists in Python and R. See https://cran.r-project.org/web/packages/diffeqr/index.html https://cran.r-project.org/web/packages/diffeqr/index.html and https://anaconda.org/conda-forge/diffeqpy https://anaconda.org/conda-forge/diffeqpy. Julia has language bindings that make it simple enough to bind to other high level languages that these projects are maintained by the core team and supports lots of the library, including forms of automatic differentiation and GPU kernel generation. See for example the bindings that allow for usage within PyTorch https://github.com/SciML/juliatorch https://github.com/SciML/juliatorch and the Python-based Collab notebooks showing the GPU usage https://colab.research.google.com/drive/1bnQMdNvg0AL-LyPcXBiH10jBij5QUmtY?usp=sharing https://colab.research.google.com/drive/1bnQMdNvg0AL-LyPcXBi.... With Julia v1.12's small binary generation, we plan to release forms via binaries with C ABIs over the next year as well. Are those sufficient or should we consider supporting other deployments?