Hi all,
Does anybody know why the following code produces different (incorrect) jacobians using ForwardDiff? The function being differentiated contains sparse matrices, but if I don’t expand them to dense matrices the jacobian from the sparse function is incorrect…
using SparseArrays, ForwardDiff, LinearAlgebra
# setup data
θ = [30.0, 70.0, 90.0]
x = [1.9873817034700316, 3.9747634069400632, 1.9873817034700316, 1.9873817034700316, 1.9873817034700316, 3.9747634069400632, 1.9873817034700316, 0.028391167192429023, 0.022082018927444796, 0.13249211356466878, 0.13249211356466878]
ν = [ 0.0, 0.0, 0.028391167192429023, 0.24290220820189273,0.7287066246056783, 0.0, 1.9873817034700316, 3.9747634069400632, 5.962145110410095, 7.9495268138801265,]
λ = zeros(24)
λ[end] = -1.9873817034700316
S = [
1.0 0.0 0.0 0.0 0.0 0.0 -1.0
0.0 1.0 0.0 0.0 0.0 0.0 -2.0
0.0 0.0 1.0 0.0 0.0 0.0 -1.0
0.0 0.0 0.0 1.0 0.0 0.0 -1.0
0.0 0.0 0.0 0.0 1.0 0.0 -1.0
0.0 0.0 0.0 0.0 0.0 1.0 -2.0
]
se_vars = ((2, -1.0), (3, -1.0),(3, -3.0), (1, -2.0), (1, -1.0),)
Cp = spzeros(11)
Cp[8] = 1; Cp[9] = 2; Cp[10] = 3; Cp[11] = 4
clb, cub = ([0.0], [1.0])
xlb = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
xub = [100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 10.0, 10.0, 10.0, 10.0]
# create functions
Se(θ) = sparse(
[1, 2, 3, 4, 3],
[3, 4, 4, 5, 5],
[k / θ[idx] for (idx, k) in se_vars],
4,
7,
)
E(θ) = [
S spzeros(6, 4)
Se(θ) I(4)
]
d = _ -> spzeros(10)
_c = spzeros(11)
_c[7] = -1.0
c = _ -> _c
M = _ -> [
-I(11)
I(11)
-Cp'
Cp'
]
h = _ -> [-xlb; xub; -clb; cub]
Q = _ -> spzeros(11, 11)
xidxs = 1:length(x)
νidxs = last(xidxs) .+ (1:length(ν))
λidxs = last(νidxs) .+ (1:length(λ))
θidxs = last(λidxs) .+ (1:length(θ))
sparse_F(z) = [
Q(z[θidxs]) * z[xidxs] + c(z[θidxs]) -
E(z[θidxs])' * z[νidxs] - M(z[θidxs])' * z[λidxs]
E(z[θidxs]) * z[xidxs] - d(z[θidxs])
spdiagm(z[λidxs]) * (M(z[θidxs]) * z[xidxs] - h(z[θidxs]))
]
dense_F(z) = [
Array(Q(z[θidxs])) * z[xidxs] + Array(c(z[θidxs])) -
Array(E(z[θidxs]))' * z[νidxs] - Array(M(z[θidxs]))' * z[λidxs]
Array(E(z[θidxs])) * z[xidxs] - Array(d(z[θidxs]))
diagm(z[λidxs]) * (Array(M(z[θidxs])) * z[xidxs] - Array(h(z[θidxs])))
]
z = [x; ν; λ; θ]
dense_jac = ForwardDiff.jacobian(dense_F, z)
sparse_jac = ForwardDiff.jacobian(sparse_F, z)
# large differences
maximum(dense_jac - sparse_jac)
unique(dense_jac)
unique(sparse_jac)
Casting everything to dense arrays works fine, but it is quite slow, since the dimensions of my problem are large. The way I am constructing these sparse functions may seem contrived, but it is part of a package I am working on 
Does anybody have any ideas?
