Incorrect computation of Hessian of squared euclidean norm

[jomi2801]

The computation of the Hessian of the squared euclidean norm results in an interesting output:

using ForwardDiff;
f(x::Vector) = .5*norm(x, 2)^2;
x = zeros(2);
ForwardDiff.hessian(f, x)

which is

 0.0  0.0
 0.0  1.0

However, the codes

f(x::Vector) = .5*x[1]^2 + .5*x[2]^2;
ForwardDiff.hessian(f, x)

and

f(x::Vector) = .5*dot(x, x);
ForwardDiff.hessian(f, x)

provide the correct result. Why is that?

[jrevels]

Possibly similar to #197 i.e. Base's linear algebra code might be optimized in an "AD unsafe" way.

[JoshuaLampert]

Is there any chance this gets fixed? Unfortunately, in my case it is not really an option to replace the norm(x)^2 by sum(x.^2) since I have a function calling norm(x) and later another function squaring the result and because the concatenation gives NaN:

julia> using ForwardDiff

julia> f(x) = 0.5*sqrt(sum(x.^2))^2
f (generic function with 1 method)

julia> ForwardDiff.hessian(f, [0.0, 0.0])
2×2 Matrix{Float64}:
 NaN  NaN
 NaN  NaN

Enabling nansafe_mode also doesn't give the correct result:

julia> ForwardDiff.hessian(f, [0.0, 0.0])
2×2 Matrix{Float64}:
 0.0  0.0
 0.0  0.0

So this would require to rewrite my code base, which I would really like to avoid.

[JoshuaLampert]

This issue happened to me again when trying to upgrade to ForwardDiff.jl v1. Even worse: Upgrading to ForwardDiff.jl@1, this also concerns gradient.
On ForwardDiff.jl v1.0.1:

julia> using ForwardDiff, LinearAlgebra

julia> f(x) = 0.5 * norm(x)^2
f (generic function with 1 method)

julia> ForwardDiff.hessian(f, [0.0, 0.0])
2×2 Matrix{Float64}:
 NaN  NaN
 NaN  NaN

julia> ForwardDiff.gradient(f, [0.0, 0.0])
2-element Vector{Float64}:
 NaN
 NaN

On ForwardDiff.jl v0.10.38:

julia> using ForwardDiff, using LinearAlgebra

julia> f(x) = 0.5 * norm(x)^2
f (generic function with 1 method)

julia> ForwardDiff.hessian(f, [0.0, 0.0])
2×2 Matrix{Float64}:
 0.0  0.0
 0.0  1.0

julia> ForwardDiff.gradient(f, [0.0, 0.0])
2-element Vector{Float64}:
 0.0
 0.0

Is there really no way to fix this or is there at least a workaround? Again, replacing norm(x)^2 by something like sum(abs2, x) or sum(x.^2) is not an option for me unfortunately because my code is a little bit more complex where norm and ^2 appear in different locations of the code.

[andreasnoack]

I think the problem here is the mathematical function, not the code. Derivates of the euclidian norm at zero are just ambiguous. The limit depends on how you get to zero. I.e. with

$$f(x,y) = \sqrt{x^2 + y^2}$$

then

$$\lim_{t \overset{+}{\rightarrow} 0} \left. \frac{\partial f}{\partial x} \right|_{x=t,y=0} = 1$$ $$\lim_{t \rightarrow 0} \left. \frac{\partial f}{\partial x} \right|_{x=0,y=t} = 0$$ $$\lim_{t \overset{+}{\rightarrow} 0} \left. \frac{\partial f}{\partial x} \right|_{x=t,y=t} = \frac{1}{\sqrt{2}}$$

See also https://math.stackexchange.com/questions/310325/is-the-euclidean-norm-differentiable-at-0. Hence, I think returning NaN makes sense. In life, it is just better to avoid roots and power functions if you want to compute derivatives. If you can't, then you'd have to stay away from zero.

[JoshuaLampert]

I am aware that the norm is not differentiable in zero, but I am not interested in the derivative of the norm in zero. The thing is that after applying the square, the function is differentiable everywhere. I also know that ForwardDiff.jl uses the chain rule to apply derivatives and therefore also stumbles upon the derivative of the norm, which leads to the issue. However, I am surprised that this seemingly general issue of forward mode AD is not discussed more (or is it?) since IMHO it's totally common that differentiable functions are concatenations of non-differentiable functions. If forward mode AD is not able to differentiate these functions, I consider it as a major drawback. So I was wondering how others tackle this problem and if there are workarounds. Of course one workaround would be to perturb zero by some eps, but this doesn't seem optimal.

[andreasnoack]

I believe that the solution is to avoid the composition with a non-differentiable intermediate solution, i.e. to use something like a norm² function if sum(abs2, ...) for some reason isn't working for your application. To make norm(x)^2 avoid the problem here would require some compiler magic that I don't think we'll ever get. For the compiler to make the square and the square root cancel, it would have to figure out that the argument to the square root is non-negative.

I my application domain, the power function often appears in ecpressions, we'd like to compute derivates for with this package. The only way to make it work seems to be to add a tiny amount to avoid hitting zeros (or go negative). If anybody has better solutions, I'd also love to hear more.

[JoshuaLampert]

Thanks for your response!
Let me briefly explain my use case to clarify my situation. Maybe this also helps someone to come up with a solution. The problem is that in the package I'm developing (KernelInterpolation.jl if you are interested), I have a general function $\Phi$, which gets a vector $x$ and passes the norm of this vector $r = \Vert x \Vert$ to another function $\varphi$. The definition of $\varphi: \mathbb{R}^+\to\mathbb{R}$ depends on the specific radial basis function (and there are potentially many of them, which could also be user-defined). Some of them square the input (which is the norm of $x$) others apply some other combination to $r$ (e.g. $$(1 + r)\exp(-r)$$), but all of them are differentiable. I do not want to pass $x$ to these functions $\varphi$ because they only depend on the norm of $x$ and I want this to be reflected in the code (and users who want to extend this function should also only write a function, which obtains $r$ and not $x$). Actually, thinking about it, I would also not be 100% how to even rewrite (the differntiable!) $f(x) = (1 + \Vert x\Vert)\exp(-\Vert x\Vert)$ without using the norm since we never apply a square (or any power) to the norm of $x$. It is also quite common to evaluate gradients of $\Phi$ in zero in the context of RBFs.

[devmotion]

which gets a vector x and passes the norm of this vector r = ‖x‖ to another function φ.

Couldn't you pass the squared norm (ie. r^2) to φ?

In KernelFunctions.jl (which uses a similar separation of functions it seems), we avoid taking the square root in the first place (as suggested by @andreasnoack) in such cases by passing the squared Euclidean norm when working with the Euclidean norm, e.g., in https://github.com/JuliaGaussianProcesses/KernelFunctions.jl/blob/60cc391cf24a65f6de48776153bbfe4ff20c2dcb/src/basekernels/exponential.jl#L23-L27.

[JoshuaLampert]

Thanks for the suggestion. This would solve the problem for all kernels depending on the squared norm. However, as I noted above that is not the case for all of the kernels. So I cannot always pass the squared norm. Let's again take a look at the Matern32Kernel with Euclidean distance, for which we have (up to constants) $$\varphi(r) = (1 + r)\exp(-r)$$. It seems like also KernelFunctions.jl cannot use ForwardDiff.jl to compute the gradient of $$x\mapsto\varphi(\Vert x\Vert)$$:

julia> using KernelFunctions, ForwardDiff

julia> k = Matern32Kernel()
Matern 3/2 Kernel (metric = Distances.Euclidean(0.0))

julia> f(x) = k(x, zero(x))
f (generic function with 1 method)

julia> ForwardDiff.gradient(f, [0.0, 0.0])
2-element Vector{Float64}:
 NaN
 NaN

[JoshuaLampert]

Oh, I just saw this: JuliaGaussianProcesses/KernelFunctions.jl#571.

[andreasnoack]

Not sure how things are implemented in KernelFunctions but if you defined $\phi$ as a function of the squared norm $$\phi(r) = (1 + \sqrt{r})\exp(-\sqrt{r})$$ and provided the derivate rule $$\phi'(r) = \frac{-1}{2}\exp(-\sqrt{r})$$ then it should work but it seems difficult to make it work if you want to use ForwardDiff for the derivative.

[JoshuaLampert]

Yeah, I would have hoped not needing to define my own derivatives.

[devmotion]

I don't think you can avoid defining an AD rule in this case.

it seems difficult to make it work if you want to use ForwardDiff for the derivative.

In this case, it should be sufficient to define something like

function matern32_phi(r2::Real)
    r = sqrt(r2)
    return (1 + r) * exp(-r)
end

function matern32_phi(r2::ForwardDiff.Dual{T}) where {T}
    r = sqrt(ForwardDiff.value(r2))
    z = exp(-r)
    return ForwardDiff.Dual{T}((1 + r) * z, (-z/2) * ForwardDiff.partials(r2))
end

[JoshuaLampert]

Ok, I see. Thanks for the clarification. Would it make sense to document this under "Limitations of ForwardDiff"? If yes, I can create a PR tomorrow.

[andreasnoack]

I think it would be good to document this and explain how you can sometimes work around it by using helper functions.