Using ForwardDiff results in a stack overflow

General Usage
1

Hi, im using the ForwardDiff package. The following code is only a simplification of the real problem.
But I am getting there as well a stack overflow in the line of the LU decomposition.

using ForwardDiff, LinearAlgebra, SparseArrays


function foo(x)
  r = [1,2,1];
  c = [1,2,2];
  v = [x[1],x[2],x[1]];
  rhs = [1.0,2.0];
  A = sparse(r,c,v);
  LU = lu(A)
  y = LU \ rhs;
  return y
end
J = ForwardDiff.jacobian(foo,[1.,2.])
2

Interesting that this pops up as a stack overflow error…

The point is that \ for SparseMatrixCSC is channeled to UMFPACK from SuiteSparse which is implemented in C and thus not able to work with ForwardDiff.Dual numbers.

That said, recently, the package Sparspak.jl was announced which is a pure Julia implementation of the classical Sparspak by George and Liu and thus provides an LU factorization working with general number types. The ExtendableSparse.jl package (shameless plug: I am the author…) provides the ExtendableSparseMatrix type which has a \ which falls back to Sparspak on number types not supported by SuiteSparse/UMFPACK .

So now we have:

using ExtendableSparse, ForwardDiff, LinearAlgebra, SparseArrays
function foo(x)
         r = [1,2,1];
         c = [1,2,2];
         v = [x[1],x[2],x[1]];
         rhs = [1.0,2.0];
         A = ExtendableSparseMatrix(sparse(r,c,v));
         LU = lu(A)
         y = LU \ rhs;
         return y
end

J = ForwardDiff.jacobian(foo,[1.,2.])
2Ă—2 Matrix{Float64}:
 -1.0   0.5
  0.0  -0.5
1 Like
3

Hi,
thank you for your quick reply. Somehow I am still getting an Error with your code (I copied your version).

ERROR: MethodError: Cannot convert an object of type
Sparspak.SpkProblem.Problem{Int64{},Float64} to an object of type
Sparspak.SpkProblem.Problem{Int64{},ForwardDiff.Dual{ForwardDiff.Tag{typeof(foo), Float64}, Float64, 2}}

4

Strange, I tried this even on Julia 1.6 in a clean environment, and it all works.

The error message indicates that you somehow tried to differentiate through a sparse solve with a sparse matrix already created with Float64 entries instead of ForwardDiff.Dual.

5

After i restarted my enviroment everything worked fine.
Thanks for your professional help :slight_smile: