Differentiation matrix stack overflow

Specific Domains Numerics
1

I’m trying to follow along with @stevengj’s lecture notes on 2D differentiation matrices, but I can’t get past the first example without Julia reporting a stack-overflow error. The piece of code in question is

# construct the (M+1)xM matrix D, not including the 1/dx factor
diff1(M) = [ [1.0 zeros(1,M-1)]; diagm(ones(M-1),1) - eye(M) ]
diff1(5)

which on my machine results in

ERROR: StackOverflowError:
Stacktrace:
 [1] diagm(::Vector{Float64}, ::Int64) (repeats 79984 times)
   @ ChainRulesCore ~/.julia/packages/ChainRulesCore/8vlYQ/src/tangent_types/thunks.jl:66

Any ideas what could be the problem here? I’m using the 1.6.3 generic binary for 64-bit Linux.

2

Those notes are from quite an old version of Julia (0.4) — that’s not how the diagm function is called anymore. An updated version of that function would be simply:

using LinearAlgebra
diff1(M) = let o = ones(M); diagm(M+1, M, 0=>o, -1=>-o); end

(And here is a sparse version using spdiagm.)

2 Likes
3

A separate problem is that ChainRulesCore is engaging in “type piracy” of the diagm function here, which is why you are getting a mysterious error message (instead of a MethodError, which would be more sensible) here. This should be fixed in ChainRulesCore.jl: https://github.com/JuliaDiff/ChainRulesCore.jl/issues/472

2 Likes
4

Thanks for the help, and for making these notes available!