I’m trying to pre-allocate a cache to be used with SparseDiffTools and PolyesterForwardDiff, but I’m struggling to get a working example.
My attempt is based off of the jacobian tests:
difftype = ADTypes.AutoSparse(ADTypes.AutoPolyesterForwardDiff(; chunksize = 1))
J_sparsity = ForwardDiff.jacobian!(dx, obj, y, x0)
# row_colorvec = spd.matrix_colors(J_sparsity; partition_by_rows = true)
# col_colorvec = spd.matrix_colors(J_sparsity; partition_by_rows = false)
# sd = spd.PrecomputedJacobianColorvec(; jac_prototype = J_sparsity, row_colorvec, col_colorvec)
colors = spd.matrix_colors(J_sparsity)
sd = spd.PrecomputedJacobianColorvec(; jac_prototype = J_sparsity, colorvec=colors)
cache = spd.sparse_jacobian_cache(difftype, sd, obj, x0)
The final line there fails. The error message is:
ERROR: MethodError: no method matching findstructralnz(::Matrix{Float64})
The function `findstructralnz` exists, but no method is defined for this combination of argument types.
Closest candidates are:
findstructralnz(::SparseArrays.SparseMatrixCSC)
@ ArrayInterfaceSparseArraysExt ~/.julia/packages/ArrayInterface/9um5f/ext/ArrayInterfaceSparseArraysExt.jl:12
findstructralnz(::Union{LinearAlgebra.SymTridiagonal, LinearAlgebra.Tridiagonal})
@ ArrayInterface ~/.julia/packages/ArrayInterface/9um5f/src/ArrayInterface.jl:349
findstructralnz(::LinearAlgebra.Bidiagonal)
@ ArrayInterface ~/.julia/packages/ArrayInterface/9um5f/src/ArrayInterface.jl:341
I think I must be misunderstanding the arguments to PrecomputedJacobianColorvec
. It seems that jac_prototype
should be accepting an example Jacobian of your base type (so a Matrix{Float64}
in my case) and the colors that you want. The example uses Symbolics to calculate the prototype jacobian, but I don’t think that’ll work with the iterative solvers in the objective function I want this to work with. My thought is to get a rough idea of the sparsity of my function by running ForwardDiff (or PolyesterForwardDiff) a couple of times. Right now, I know the sparsity for the problem I’m trying to get working, but I want to use the method I’m going to use.
For some background, obj
is a functor with methods obj(y, x)
, and obj(x)
which calculate my constraints and objective respectively. i.e.
function (obj::uo.ObjectiveFunction)(g, x)
f = x[1]^2 - x[2] #Todo: I don't think this is needed..
# g[1] = x[2] - 2*x[1]
g[1] = - 2*x[1] #Removing x[2] to test sparsity
g[2] = -x[2]
g[3] = x[1]^2
return f
end
function (obj::uo.ObjectiveFunction)(x)
f = x[1]^2 - x[2]
return f
end
I have those two functions because I’m going to use reverse mode for the gradient of the objective and then I want to use sparse PolyesterForwardDiff for the jacobian of the constraints.