-
-
Save tyleransom/243d423290f845dcc2ea306d54f7f56e to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
using BenchmarkTools | |
using Optim, NLSolversBase | |
function datagen(n=5500,k=13) | |
b = rand(k) | |
l = rand(2) | |
x = randn(n,k) | |
abil = randn(n,1) | |
vabil = randn(n,1) | |
sector = 1.0.*(rand(n,1).>0.8) | |
weight = rand(n,1) | |
y = l[1] .+ l[2].*(x*b .+ abil .+ vabil.^2) | |
return y,x,abil,vabil,sector,weight | |
end | |
# Objective function A | |
function wolslambdaNoAssignCorrFusion(b::Vector,y::Array,x::Array,abil::Array,vabil::Array,sector::Array,weight::Array) | |
# λ0 = b[end-1] | |
# λ1 = b[end] | |
# β = b[1:end-2] | |
# sector=1 corresponds to subset of observations that get hit by lambda | |
# sector=0 corresponds to all other observations | |
if ~isequal(size(y),size(sector)) | |
error("sector variable is wrong") | |
end | |
SSE = sum(weight.*(vabil.*(((sector.==0).*1).+((sector.==1).*b[end-1])).^2 .+ (y.-((sector.==0).*0).-((sector.==1).*b[end-1]).-(((sector.==0).*1).+((sector.==1).*b[end])).*(x*b[1:end-2].+abil)).^2)) | |
return SSE | |
end | |
# Objective function B | |
function wolslambdaNoAssign(b::Vector,y::Array,x::Array,abil::Array,vabil::Array,sector::Array,weight::Array) | |
# λ0 = b[end-1] | |
# λ1 = b[end] | |
# β = b[1:end-2] | |
# sector=1 corresponds to subset of observations that get hit by lambda | |
# sector=0 corresponds to all other observations | |
if ~isequal(size(y),size(sector)) | |
error("sector variable is wrong") | |
end | |
SSE = sum(weight.*(vabil.*(((sector.==0)*1).+((sector.==1)*b[end-1])).^2 .+ (y.-((sector.==0)*0).-((sector.==1)*b[end-1]).-(((sector.==0)*1).+((sector.==1)*b[end])).*(x*b[1:end-2].+abil)).^2)) | |
return SSE | |
end | |
# Objective function C | |
function wolslambda(b::Vector,y::Array,x::Array,abil::Array,vabil::Array,sector::Array,weight::Array) | |
λ0 = b[end-1] | |
λ1 = b[end] | |
β = b[1:end-2] | |
# sector=1 corresponds to subset of observations that get hit by lambda | |
# sector=0 corresponds to all other observations | |
if ~isequal(size(y),size(sector)) | |
error("sector variable is wrong") | |
end | |
SSE = sum(weight.*(vabil.*(((sector.==0)*1).+((sector.==1)*λ1)).^2 .+ (y.-((sector.==0)*0).-((sector.==1)*λ0).-(((sector.==0)*1).+((sector.==1)*λ1)).*(x*β.+abil)).^2)) | |
return SSE | |
end | |
# Objective function D | |
function wolslambdaNoAssignLooping(b::Vector,y::Array,x::Array,abil::Array,vabil::Array,sector::Array,weight::Array) | |
# λ0 = b[end-1] | |
# λ1 = b[end] | |
# β = b[1:end-2] | |
# sector=1 corresponds to subset of observations that get hit by lambda | |
# sector=0 corresponds to all other observations | |
if ~isequal(size(y),size(sector)) | |
error("sector variable is wrong") | |
end | |
T = promote_type(eltype(b), eltype(x)) | |
SSE = zero(T) | |
for i=1:size(y,1) | |
SSE+=weight[i]*(vabil[i]*(((sector[i]==0)*1)+((sector[i]==1)*b[end-1]))^2 + (y[i]-((sector[i]==0)*0)-((sector[i]==1)*b[end-1])-(((sector[i]==0)*1)+((sector[i]==1)*b[end]))*(dot(x[i,:],b[1:end-2])+abil[i]))^2) | |
end | |
return SSE | |
end | |
# Objective function D | |
function wolsNoAlloc(b::Vector,y::Array,x::Array,abil::Array,vabil::Array,sector::Array,weight::Array) | |
# λ0 = b[end-1] | |
# λ1 = b[end] | |
# β = b[1:end-2] | |
# sector=1 corresponds to subset of observations that get hit by lambda | |
# sector=0 corresponds to all other observations | |
if ~isequal(size(y),size(sector)) | |
error("sector variable is wrong") | |
end | |
T = promote_type(eltype(b), eltype(x)) | |
SSE = zero(T) | |
@inbounds for i=1:size(y,1) | |
b_dot_x = zero(T) | |
for j in 1:size(x, 1) | |
b_dot_x += x[j, i] * b[j] | |
end | |
SSE+=weight[i]*(vabil[i]*(((sector[i]==0)*1)+((sector[i]==1)*b[end-1]))^2 + (y[i]-((sector[i]==0)*0)-((sector[i]==1)*b[end-1])-(((sector[i]==0)*1)+((sector[i]==1)*b[end]))*(b_dot_x+abil[i]))^2) | |
end | |
return SSE | |
end | |
function runtests() | |
# data (arrays): | |
wageg,xg,abilg,vabilg,gflg,qg = datagen() | |
bstartg = rand(size(xg,2)) | |
lambdag0start = 0.5 | |
lambdag1start = 0.75 | |
# starting value | |
x0 = vec([bstartg;lambdag0start;lambdag1start]) | |
fung = TwiceDifferentiable((arg)->wolslambda(arg,wageg,xg,abilg,vabilg,gflg,qg), x0; autodiff = :forward); | |
println("\n original") | |
res = optimize(fung, x0, LBFGS(), Optim.Options(show_trace=false,iterations=100_000,g_tol=1e-6,f_tol=1e-20)) | |
@time res = optimize(fung, x0, LBFGS(), Optim.Options(show_trace=false,iterations=100_000,g_tol=1e-6,f_tol=1e-20)) | |
println(res.minimizer) | |
fung = TwiceDifferentiable((arg)->wolslambdaNoAssign(arg,wageg,xg,abilg,vabilg,gflg,qg), x0; autodiff = :forward); | |
println("\n get rid of mutliple copying") | |
res = optimize(fung, x0, LBFGS(), Optim.Options(show_trace=false,iterations=100_000,g_tol=1e-6,f_tol=1e-20)) | |
@time res = optimize(fung, x0, LBFGS(), Optim.Options(show_trace=false,iterations=100_000,g_tol=1e-6,f_tol=1e-20)) | |
println(res.minimizer) | |
fung = TwiceDifferentiable((arg)->wolslambdaNoAssignCorrFusion(arg,wageg,xg,abilg,vabilg,gflg,qg), x0; autodiff = :forward); | |
println("\n get rid of multiple copying and use correct fusion operators") | |
res = optimize(fung, x0, LBFGS(), Optim.Options(show_trace=false,iterations=100_000,g_tol=1e-6,f_tol=1e-20)) | |
@time res = optimize(fung, x0, LBFGS(), Optim.Options(show_trace=false,iterations=100_000,g_tol=1e-6,f_tol=1e-20)) | |
println(res.minimizer) | |
fung = TwiceDifferentiable((arg)->wolslambdaNoAssignLooping(arg,wageg,xg,abilg,vabilg,gflg,qg), x0; autodiff = :forward); | |
println("\n loop over observations") | |
res = optimize(fung, x0, LBFGS(), Optim.Options(show_trace=false,iterations=100_000,g_tol=1e-6,f_tol=1e-20)) | |
@time res = optimize(fung, x0, LBFGS(), Optim.Options(show_trace=false,iterations=100_000,g_tol=1e-6,f_tol=1e-20)) | |
println(res.minimizer) | |
xg_transpose = Matrix(xg') | |
fung = TwiceDifferentiable((arg)->wolsNoAlloc(arg,wageg,xg_transpose,abilg,vabilg,gflg,qg), x0; autodiff = :forward); | |
println("\n no-alloc") | |
res = optimize(fung, x0, LBFGS(), Optim.Options(show_trace=false,iterations=100_000,g_tol=1e-6,f_tol=1e-20)) | |
@time res = optimize(fung, x0, LBFGS(), Optim.Options(show_trace=false,iterations=100_000,g_tol=1e-6,f_tol=1e-20)) | |
println(res.minimizer) | |
end | |
runtests() | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment