I'm trying to estimate a weighted logistic regression as part of a bigger project. I have an implementation in Matlab2015b, but I wanted to give Julia a try to see if I could speed up the estimation and learn Julia better. The estimation time is on the project level already running into hours as this function is repeatedly called. I'm suspecting that due to my Matlab background, I'm missing out on some optimizations available in Julia.

function logit(y,x,w)
b = x\y;  # use ols values as start
maxit = 100;
tol = 0.000001;

#initialize breaks
crit = 1.0;
iter = 1;

  # loop
  while (iter < maxit) && (crit > tol)
          # Gradient & Hessian for logit
          (g,H) = LogitGradientAndHessian(y,x,b,w);

          # proposed change in coefficients
          db = -H\g

          # stepsize determination. Try steps of 1 or smaller
          s = 1.0
          L1 = LogitLogLikelihood(b+s*db,y,x,w)
          while s>tol
              s = s/2
              L2 = LogitLogLikelihood(b+s*db,y,x,w)
              if (L2-L1)<0 # log likelihood increases by less than 0/tol, abort
                  s = 2*s # take previous step size, with step size the likelihood is declining
              L1 = L2; # reduce step-size, this is the new log likelihood

      # take step
      b = b + s*db; # update coefficients
      crit = maximum(abs.(db)); # maximum absolute change in coefficient values (without step size like in original code...)
      iter = iter + 1;
end # end of while

return b,iter

function sigmoid(z::Float64)
z = 1.0/(1.0+exp(-z))
if z<0.0000001
  z = 0.0000001
elseif z>0.99999999
  z = 0.99999999
return z

function LogitLogLikelihood(b,y,x,w)
xb = x*b
L = sum(w .* (y.*xb-log.(1+exp.(xb))))
return L

function LogitGradientAndHessian(y,x,b,w)
# gradient
delta = sigmoid.(x*b)
g = x'*(w.*(y.-delta))

# pre-allocate Hessian
k = size(x,2)
H = zeros(eltype(g), (k,k))

# compute Hessian (could also do only upper-right and copy to bottom-left)
for n=1:size(x,1)
tmp = w[n]*delta[n]*(1.0-delta[n])
for kk1=1:k
    for kk2=1:k
          H[kk1,kk2] = H[kk1,kk2]-tmp*x[n,kk1]*x[n,kk2]

return g,H

The code to benchmark this (of course first do the warm up as per the instructions)

N = Int(1e6);
y = Int.(randn(N,1).>0)
x = [ones(N,1) randn(N,1)]
w = ones(N,1)
w = w./sum(w)

@time logit(y,x,w)

I find that at the moment the Julia code is about as fast as the Matlab code which follows a very similar style. Can this code be (substantially) improved?

  • \$\begingroup\$ How much of a difference does giving specific types to the variables for your functions? \$\endgroup\$ – Oscar Smith Aug 28 '17 at 4:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.