# How to implement AD for this function?

**URL:** https://discourse.julialang.org/t/how-to-implement-ad-for-this-function/60623
**Category:** General Usage
**Tags:** differentiation, autodiff
**Created:** [May 6, 2021, 6:13am UTC](https://discourse.julialang.org/t/how-to-implement-ad-for-this-function/60623 "2021-05-06T06:13:18Z")
**Posts on this page:** 10
**Page:** 1

<div class="post-metadata">

### Author: ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [May 6, 2021, 6:13am UTC](https://discourse.julialang.org/t/how-to-implement-ad-for-this-function/60623/1 "2021-05-06T06:13:19Z")

</div>

I am solving a system of equations with the [iterative proportional fitting procedure](https://en.wikipedia.org/wiki/Iterative_proportional_fitting). That means that at each iteration I solve for a subset of the variables `μ0m`, then hold those fixed and solve for the other subset of variables `μ0f`. I am storing the residuals of the system in preallocated arrays `resm` and `resf`. I programmed a function `sys!` that computes the residuals of the system for each of those subsets:

```julia
const szA = 9
const szE = 2
const szG = 2
function sys!(resm, resf, μ0m, μ0f; Sm=Sm, Sf=Sf, Π=Π)
    μ0 = reshape(vcat(vec(μ0m), vec(μ0f)), (szA, szE, szG))
    μ0m′ = μ0[:, :, 1]
    μ0f′ = μ0[:, :, 2]

    μ0m′[1, szE] = 0.0
    μ0f′[1, szE] = 0.0

    prodsummand = similar(μ0, (szA, szE, szA, szE))
    prodsummand!(prodsummand, μ0m′, μ0f′; Sm=Sm, Sf=Sf, Π=Π)

    for idx in CartesianIndices(resm)
        am, em = Tuple(idx)
        resm[am, em] = Sm[am, em] - abs(μ0m′[am,em]) - sum( prodsummand[am,em,af,ef] for af in 1:szA, ef in 1:szE )
    end
    for idx in CartesianIndices(resf)
        af, ef = Tuple(idx)
        resf[af, ef] = Sf[af, ef] - abs(μ0f′[af,ef]) - sum( prodsummand[am,em,af,ef] for am in 1:szA, em in 1:szE )
    end
    return nothing
end

function prodsummand!(prodsum, μ0m′, μ0f′; Sm=Sm, Sf=Sf, Π=Π)
    for idx in CartesianIndices(prodsum)
        am, em, af, ef = Tuple(idx)
        prodsum[am, em, af, ef] = Π[am,em,af,ef]*sqrt(Sm[am,em]*Sf[af,ef])*prod( auxprod((μ0m′[am+k,em]*μ0f′[af+k,ef])/(Sm[am+k,em]*Sf[af+k,ef]))^(0.5*(β*(1-δ))^k) for k in 0:(maxT(am,af; Z=szA)-1) )
    end
    return prodsum
end

function maxT(am, af; Z=szA)
    Z - max(am, af) + 1
end

auxprod(x) = isfinite(abs(x)) ? abs(x) : 1.0

```

To solve the system, I am using `NLsolve` with finite differences, and it finds a solution. I am now wanting to speed up the computation by using automatic differentiation. I am having trouble specifying the Jacobian with AD. This is the algorithm that solves the system with finite diff:

```julia
using NLsolve
using LinearAlgebra

function ipfp(init, sys!; Sm=Sm, Sf=Sf, Π=Π, method=:trust_region, autodiff=:finite, outmaxiter=5_000, outtol=1e-8, opts...)
    μ0f0 = vec(init)

    μ0m1 = similar(init)
    μ0f1 = similar(init)
    d = Inf
    i = 0
    while (i <= outmaxiter) && (outtol < d)
        solm = nlsolve( (resm, μ0m) -> sys!(resm, resf, μ0m, μ0f0; Sm=Sm, Sf=Sf, Π=Π), init; method=method, autodiff=autodiff, opts... )
        if converged(solm)
            μ0m1 = abs.(vec(solm.zero))
        else
            println("solving for μ0m given μ0f did not converge at iteration ", i)
            return nothing
        end

        solf = nlsolve( (resf, μ0f) -> sys!(resm, resf, μ0m1, μ0f; Sm=Sm, Sf=Sf, Π=Π), init; method=method, autodiff=autodiff, opts... )
        if converged(solf)
            μ0f1 = abs.(vec(solf.zero))
        else
            println("solving for μ0f given μ0m did not converge at iteration ", i)
            return nothing
        end
        d = norm(μ0f1 - μ0f0)
        μ0f0 = copy(μ0f1)
        i += 1
    end
    ret = reshape(vcat(μ0m1, μ0f1), (szA, szE, szG))
    ret[1, szE, 1] = 0.0
    ret[1, szE, 2] = 0.0
    return ret
end

```

and here are some objects that complete a working example:

```julia
Sm = rand(szA, szE)
Sf = rand(szA, szE)
Π = rand(szA, szE, szA, szE)
resm = similar(Sm)
resf = similar(Sf)
init = rand(szA, szE)

```

I tried adding AD-able functions to the algorithm like

```julia
    lfunm! = (resm, μ0m) -> sys!(resm, resf, μ0m, μ0f0; Sm=Sm, Sf=Sf, Π=Π)
        sys_jacm! = ForwardDiff.jacobian!(lfunm!, resm, init)

```

but I get an error:

```julia
ERROR: MethodError: no method matching jacobian!(::var"#94#96"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, ::Matrix{Float64}, ::Matrix{Float64})

```

How can I specify the (partial) Jacobian of `sys!` using AD to pass it to `NLsolve`?

---

<div class="post-metadata">

### Author: ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [May 6, 2021, 6:23am UTC](https://discourse.julialang.org/t/how-to-implement-ad-for-this-function/60623/2 "2021-05-06T06:23:49Z")

</div>

I believe you accidentally swapped the first and second argument of [`jacobian!()`](https://juliadiff.org/ForwardDiff.jl/latest/user/api.html#ForwardDiff.jacobian!).

---

<div class="post-metadata">

### Author: ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [May 6, 2021, 6:31am UTC](https://discourse.julialang.org/t/how-to-implement-ad-for-this-function/60623/3 "2021-05-06T06:31:26Z")

</div>

With this algorithm

```julia
function ipfp_ad(init, sys!; Sm=Sm, Sf=Sf, Π=Π, method=:trust_region, outmaxiter=5_000, outtol=1e-8, opts...)
    μ0f0 = vec(init)

    μ0m1 = similar(init)
    μ0f1 = similar(init)
    d = Inf
    i = 0
    while (i <= outmaxiter) && (outtol < d)
        lfunm! = (resm, μ0m) -> sys!(resm, resf, μ0m, μ0f0; Sm=Sm, Sf=Sf, Π=Π)
        sys_jacm! = ForwardDiff.jacobian!(resm, lfunm!, init)
        solm = nlsolve( lfunm!, sys_jacm!, init; method=method, opts... )
        if converged(solm)
            μ0m1 = abs.(vec(solm.zero))
        else
            println("solving for μ0m given μ0f did not converge at iteration ", i)
            return nothing
        end

        lfunm! = (resf, μ0f) -> sys!(resm, resf, μ0m1, μ0f; Sm=Sm, Sf=Sf, Π=Π)
        sys_jacf! = ForwardDiff.jacobian!(resf, lfunf!, μ0f)
        solf = nlsolve( lfunf!, sys_jacf!, init; method=method, opts... )
        if converged(solf)
            μ0f1 = abs.(vec(solf.zero))
        else
            println("solving for μ0f given μ0m did not converge at iteration ", i)
            return nothing
        end
        d = norm(μ0f1 - μ0f0)
        μ0f0 = copy(μ0f1)
        i += 1
    end
    ret = reshape(vcat(μ0m1, μ0f1), (szA, szE, szG))
    ret[1, szE, 1] = 0.0
    ret[1, szE, 2] = 0.0
    return ret
end

```

I still get an error:

```julia
ERROR: MethodError: no method matching (::var"#71#73"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)})(::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{var"#71#73"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9}})

```

---

<div class="post-metadata">

### Author: ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [May 6, 2021, 10:35am UTC](https://discourse.julialang.org/t/how-to-implement-ad-for-this-function/60623/4 "2021-05-06T10:35:41Z")

</div>

` sys_jacf!` is not a function, it is an array. However, nlsolve expects a function to calculate the jacobian. Please also provide the full stack trace. From the error it impossible to figure out on which line of your code the problem is

---

<div class="post-metadata">

### Author: ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [May 6, 2021, 11:00am UTC](https://discourse.julialang.org/t/how-to-implement-ad-for-this-function/60623/5 "2021-05-06T11:00:51Z")

</div>

Here I post the full script:

```julia
using NLsolve
using LinearAlgebra
using ForwardDiff

const szA = 9
const szE = 2
const szG = 2

function sys!(resm, resf, μ0m, μ0f; Sm=Sm, Sf=Sf, Π=Π)
    μ0 = reshape(vcat(vec(μ0m), vec(μ0f)), (szA, szE, szG))
    μ0m′ = μ0[:, :, 1]
    μ0f′ = μ0[:, :, 2]

    μ0m′[1, szE] = 0.0
    μ0f′[1, szE] = 0.0

    prodsummand = similar(μ0, (szA, szE, szA, szE))
    prodsummand!(prodsummand, μ0m′, μ0f′; Sm=Sm, Sf=Sf, Π=Π)

    for idx in CartesianIndices(resm)
        am, em = Tuple(idx)
        resm[am, em] = Sm[am, em] - abs(μ0m′[am,em]) - sum( prodsummand[am,em,af,ef] for af in 1:szA, ef in 1:szE )
    end
    for idx in CartesianIndices(resf)
        af, ef = Tuple(idx)
        resf[af, ef] = Sf[af, ef] - abs(μ0f′[af,ef]) - sum( prodsummand[am,em,af,ef] for am in 1:szA, em in 1:szE )
    end
    return nothing
end

function prodsummand!(prodsum, μ0m′, μ0f′; Sm=Sm, Sf=Sf, Π=Π)
    for idx in CartesianIndices(prodsum)
        am, em, af, ef = Tuple(idx)
        prodsum[am, em, af, ef] = Π[am,em,af,ef]*sqrt(Sm[am,em]*Sf[af,ef])*prod( auxprod((μ0m′[am+k,em]*μ0f′[af+k,ef])/(Sm[am+k,em]*Sf[af+k,ef]))^(0.5*(β*(1-δ))^k) for k in 0:(maxT(am,af; Z=szA)-1) )
    end
    return prodsum
end

function maxT(am, af; Z=szA)
    Z - max(am, af) + 1
end

auxprod(x) = isfinite(abs(x)) ? abs(x) : 1.0

function ipfp_ad(init, sys!; Sm=Sm, Sf=Sf, Π=Π, method=:trust_region, autodiff=:finite, outmaxiter=5_000, outtol=1e-8, opts...)
    μ0f0 = vec(init)

    μ0m1 = similar(init)
    μ0f1 = similar(init)
    d = Inf
    i = 0
    while (i <= outmaxiter) && (outtol < d)
        lfunm! = (resm, μ0m) -> sys!(resm, resf, μ0m, μ0f0; Sm=Sm, Sf=Sf, Π=Π)
        sys_jacm! = ForwardDiff.jacobian!(resm, lfunm!, init)
        solm = nlsolve( (resm, μ0m) -> sys!(resm, resf, μ0m, μ0f0; Sm=Sm, Sf=Sf, Π=Π), init; method=method, autodiff=autodiff, opts... )
        if converged(solm)
            μ0m1 = abs.(vec(solm.zero))
        else
            println("solving for μ0m given μ0f did not converge at iteration ", i)
            return nothing
        end

        lfunf! = (resf, μ0f) -> sys!(resm, resf, μ0m0, μ0f; Sm=Sm, Sf=Sf, Π=Π)
        sys_jacf! = ForwardDiff.jacobian!(resf, lfunf!, μ0m1)
        solf = nlsolve( (resf, μ0f) -> sys!(resm, resf, μ0m1, μ0f; Sm=Sm, Sf=Sf, Π=Π), init; method=method, autodiff=autodiff, opts... )
        if converged(solf)
            μ0f1 = abs.(vec(solf.zero))
        else
            println("solving for μ0f given μ0m did not converge at iteration ", i)
            return nothing
        end
        d = norm(μ0f1 - μ0f0)
        μ0f0 = copy(μ0f1)
        i += 1
    end
    ret = reshape(vcat(μ0m1, μ0f1), (szA, szE, szG))
    ret[1, szE, 1] = 0.0
    ret[1, szE, 2] = 0.0
    return ret
end

Sm = rand(szA, szE)
Sf = rand(szA, szE)
Π = rand(szA, szE, szA, szE)
resm = similar(Sm)
resf = similar(Sf)
init = rand(szA, szE)
const β = 0.92
const δ = 0.10

```

Here is the error

```julia
julia> ipfp_sol = ipfp_ad(init, sys!; Sm=Sm, Sf=Sf, store_trace=true, extended_trace=true, method=:trust_region)
ERROR: MethodError: no method matching (::var"#15#19"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)})(::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{var"#15#19"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9}})
Closest candidates are:
  (::var"#15#19")(::Any, ::Any) at /Users/amrods/discourseQ4.jl:53
Stacktrace:
 [1] chunk_mode_jacobian!(result::Matrix{Float64}, f::var"#15#19"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, x::Matrix{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#15#19"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{var"#15#19"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9}}})
   @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/jacobian.jl:223
 [2] jacobian!(result::Matrix{Float64}, f::Function, x::Matrix{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#15#19"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{var"#15#19"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9}}}, ::Val{true})
   @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/jacobian.jl:60
 [3] jacobian!(result::Matrix{Float64}, f::Function, x::Matrix{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#15#19"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{var"#15#19"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9}}}) (repeats 2 times)
   @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/jacobian.jl:56
 [4] ipfp_ad(init::Matrix{Float64}, sys!::typeof(sys!); Sm::Matrix{Float64}, Sf::Matrix{Float64}, Π::Array{Float64, 4}, method::Symbol, autodiff::Symbol, outmaxiter::Int64, outtol::Float64, opts::Base.Iterators.Pairs{Symbol, Bool, Tuple{Symbol, Symbol}, NamedTuple{(:store_trace, :extended_trace), Tuple{Bool, Bool}}})
   @ Main ~/discourseQ4.jl:54
 [5] top-level scope
   @ REPL[1]:1

```

Line 53 is where I define `lfunm!`.

---

<div class="post-metadata">

### Author: ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [May 6, 2021, 11:08am UTC](https://discourse.julialang.org/t/how-to-implement-ad-for-this-function/60623/6 "2021-05-06T11:08:58Z")

</div>

> [@josuagrw](#):
>
> `sys_jacf!` is not a function, it is an array.

so this should do it?

```julia
sys_jacm! = (resm, x) -> ForwardDiff.jacobian!(resm, lfunm!, x)

```

---

<div class="post-metadata">

### Author: ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [May 6, 2021, 11:39am UTC](https://discourse.julialang.org/t/how-to-implement-ad-for-this-function/60623/7 "2021-05-06T11:39:36Z")

</div>

I changed the lines defining the Jacobian:

```julia
using NLsolve
using LinearAlgebra
using ForwardDiff

const szA = 9
const szE = 2
const szG = 2
const β = 0.92
const δ = 0.10

function sys!(resm, resf, μ0m, μ0f; Sm=Sm, Sf=Sf, Π=Π)
    μ0 = reshape(vcat(vec(μ0m), vec(μ0f)), (szA, szE, szG))
    μ0m′ = μ0[:, :, 1]
    μ0f′ = μ0[:, :, 2]

    μ0m′[1, szE] = 0.0
    μ0f′[1, szE] = 0.0

    prodsummand = similar(μ0, (szA, szE, szA, szE))
    prodsummand!(prodsummand, μ0m′, μ0f′; Sm=Sm, Sf=Sf, Π=Π)

    for idx in CartesianIndices(resm)
        am, em = Tuple(idx)
        resm[am, em] = Sm[am, em] - abs(μ0m′[am,em]) - sum( prodsummand[am,em,af,ef] for af in 1:szA, ef in 1:szE )
    end
    for idx in CartesianIndices(resf)
        af, ef = Tuple(idx)
        resf[af, ef] = Sf[af, ef] - abs(μ0f′[af,ef]) - sum( prodsummand[am,em,af,ef] for am in 1:szA, em in 1:szE )
    end
    return nothing
end

function prodsummand!(prodsum, μ0m′, μ0f′; Sm=Sm, Sf=Sf, Π=Π)
    for idx in CartesianIndices(prodsum)
        am, em, af, ef = Tuple(idx)
        prodsum[am, em, af, ef] = Π[am,em,af,ef]*sqrt(Sm[am,em]*Sf[af,ef])*prod( auxprod((μ0m′[am+k,em]*μ0f′[af+k,ef])/(Sm[am+k,em]*Sf[af+k,ef]))^(0.5*(β*(1-δ))^k) for k in 0:(maxT(am,af; Z=szA)-1) )
    end
    return prodsum
end

function maxT(am, af; Z=szA)
    Z - max(am, af) + 1
end

auxprod(x) = isfinite(abs(x)) ? abs(x) : 1.0

function ipfp_ad(init, sys!; Sm=Sm, Sf=Sf, Π=Π, method=:trust_region, outmaxiter=5_000, outtol=1e-8, opts...)
    μ0f0 = vec(init)

    μ0m1 = similar(init)
    μ0f1 = similar(init)
    d = Inf
    i = 0
    while (i <= outmaxiter) && (outtol < d)
        lfunm!(resx, x) = sys!(resx, resf, x, μ0f0; Sm=Sm, Sf=Sf, Π=Π)
        sys_jacm!(resx, x) = ForwardDiff.jacobian!(resx, lfunm!, x)
        solm = nlsolve( lfunm!, sys_jacm!, init; method=method, opts... )
        if converged(solm)
            μ0m1 = abs.(vec(solm.zero))
        else
            println("solving for μ0m given μ0f did not converge at iteration ", i)
            return nothing
        end

        lfunf!(resx, x) = sys!(resm, resx, μ0m1, x; Sm=Sm, Sf=Sf, Π=Π)
        sys_jacf!(resx, x) = ForwardDiff.jacobian!(resx, lfunf!, x)
        solf = nlsolve( lfunf!, sys_jacf!, init; method=method, opts... )
        if converged(solf)
            μ0f1 = abs.(vec(solf.zero))
        else
            println("solving for μ0f given μ0m did not converge at iteration ", i)
            return nothing
        end
        d = norm(μ0f1 - μ0f0)
        μ0f0 = copy(μ0f1)
        i += 1
    end
    ret = reshape(vcat(μ0m1, μ0f1), (szA, szE, szG))
    ret[1, szE, 1] = 0.0
    ret[1, szE, 2] = 0.0
    return ret
end

Sm = rand(szA, szE)
Sf = rand(szA, szE)
Π = rand(szA, szE, szA, szE)
resm = similar(Sm)
resf = similar(Sf)
init = rand(szA, szE)

sol = ipfp_ad(init, sys!; Sm=Sm, Sf=Sf, Π=Π, method=:trust_region, outmaxiter=5_000, outtol=1e-8)

```

I get this error:

```julia
ERROR: LoadError: MethodError: no method matching (::var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)})(::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9}})
Closest candidates are:
  (::var"#lfunm!#15")(::Any, ::Any) at /Users/amrods/discourseQ4.jl:55
Stacktrace:
  [1] chunk_mode_jacobian!(result::Matrix{Float64}, f::var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, x::Matrix{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/jacobian.jl:223
  [2] jacobian!(result::Matrix{Float64}, f::Function, x::Matrix{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9}}}, ::Val{true})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/jacobian.jl:60
  [3] jacobian!(result::Matrix{Float64}, f::Function, x::Matrix{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9}}}) (repeats 2 times)
    @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/jacobian.jl:56
  [4] (::var"#sys_jacm!#16"{var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}})(resx::Matrix{Float64}, x::Matrix{Float64})
    @ Main ~/discourseQ4.jl:56
  [5] (::NLSolversBase.var"#fj!#7"{var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, var"#sys_jacm!#16"{var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}}})(fx::Matrix{Float64}, jx::Matrix{Float64}, x::Matrix{Float64})
    @ NLSolversBase ~/.julia/packages/NLSolversBase/geyh3/src/objective_types/abstract.jl:7
  [6] value_jacobian!!(obj::OnceDifferentiable{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, F::Matrix{Float64}, J::Matrix{Float64}, x::Matrix{Float64})
    @ NLSolversBase ~/.julia/packages/NLSolversBase/geyh3/src/interface.jl:124
  [7] value_jacobian!!
    @ ~/.julia/packages/NLSolversBase/geyh3/src/interface.jl:122 [inlined]
  [8] trust_region_(df::OnceDifferentiable{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, initial_x::Matrix{Float64}, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, factor::Float64, autoscale::Bool, cache::NLsolve.NewtonTrustRegionCache{Matrix{Float64}})
    @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/solvers/trust_region.jl:119
  [9] trust_region (repeats 2 times)
    @ ~/.julia/packages/NLsolve/gJL1I/src/solvers/trust_region.jl:235 [inlined]
 [10] nlsolve(df::OnceDifferentiable{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, initial_x::Matrix{Float64}; method::Symbol, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::Static, linsolve::NLsolve.var"#27#29", factor::Float64, autoscale::Bool, m::Int64, beta::Int64, aa_start::Int64, droptol::Float64)
    @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:26
 [11] nlsolve(f::Function, j::Function, initial_x::Matrix{Float64}; inplace::Bool, kwargs::Base.Iterators.Pairs{Symbol, Symbol, Tuple{Symbol}, NamedTuple{(:method,), Tuple{Symbol}}})
    @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:67
 [12] ipfp_ad(init::Matrix{Float64}, sys!::typeof(sys!); Sm::Matrix{Float64}, Sf::Matrix{Float64}, Π::Array{Float64, 4}, method::Symbol, outmaxiter::Int64, outtol::Float64, opts::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Main ~/discourseQ4.jl:57
 [13] top-level scope
    @ ~/discourseQ4.jl:91
in expression starting at /Users/amrods/discourseQ4.jl:91

```

Line 55 is `lfunm!(resx, x) = sys!(resx, resf, x, μ0f0; Sm=Sm, Sf=Sf, Π=Π)`

---

<div class="post-metadata">

### Author: ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [May 7, 2021, 7:50am UTC](https://discourse.julialang.org/t/how-to-implement-ad-for-this-function/60623/8 "2021-05-07T07:50:56Z")

</div>

yes that looks good

---

<div class="post-metadata">

### Author: ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [May 7, 2021, 7:58am UTC](https://discourse.julialang.org/t/how-to-implement-ad-for-this-function/60623/9 "2021-05-07T07:58:33Z")

</div>

You’re still calling `jacobian!()` incorrectly. For your case you need to provide 4 arguments. [Please read the documentation—it is well explained there.](https://juliadiff.org/ForwardDiff.jl/stable/user/api/#ForwardDiff.jacobian!)

---

<div class="post-metadata">

### Author: ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [May 7, 2021, 9:03am UTC](https://discourse.julialang.org/t/how-to-implement-ad-for-this-function/60623/10 "2021-05-07T09:03:51Z")

</div>

yeah, i missed that bang at the end, see [A bug in ForwardDiff? - #4 by amrods](https://discourse.julialang.org/t/a-bug-in-forwarddiff/60691/4)