Hacker News Reader

Comment by yarg | Hacker News Reader

One thing that I've thought about, when using calculus in programs, you're often dealing with a very small (but finite) Δx, rather than an infinitesimal 𝛿x.

And I don't think that the common differential equation is the ideal form when dealing with that situation.

(Take for example g(x) = -f(-x), the gradients at x and -x are not equivalent for Δx, but would be for 𝛿x.)

Anyway,

    (limit(h -> 0)((f(x + h) - f(x))/h)

Why do we even need h? It's just (the infinitesimal) 𝛿x multiplied by an arbitrary finite constant.

(You're dealing with a linear equivalent infinitesimal subsection of the graph - so 𝛿x and 𝛿y scale linearly with each other.)

So remove h:

    ((f(x + 𝛿x) - f(x))/𝛿x)

Looks cleaner.

    (f(x + 𝛿x) - f(x))/𝛿x ≡ (f(x) - f(x - 𝛿x))/𝛿x ≡ (f(x + 𝛿x) - f(x - 𝛿x))/(2 * 𝛿x)

You can take the gradient from x to (x + 𝛿x) or from (x - 𝛿x) to x, or from (x - 𝛿x) to (x + 𝛿x).

(The cube-root of the other three forms is also a valid differential equation, if you want to be evil about it.)

    (f(x + 𝛿x) - f(x - 𝛿x))/(2 * 𝛿x)

(That's a nicer form for programmatic use; resolving the earlier mentioned gradient issue when dealing with finite deltas.)

Another thing I noticed: the standard (h -> 0) form eliminates all parts of the gradient containing 𝛿x values - which is fine for infinitesimal 𝛿x, but is less ideal for finite Δx.

View on HN

You really just derived the central derivative. Yes, this derivative is more accurate. The issue is that you can't use it on the edges. For the boundaries you're stuck with lower precision.

Here's a resource you might be interested in http://www2.math.umd.edu/~dlevy/classes/amsc466/lecture-note...

The central derivative is equal to Im(f(x + j h) / h) where j is the split-complex number which squares to 1, "Im" is the imaginary-part function, and h is a small real number. When j is replaced with the dual number epsilon, this instead equals the actual derivative, as is shown in the article. You can also replace j with the complex number i, producing the complex-step method.

Central difference and complex-step are both bad approximations to the actual derivative. Of these, the complex-step method is the most ridiculous because it requires extending a function to act on complex numbers, while not extending it to the closely-related dual numbers.

"Ridiculous" is not a word I'd use... sometimes symbolic derivatives are intractable, in which case complex step derivatives provide high-precision numerical derivatives which are otherwise unattainable.

https://vladfeinberg.com/2020/07/05/discrete-residues.html

No. The dual numbers are very similar to the complex numbers and produce exact results. Ridiculous is true and fair.

[edit] Also, I wouldn't count using Cauchy's Integral Formula as the complex-step method.

This!

I always find it curious how computer scientists end up rediscovering things known for decades, and then get all hyped about it ...

I figure it's because math people aren't great at communication with people outside their domain. Point in case being Wikipedia articles on math being absolute gobledygook to non-math people.

To be fair, I can't get computer scientists to read any math books. But I don't think this is quite odd. The only people I know that go out and read math books are mostly (current or former) mathematicians or physicists (which a random sprinkling of others). I'm not sure it is exactly mathematicians faults that your third grade teacher impressed upon you that math is hard and useless.

I don't think math is hard or useless, thank you very much. However if I'm programming a game and want to project a point onto a line I'd prefer to not get told by math folks that I first need to learn all of linear algebra before they can get their solution explained to me. I'm not stupid, I understand plenty of math, I'm just not steeped in the jargon.

I feel like you're responding as if I attacked you. I'm not sure why. Maybe it was the comment about not being able to get CS people to read math books but you do? I don't mean that statement in an absolute way (there's no absolutes) but rather as a general case. Of course there are CS people that read math books. It is just my experience that the majority aren't interested in it and I don't think this is an uncommon experience. If you're different, good for you, but don't take general statements personally. I'm not calling you stupid nor am I saying you don't have an interest in math. You're an individual person.

"I'm not sure it is exactly mathematicians faults that your third grade teacher impressed upon you that math is hard and useless." doesn't exactly strike me as a general statement, but perhaps I overreacted.

If you're interested in how to best calculate derivatives numerically using a finite difference rather than an infinitesimal, try having a look at the Wikipedia page on numerical differentiation. The punchline is that you can control the error in f'(x) much better if you use many displacements, ie. f'(x) = a1 f(x + b1 h) + a2 f(x + b2 h) + a3 f(x + b3 h)..., where the bi are chosen and the ai depend on the bi and h.

Such methods are doomed. Consider the function

  f(x) = a sin(x / a)

for very small a. Finite differencing methods would approximate f'(0) as 0 as f(x) is very close to 0, but the actual value of f'(0) is 1.

Differentiation is a discontinuous operator, and cannot be accurately estimated using numerical methods. Only exact methods, which include the dual numbers, can produce accurate results.

This is not true. Doomed is sensational. With proper perturbation analysis, one can choose an appropriate stencil and formula that has desired error bounds.

https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.72...

https://www.google.com/search?q=finite+difference+stencil+de...

Won't work on my example. Given any differentiable function f:R -> R, you can add a very small function to it that dramatically changes its derivative. But maybe this is too pathological to bother about.

Over the complex numbers, this phenomenon doesn't happen. This can be shown, for instance, using the maximum modulus principle. You can therefore use Cauchy's Integral Formula to accurately estimate derivatives. But why not then just use the dual numbers?

I see your edit - fair point, analyticity does make complex functions much more well-behaved than their real counterparts. Still, this isn't a death knell like you previously claimed.

Let's put up a bet. I win the bet if I provide a polynomial approximation of f'(x), R[f(x)], bounded by a maximum of 5% error, on the interval (-1, 1), by November 1st, 2022. If I fail to, or resign by then, you win.

If I win, you agree to silence yourself for 1 year, regarding the split-complex and dual numbers, no longer touting them as a panacea.

If I lose, I'll agree to make a YouTube video extolling the virtues of a mathematical subject of your choice, such as these types of numbers.

I don't fully understand your bet (and I don't take kindly to being asked to silence myself). I don't have time at the moment to craft an example. I might play around a bit with the edges of floating point. Differentiation is indeed discontinuous, so with enough determination I should be able to win such a bet. My function would probably be a sigmoid function which would go from -a to +a (for small a) over a very small interval [-b,b], where b can equal a.

> If I win, you agree to silence yourself for 1 year, regarding the split-complex and dual numbers, no longer touting them as a panacea.

I never said these numbers were a panacea. I said the dual numbers were an exact method. In this comment [1], which is what you're referring to, I said that the connection with other planar algebras was cute, but you mischaracterised what I said, which was rude. The algebras are pertinent because the article mentioned the dual numbers.

[1] - https://news.ycombinator.com/item?id=32305644

Fair enough. The bet I had in mind was the function you'd already stated as intractable: f(x) = a sin(x / a), the value of a is up to you, within reason (no IEEE fp shenanigans.)

If you'd rather some other consequence than silencing, please suggest one. I am tired of hearing your panaceas. Symbolic differentiation is also exact but intractible. What makes the duals yield exactness in any way that isn't prone to the same finitude as other approximations?

I see a lot of beasts but only one individual claiming to tame them with ease and exactness.

If it comes across as rude, that's not my intent. Moreso to have some clarity and to have your claims be validated.

I'm not going to accept your bet because I don't understand its conditions (and what's the point?), but I accept your challenge. Also, I'm going to f*** around bigly with floating point.

  def f(x,a):
    return 1 + a * sin(x / a)

Now let epsilon = 2.2204460492503133e-36. Consider f(x, epsilon).

The function should always return 1 in floating point, so any finite differencing method will estimate its derivative as 0. But over the dual numbers, the function should have f(e, epsilon) = e, where e is the dual number imaginary. In other words, the dual numbers return the correct value of the derivative at x=0, which is 1. You should be able to use any Python implementation of the dual numbers; my own one uses Sympy. Actually, you might be able to use Scipy and represent dual numbers as matrices.

Look up automatic differentiation because you're not aware of it, and it is the subject of the article. You seem to misunderstand symbol methods. The problem of finding an exact derivative at a single point is not as hard as finding the symbolic derivative everywhere. You're referring to the product rule and chain rule, but these are not needed to find derivatives at single points.

I am aware of AD. All AD does is use the chain rule on the Taylor (or otherwise) approximations of all arithmetic functions used in the calculation of the target function. It's symbolic but only at the single arithmetic operation level. The chain rule is concealed under rules on epsilon. The reason epsilon's square equals 0 is to eliminate any higher order derivatives as soon as they appear. If the derivative is abnormal in a way such that the normal truncated calculative formulas would conceal it, then AD would fail to give an accurate result as well.

You've just described the dual numbers, which provide a way of implementing forward-mode autodiff.

  > If the derivative is abnormal in a way such that the normal
  > truncated calculative formulas would conceal it, then AD
  > would fail to give an accurate result as well.

What I'm seeing directly contradicts that. I've tested the example above (the one I called f(x,a) with a=1e-36, trying to find its derivative at x=0), and it gave the right answer with dual numbers, but not with finite differencing. So what you're saying isn't true. I'll post a minimal implementation here later.

[edit]

I've used Sympy out of laziness, which is inelegant. I don't know how clear this will be.

First of all, I have to change the number 1 in the definition of f(x,a) into a matrix. We need this because we'll be representing dual numbers as matrices, and you can't add the scalar 1 to a matrix:

  one = eye(2)
  
  def f(x, a=1e-36):
    return one + a * sin(x / a)

Sympy doesn't have its own implementation of `sin`, so we'll need to provide one:

  def sin(M):
    return im(exp(I*M).n())

At the dual number ε, we have f(ε) = f(0) + ε f'(0). The dual number ε will be represented as the matrix

  ⎡0  1⎤
  ⎢    ⎥
  ⎣0  0⎦

As code:

  epsilon = Matrix([[0,1],[0,0]])

Evaluating f on the above matrix gives:

  In: f(epsilon)
  Out: 
  ⎡1  1.0⎤
  ⎢      ⎥
  ⎣0   1 ⎦

Which tells us that f(0) ≈ 1 and f'(0) ≈ 1. Both are correct.

If we use finite differencing instead, we get:

  In: (f(h * eye(2)) - f(0 * eye(2)))/h
  Out: 
  ⎡0  0⎤
  ⎢    ⎥
  ⎣0  0⎦

So it's claiming incorrectly that the derivative at 0 is 0. No finite differencing scheme can fix this because floating point causes the function f to become constant.

Autodiff using the dual numbers has pulled off the seemingly impossible.

[edit: An early typo caused the code to produce an incorrect result. Now fixed, and my claim holds.]

> Autodiff using the dual numbers has pulled off the seemingly impossible. It's not seemingly impossible, if you understand that the chain rule for f' is just being executed at the same time as the calculation of f by having derivatives for basic operations already defined.

However, like I said, if you have a calcuation method that hides derivatives in terms that have been truncated, then this will not save you.

(-1)*i * x*(2*i + 1) / factorial(i) is not one of those methods -- and sin is rather regular in regard to its Taylor series. So of course it works out in this case, with a normal power series calculation method.

Try a different series or calculation method, and dual numbers will get you a wildly different result. Understand, dual numbers only work well, when you use a method of calculation that front-loads terms that have high bearing on the derivative. Otherwise the missing terms/truncation, causes severe inaccuracy.

However, stenciling might actually perform better in these scenarios.

> Try a different series or calculation method, and dual numbers will get you a wildly different result

No. The example works because:

  While the exact value of f(x,a) isn't 1, given any inexact representation of real numbers like floating point or fixed point, the value of "a" can be chosen so that f(x,a) has 1 as its closest representation.

Trying to compute f(x,a) differently isn't going to change that, so stencilling methods are never going to work here. But autodiff will always work. This means I win your challenge.

Your other claims are probably gibberish. You need to provide an example.

What don't you understand -- the method of calculating the function is an approximation, thus the AD derivative is dependent on it, and the AD derivative is a derivative of the approximation, not the actual function. Whereas an approximation of the actual derivative is what we are truly after.

> The key take away here is that the map is not the territory. Most nontrivial functions on computers are implemented as some function that that approximates (the map) the mathematical ideal (the territory). Automatic differentiation gives back a completely accurate derivative of the that function (the map) doing the approximation. Furthermore, the accurate derivative of an approximation to the idea (e.g d_my_sin), is less accurate than and approximation to the (ideal) derivative of the ideal (e.g. my_cos). There is no truncation error in the work the AD did; but there is a truncation error in the sense that we are now using a more truncated approximation that we would write ourselves.

https://www.oxinabox.net/2021/02/08/AD-truncation-error.html

AD is great but if you have a calculation method ill-suited for AD, then you'll get shite results. Why is this surprising?

And yeah, stenciling is mostly for PDEs and other state spaces that we don't have a closed form for. It's generally not used for an analytic function. But you can use it for an analytic function if you tailor the stencil to the function. In fact, you'll just yield a truncated Taylor polynomial if you provide a perfect function-specific stencil.

> If the derivative is abnormal in a way such that the normal truncated calculative formulas would conceal it, then AD would fail to give an accurate result as well.

The following is proof of otherwise. First implement the dual numbers by borrowing the code from here [1]. Then define the sine function in a crude way:

  from math import factorial
  
  def sin(x):
      return sum((-1)**i * x**(2*i + 1) / factorial(i) for i in range(100))

We then define f:

  def f(x, a=1e-36):
      return 1 + a * sin(x / a)

and use the dual numbers to estimate the derivative at 0:

  In [49]: print(f(Dual(0,{'epsilon':1})))
  f = 1.0
  fepsilon = 1.0

We get 1 as the estimate, which is correct! For comparison, we may try using finite differencing to estimate the derivative at 0:

  In [50]: h = 1e-36; (f(h) - f(0))/h
  Out[50]: 0.0

The estimate for the derivative is 0 here, which is wrong. The finite differencing cannot be improved using stencil theory because the estimate of f(x) is always exactly 1 under floating point, and therefore the estimate for f'(x) is always exactly 0. Dual numbers win.

[1] - https://github.com/ujjwalkhandelwal/Dual-numbers-and-automat...

Finite difference derivatives are the basis for a range of numerical methods for solving differential equations. Given that these techniques underpin various areas of modern engineering, I'd say your statement that 'differentiation [...] cannot be accurately estimated using numerical methods' is demonstrably incorrect. I think you're fishing for something like 'there exist functions whose derivatives cannot be accurately estimated numerically'. Presumably some assumptions have been made which rule out such functions in using finite difference schemes in practical applications.

> Finite difference derivatives are the basis for a range of numerical methods for solving differential equations.

I thought the basis usually involved rewriting the equation in terms of integrals, which can be estimated using quadrature. Numerical quadrature does not have the same problems as finite differencing because integration is a continuous operator. How accurate are these methods you're proposing? Could they be improved using autodiff (via the dual numbers or via reverse-mode autodiff)? Do you have a reference, so I can understand your point better? I'm sceptical and I'm not inclined to believe you yet.

Well, there are many different methods and certainly some are based on numerical integration as you say. The most simple method based on finite difference derivatives is Euler's method which uses the scheme given by OP, with error O(h). This can be extended to Runge-Kutta methods which use a finite difference with n terms and have error O(h^(n-1)), and which to my knowledge are common and widely used to solve ODEs. You can read about Runge-Kutta methods on Wikipedia, and they also have a page on numerical methods for ODEs which gives an overview of different methods.

I don't know what autodiff is, or how the dual numbers fit in here. I'd be interested to learn if you want to explain

> I don't know what autodiff is, or how the dual numbers fit in here. I'd be interested to learn if you want to explain

Autodiff is an algorithm for computing derivatives in areas like Deep Learning. It's an exact method with two main variants -- forward mode and reverse mode -- differing in terms of when each of them is more efficient. It's implemented in widely-used numerical libraries like PyTorch and TensorFlow. The simplest implementations use the dual numbers to implement the forward mode variant. More sophisticated implementations implement the reverse mode variant. Wikipedia has plenty on it.

Ok I had a look at wiki and it's clear to me how the dual numbers encode differentiation and the theory behind how autodiff works. How do the two fit together? I didn't follow that part

You define a dual number type in your programming language of choice, and overload the operators +, -, *, /, <, >, sin, cos, etc. Then you take a function which computes a numerical function f which is defined using the previous operations, and simply evaluate them on a dual number. It's "automatic" because, assuming the right programming language, a user-defined numerical function will often be able to take in many different kinds of numbers, like complex numbers, real numbers, arbitrary precision floats, and even user-defined classes like the dual numbers.

Ah thankyou, it was so neat and automatic that I missed the punchline. It looks like you'll always calculate f and f' together and it's not possible to calculate only f' though? Seems like this could waste time in some applications?

I can see how this is the solution to the problem of calculating the numerical derivative of a known function composed of finite combinations of rational functions and a set of specified functions, elementary or special. I don't see how it's relevant to solving ODEs numerically though, given that in that case you don't know the functional form of f or f'

> I can see how this is the solution to the problem of calculating the numerical derivative of a known function composed of finite combinations of rational functions and a set of specified functions, elementary or special.

Yes and no. It's worth taking into account that you can apply it to some functions defined using conditions and loops. For instance, if f is the sqrt function defined using Newton's method, then autodiff via the dual numbers will compute an estimate for f'(x) = 1/(2 * sqrt(x)). More explicitly:

  def f(a):
    x = 1
    oldx = 2
    threshold = 1e-10
    while abs(x - oldx) < threshold:
      oldx, x = x, (x + a/x) / 2
    return x

^ You can evaluate f on a dual number. This algorithm contains conditions and loops.

This was an interesting comment and I learnt something new, thank you. It does feel a little bit like you are proselytising the good news of the dual number autodiff method though, when you don't respond to questions on potential downsides and only add information about positives.

I'm more familiar with using differentiation in optimisation, where autodiff seems like the best method. I'm also aware of the fact that the differentiation operator is discontinuous, which makes me sceptical of finite-differencing methods. I'm not familiar with numerically solving differential equations, but the natural connection between Runge-Kutta and finite differencing seems to be an interesting one.

An entire subfield of analysis called Stenciling exists, for this purpose. Depending on the function or differential equation / system, different stencils are used.

Stenciling doesn't just deal with the derivative but it tries to come up with approximations involving a fixed number of sample points, the stencil, for any differential operator, ie. the Laplacian, higher order derivatives, etc.

> The term "stencil" was coined for such patterns to reflect the concept of laying out a stencil in the usual sense over a computational grid to reveal just the numbers needed at a particular step. [2]

[0] https://en.wikipedia.org/wiki/Five-point_stencil

[1] https://en.m.wikipedia.org/wiki/Finite_difference_coefficien...

[2] https://en.wikipedia.org/wiki/Stencil_(numerical_analysis)

Thanks for the links (and the flame war).

Another thing that I'm wondering about is how to calculate the length of a function.

What I'm thinking, is that you can treat a function like a piece of string, you take f(x) and generate a pair of functions, one above and one below f(x).

You generate the functions at a fixed (infinitesimal) tangential distance 𝛿t and generate deltas for {x, y}

    mul = 𝛿t/((𝛿x)^2 + (𝛿y)^2)^.5
    Δx = 𝛿y * mul
    Δy = 𝛿x * mul

    u(x - Δx) = f(x) + Δy
    l(x + Δx) = f(x) - Δy

The idea being that the function length is the area ∫(u(x) - l(x)) divided by the thickness 2 * 𝛿t.

I have no idea how wrong this is, but if you can point me toward any resources that'll help when I (inevitably) get stuck, it would be very much appreciated.

Not exactly sure what you are getting at but it does remind me of the Cauchy's Integral Formula which shows that differentiation is equivalent to integration.

https://en.wikipedia.org/wiki/Cauchy%27s_integral_formula

Trying to measure the length of a curve f(x).

I turn the curve into a fixed infinitesimal thickness (2 * 𝛿t) 'rope' and measure the area.

Then I divide the area by the thickness to get the length.

The rope is defined as the area between two curves, both of them manipulations of f(x).

The curves are generated by moving each point (in f(x)) a constant infinitesimal distance at a tangent to the function (this gives the rope its fixed thickness), one in either direction.

For the upper curve {Δx, Δy} = {(n * 𝛿y), -(n * 𝛿x)}, for the lower {(n * 𝛿y), -(n * 𝛿x)}, where n is a normalisation constant 𝛿t/((𝛿x)^2 + (𝛿y)^2)^.5

Once you have these two functions, you can use the integral of the difference to get the area.

That's what I'm thinking anyway.

I'm a brainlet or maybe it's the 10+ hours of work I did today but I just can't grasp it right now.

There are a lot of ways to take integrals though, so I wouldn't doubt your method could work. Generally though the holes in a lot of these integration methods is their inability to work with functions that are pathological: https://en.wikipedia.org/wiki/Pathological_(mathematics)

|group1 = Types of integrals |list1 = * [[Riemann integral]] * [[Lebesgue integration|Lebesgue integral]] * [[Burkill integral]] * [[Bochner integral]] * [[Daniell integral]] * [[Darboux integral]] * [[Henstock–Kurzweil integral]] * [[Haar measure|Haar integral]] * [[Hellinger integral]] * [[Khinchin integral]] * [[Kolmogorov integral]] * [[Lebesgue–Stieltjes integration|Lebesgue–Stieltjes integral]] * [[Pettis integral]] * [[Pfeffer integral]] * [[Riemann–Stieltjes integral]] * [[Regulated integral]]

|group2 = Integration techniques |list2 = * [[Integration by substitution|Substitution]] * [[Trigonometric substitution|Trigonometric]] * [[Euler substitution|Euler]] * [[Weierstrass substitution|Weierstrass]] * [[Integration by parts|By parts]] * [[Integration by partial fractions|Partial fractions]] * [[Integration using Euler's formula|Euler's formula]] * [[Integral of inverse functions|Inverse functions]] * [[Order of integration (calculus)|Changing order]] * [[Integration by reduction formulae|Reduction formulas]] * [[Integration using parametric derivatives|Parametric derivatives]] * [[Leibniz integral rule#Evaluating definite integrals|Differentiation under the integral sign]] * [[Laplace transform#Evaluating improper integrals|Laplace transform]] * [[Contour integration]] * [[Laplace's method]] * [[Numerical integration]] * [[Simpson's rule]] * [[Trapezoidal rule]] * [[Risch algorithm]]

(No idea if you'll ever see this, I also have no damned idea what I'm talking about.)

Looking at the Weierstrass function, it's a chaotically convergent series with a non-convergent gradient series.

(The standard differential forms don't apply, because the assumption of linearisation at the scale of 𝛿x is invalid, due to the scale invariant properties of the fractal.)

In that case the integral function would be even more convergent since the 'noise' (pretty much) adds to zero.

So can you differentiate via integration, and get a gradient for a related function (with less noise)?

I messed around with this, taking the areas of a pair of little triangle approximations and adding them to get a ~rectangle.

Then divide that area by 𝛿x to get 𝛿y, and divide that by 𝛿x to get 𝛿y/𝛿x.

This was the form that I reached:

    let y = f(x)
    let F(x) = ∫f(x)𝛿x
    𝛿y/𝛿x = (F(x + 𝛿x) + F(x - 𝛿x) - 2 * F(x))/((𝛿x) ^ 2)

At very least, it passed my polynomial sanity check:

    let F(x) = x^6
    𝛿y/𝛿x = 30(x^4) + 30(x^2)(𝛿x^2) + 2(𝛿x^4)

That's awesome. Glad to see your method worked. If you ever have a chance to encounter a serious mathematician, they can probably give you the proper name of it. Or perhaps we'll have to name it the Yarg integral :-)

I've just been thinking of it as differentiation via integration.

It took me until the end to realise that all I'd done was derive a form for the second derivative of F(x).

Ah yes, this is why epsilon squared equals 0 in the dual numbers -- to automatically remove all higher order derivatives when they appear. Exponents on derivatives are indicate of their rank/order.

It's the second symmetric derivative.

https://en.wikipedia.org/wiki/Symmetric_derivative#The_secon...

https://en.wikipedia.org/wiki/Second_derivative#Limit

> The limit is called the second symmetric derivative. Note that the second symmetric derivative may exist even when the (usual) second derivative does not.

> This limit can be viewed as a continuous version of the second difference for sequences.

So it actually can do what I intended it to do, nothing original - but there was a degree of satisfaction in deriving it myself.

Thanks yarg, I just learned something new.

Is h use instead of 𝛿 or Δ, because it is applicable to both partial and normal differential. 𝛿 Is always a partial symbol to me.

d in dx whilst is not the same as h, as it stands its own way now as “operator”.

No, the h here is a number, not an operator.

I just find that it buries the lede.

Reply on news.ycombinator.com