MethodError: no method matching Float64 when using ForwardDiff.jl despite me pre-allocating arrays to take floating point

Specific Domains Optimization (Mathematical)
1 
using Optimization, OptimizationOptimJL, OptimizationBBO, OptimizationMOI, ForwardDiff, ModelingToolkit, Random, OptimizationOptimisers, Statistics, Plots

function MSEg(u, p)
  X = p.X
  y = p.y
  h = u[1]
  beta = u[2:end]
  N = length(y)
  Xbeta = X*beta

  krnl = u -> exp.(-(u .* u) ./ 2) ./ sqrt(2*π)

  yhat = Array{Float64, 2}(undef, N, 1)
  weight_i = Array{Float64, 2}(undef, N, 1)
  weight_sum = Array{Float64, 2}(undef, N, 1)

  for i = 1:N
    u = (Xbeta .- Xbeta[i]) ./ h
    Ku = krnl(u)
    weight_i[i] = Ku[i]
    weight_sum[i] = sum(Ku)
    yhat[i] = sum(Ku .* y) ./ sum(Ku)
  end

  residual = (y .- yhat) ./ (1 .- (weight_i ./ weight_sum))
  SBeta = sum(residual.^2)
end

Random.seed!(1)
N = 1000
X = -1 .+ 2 .* rand(N, 2)
beta = [1;3]
Xbeta = X*beta
y = Xbeta.^2 .* (1 .+ sin.(Xbeta)) .+ randn(N, 1)

y_ = y .- X[:,1]
X_ = X[:,2]
beta0 = X_'*X_\X_'*y_ 
beta0 = [1; beta0]

x_ = X*beta0
hx = median(abs.(x_ .- median(x_)))/0.6745*(4/3/N)^0.2
hy = median(abs.(y .- median(y)))/0.6745*(4/3/N)^0.2
h = sqrt(hy*hx)

param0 = [h;beta0]
u = param0

struct smol
  X
  y
end

p = smol(X, y)

cons(res, u, p) = (res .= u[2])

lb = [0;-Inf;-Inf]
ub = [Inf;Inf;Inf]

optf = OptimizationFunction(MSEg, Optimization.AutoForwardDiff(), cons = cons)
prob = OptimizationProblem(optf, u, p, lb = lb, ub = ub, lcons = 1, ucons = 1)
sol = solve(prob, NelderMead())

I had received this error in the past when I tried to use a for loop to populate elements of arrays and matrices. I had not had trouble in the past with this error in a while but now it is coming back again and I cannot understand why

2

despite me pre-allocating arrays to take floating point

This is the problem. You need:

function MSEg(u::AbstractVector{T}, p) where {T}
    # ...
     yhat = Array{T, 2}(undef, N, 1)
    #...

Your function must not assume that u is a Vector{Float64}.

3

I am still confused. I changed the function below and get the same error when running the script

function MSEg(u::AbstractVector{T}, p) where {T}
  typeof(u)
  X = p.X
  y = p.y
  h = u[1]
  beta = u[2:end]
  N = length(y)
  Xbeta = X*beta

  krnl = u -> exp.(-(u .* u) ./ 2) ./ sqrt(2*π)

  # yhat = Array{Float64, 2}(undef, N, 1)
  # weight_i = Array{Float64, 2}(undef, N, 1)
  # weight_sum = Array{Float64, 2}(undef, N, 1)

  yhat = Array{T, 2}(undef, N, 1)
  weight_i = Array{T, 2}(undef, N, 1)
  weight_sum = Array{T, 2}(undef, N, 1)

  for i = 1:N
    u = (Xbeta .- Xbeta[i]) ./ h
    Ku = krnl(u)
    weight_i[i] = Ku[i]
    weight_sum[i] = sum(Ku)
    yhat[i] = sum(Ku .* y) ./ sum(Ku)
  end
  residual = (y .- yhat) ./ (1 .- (weight_i ./ weight_sum))
  SBeta = sum(residual.^2)
end
4

These won’t be Float64, so the type is incorrect.

5

You can simplify your MSEg code a lot (and there were a bunch of other problems in your code, like the wrong optimizer). How about this:

using Optimization
import OptimizationMOI
import ForwardDiff
import Ipopt
import Random
using Statistics

kernel(u) = exp(-(u * u) / 2) / sqrt(2π)

function MSEg(u::AbstractArray{T}, p) where {T}
    X, y = p.X, p.y
    h, beta = u[1], u[2:end]
    Xbeta = X * beta
    residuals = zero(T)
    for i in 1:length(y)
        Ku = kernel.((Xbeta .- Xbeta[i]) ./ h)
        yhat = sum(Ku .* y) ./ sum(Ku)
        residuals += ((y[i] - yhat) / (1 - (Ku[i] / sum(Ku))))^2
    end
    return residuals
end

Random.seed!(1)
N = 1000
X = -1 .+ 2 .* rand(N, 2)
beta = [1;3]
Xbeta = X*beta
y = Xbeta.^2 .* (1 .+ sin.(Xbeta)) .+ randn(N, 1)

y_ = y .- X[:,1]
X_ = X[:,2]
beta0 = X_' * X_ \ X_' * y_ 
beta0 = [1; beta0]

x_ = X * beta0
hx = median(abs.(x_ .- median(x_))) / 0.6745 * (4 / 3 / N)^0.2
hy = median(abs.(y .- median(y))) / 0.6745 * (4 / 3  / N)^0.2
h = sqrt(hy*hx)

struct smol
  X
  y
end

p = smol(X, y)

cons(res, u, p) = (res .= u[2])

optf = OptimizationFunction(MSEg, Optimization.AutoForwardDiff(), cons = cons)
prob = OptimizationProblem(
  optf, 
  [h; beta0], 
  p; 
  lb = [0, -Inf, -Inf],
  ub = [Inf, Inf, Inf],
  lcons = [1.0],
  ucons = [1.0],
)
sol = solve(prob, Ipopt.Optimizer())