ForwardDiff.jl and its Dual Type question

General Usage
,
1

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)[1];
    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?

Thank you for your attention.

2

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(β))[1] 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]?

1 Like
3

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)[1];
    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
2 Likes
4

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
5

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!