Using `Optim` package for a non-linear object function is too slow

I want to improve the performance for optimizing a function.

I use Optim package for optimizing a non-linear function with BFGS algorithm.

I put an object function (non-linear due to link_approx, which generates a cubic spline)

and its gradient vector into optimize.

However, it turned out to be 4 times slower than an R-programming complement.

I set the tolerance (Criteria for convergence) as the same as R.

(I can attach the code if it is needed)

Z::Matrix{Float64}; X::Matrix{Float64}; Y::Matrix{Float64}; B::Matrix{Float64}; G::Matrix{Float64}


using Splines2
using LinearAlgebra
using Optim

function link_approx(x_v::Array)
    local est; local der
    est = bs(x_v, knots = knots, order = 4)[:, 3:end-3] * fit[:theta]
    der = bs(x_v, knots = knots, order = 3)[:, 3:end-3] * coef
    return Dict{Symbol, Array{Float64}}(:est => est, :der => der)
end

@time for j in 1:r

# for update G

    function grad!(storage, gamma)
      local linkfit
      linkfit = link_approx(Y*gamma)

      output = (transpose(Y) * ((X*B[:,j] + linkfit[:est] - Z[:,j]) .* linkfit[:der])./n - U0[:,j] - U2[:,j] - U3[:,j] 
+ rho*(pennum * gamma - C0[:,j] - C2[:,j] - C3[:,j]))

      for i in 1:size(Y)[2]
        storage[i] = output[i]
      end
    end
  
    function obj(gamma)
      return norm(Z[:,j] - X*B[:,j] - link_approx(Y*gamma)[:est], 2)^2/(2*n) - transpose(U0[:,j] + U2[:,j] + U3[:,j])*(gamma) 
+ rho*(norm(gamma - C0[:,j], 2)^2 + norm(gamma - C2[:,j], 2)^2*lowrank_G + norm(gamma - C3[:,j], 2)^2*sparse_G)/2
    end
  
    temp = optimize(obj, grad!, G[:,j], BFGS(), Optim.Options(iterations = Int(5e1)))
    G[:,j] = Optim.minimizer(temp)

end

2.419329 seconds (32.44 M allocations: 824.036 MiB, 3.52% gc time, 3.57% compilation time)

(the gradient is calculated by the formula of derivatives of a B-spline Curve)

I think there is a problem with its gradient vector or duplicated compiling.

I don't know how to put value on a storage of gradient in a high dimension case.

Since its dimension is over 100, I used for loop.