# ForwardDiff.jl and its Dual Type question

Hello fellow Julia users,

I want to ask a question using ForwardDiff.jl.

``````using LinearAlgebra, Statistics, Distributions, Random
using Optim, ForwardDiff
using Optim: converged, maximum, maximizer, minimizer, iterations

#%% DGP
Random.seed!(1)
N  = 100000;
β0 = 0.15; β1 = 0.25; β2 = -1.5;
σ = 3.0;
x1 = randn(N,1);
x2 = randn(N,1);
ε  = σ * randn(N,1);
ydata = β0 .+ β1.*x1 .+ β2.*x2 .+ ε;
X = [ones(N,1) x1 x2]; # Nx3 matrix
β_hat = inv(X'X)*X'ydata

#%%
function normaleqn(β,ydata,X)
N = size(X);
f = zeros(eltype(β),size(X));
A = zeros(eltype(β),size(β));
for i=1:N
f[i,:] = X[i,:] * (ydata[i] .- X[i,:]' * β) # 1x3 * (1x1)
end
A = 1/N*sum(f,dims=1)'
end
#%%
β = ones(3,1)
ortho(β) = normaleqn(β,ydata,X)
jaco = ForwardDiff.jacobian(β->ortho(β),zeros(3,1)) # for testing

obj_I(β) = ortho(β)' * I * ortho(β)
obj_I(ones(3,1)) # for testing

opt = optimize(obj_I,zeros(3,1),BFGS(),autodiff=:forward,Optim.Options(iterations=1200));
``````

The last line returns an error message,

``````**MethodError: no method matching extract_gradient!(::Type{ForwardDiff.Tag{typeof(obj_I),Float64}}, ::DiffResults.MutableDiffResult{1,Float64,Tuple{Array{Float64,2}}}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(obj_I),Float64},Float64,3},2})**
``````

I now at least have a sense that ForwardDiff.jl uses the Dual Type and so I need to provide `eltype()` to let Julia flexibly define the types for the returned object. So, I managed to get the jacobian of `ortho(β)`, as shown in the middle of the script. However, I still can’t get the optimization routine run. Can someone explain why this is happening, please?

For those of you who are reading this, I found that `obj_I(ones(3,1))` returns Array{Float64,2} type, so changing the script `obj_I(β) = (ortho(β)' * W * ortho(β))` solves the problem. However, the code looks ugly. Is there any way I can make it cleaner, i.e. letting the code spit out `Float64` right away without having to write `obj_I(β) = (ortho(β)' * W * ortho(β))`?

1 Like

After some hitting my head on the wall, I came across `vec()` which returns a `Vector{Float64}` type so this worked. I will leave this post for future reference. I wish there is some kind of structured tutorial on Julia so that I can understand the programming language in depth so I don’t have to be depressed with this kind of things…

``````function normaleqn2(β,ydata,X)
N = size(X);
f = zeros(eltype(β),size(X));
A = zeros(size(β));
for i=1:N
f[i,:] = (ydata[i] .- X[i,:]' * β)*X[i,:] # 1x3 * (1x1)
end
return A = vec(sum(f,dims=1)'/N)
end
``````
1 Like

Generally `vec` just converts to a vector, w/o affecting the type. Eg

``````julia> vec(fill(1, 2, 3))
6-element Array{Int64,1}:
1
1
1
1
1
1
``````

It is documented here.

1 Like

Thank you for your reply. That’s the method I used in the second reply to myself by using `vec(sum(f,dims=1)'/N)` to return a vector. I am glad I got confirmation from an expert!