Solving Logistic Regression estimation equations

Background

In logistic regression, a dichotomous variable y is modelled as follows

p(y|(Subscript[x, 1]),(Subscript[x, 2]),...,(Subscript[x, n]))= ((Exp[y <ϕ((Subscript[x, i])),θ>])/(1+Exp[<ϕ((Subscript[x, i])),θ>]))

Where ϕ is a predetermined  linear function of input and θ is a vector of parameters

This model is trained by maximum conditional likelihood, ie minimizing the following

O[θ]=-(∑_i)Log[p((y^i)|((Subscript[x, 1])^i),((Subscript[x, 2])^i),...,((Subscript[x, n])^i))]

The gradient of this objective function is

(∂_(Subscript[θ, j])) O = (∑_i)(Subscript[ϕ, j])((Subscript[x, i])) (((Exp[<ϕ((x^i)),θ>])/(1+Exp[<ϕ((x^i)),θ>]))-(y^i))

You can rewrite it as

(∂_(Subscript[θ, j])) O = N(∑_x)p(x)(Subscript[ϕ, j])(x) (((Exp[<ϕ(x),θ>])/(1+Exp[<ϕ(x),θ>]))-p(y=1|x))

Where N is the sample size, p (x) and p (y = 1 | x) are the empirical probabilities of given variable assignments. N can be dropped because it doesn't affect the solution

The Hessian will be as follows

(∂_((Subscript[θ, j]),(Subscript[θ, k])))O =(∑_x)p(x)(Subscript[ϕ, i])(x) (Subscript[ϕ, j])(x)((Exp[<ϕ(x),θ>])/(((1+Exp[<ϕ(x),θ>]))^2))

The Hessian is positive definite, so all we need is to solve for the point where gradient is 0

For simplicity, suppose p (x) is uniform

Solving the equations

Since (Subscript[ϕ, k]) (x) is a linear function, it can be represented as a matrix -- let k' th column of the matrix be the values of (Subscript[ϕ, k])(x) for different values of x in some order

If we order the x's in the same way as the binary encodings of integers go (00,01,10,11), and use features of the form (1,(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 1])(Subscript[x, 2])), we can write the design matrix for k variables as follows

k=2;
design=KroneckerProduct@@Table[{{1,0},{1,1}},{i,1,k}];
design//MatrixForm

params=(Subscript[θ, #])&/@Range[(2^k)]

We can now write i' th gradient equation as follows

eq[i_]:=Sum[Transpose[design][[i]][[x]](((Exp[design[[x]].params])/(Exp[design[[x]].params]+1))-g[x]),{x,1,(2^k)}]==0
eq[2]

solve the following system

eqs=eq[#]&/@Range[1,(2^k)]

Notice the equations are linear in terms of modelled p (y = 1 | x) probabilities so do variable substitution

odds[x_]:=((Exp[design[[x]].params])/(Exp[design[[x]].params]+1))

leftvars=odds[#]&/@Range[(2^k)]

rightvars=(Subscript[t, #])&/@Range[(2^k)]

subst=Thread[leftvars->rightvars]

gs=g[#]&/@Range[(2^k)] (* these are empirical probabilities of y=1 for each x *)

eqs/.subst//MatrixForm

You can see that the system in terms of new variables is trivial to solve

Solve[eqs/.subst,rightvars]

Now solve for θ' s in terms of t' s

Solve[subst/.Rule->Equal,params]

This works for 2 variables, but for 3 variables and more, Mathematica has trouble solving the equations. Notice however that you can solve this system by solving 1 equation at a time in the right order (solve for (Subscript[θ, 1]), substitute value into next equation, and repeat)

eqs2=subst/.Rule->Equal

tsol=MapThread[Solve,{eqs2,params}]//Flatten

Use these rules repeatedly to eliminate θ' s completely

tsol2=tsol[[All,2]]//.tsol//Simplify

solution=Thread[params->tsol2]

Below is a function that does the steps above. Note that it builds the design matrix explicitly, so it has  (2^(2k))complexity where k is the number of variables. The first set of equations (eqs) is trivial for saturated models, but it wouldn't be trivial for unsaturated models. For unsaturated models, the design matrix could be obtained by removing columns from the corresponding saturated design matrix. For instance for 2 variables, the "no 2 way interactions" matrix will look like this

saturated=KroneckerProduct@@Table[{{1,0},{1,1}},{i,1,2}];
saturated//MatrixForm

Transpose[Take[Transpose[saturated],{1,3}]]//MatrixForm

solvelog[k_]:=Module[{design,params,eq,eqs,odd,odds,ts,subst,gs,sol1,tsol,sol2},
design=KroneckerProduct@@Table[{{1,0},{1,1}},{i,1,k}]; (* form design matrix corresponding to conjuction features *)
params=(Subscript[θ, #])&/@Range[(2^k)];
(* equation for i'th feature *)
eq[i_]:=Sum[Transpose[design][[i]][[x]](((Exp[design[[x]].params])/(Exp[design[[x]].params]+1))-g[x]),{x,1,(2^k)}]==0;
eqs=eq[#]&/@Range[1,(2^k)];
(* predicted odds for given x *)
odd[x_]:=((Exp[design[[x]].params])/(Exp[design[[x]].params]+1));
odds=odd[#]&/@Range[(2^k)];
ts=(Subscript[t, #])&/@Range[(2^k)];
subst=Thread[odds->ts];
gs=g[#]&/@Range[(2^k)] ;(* these are empirical probabilities of y=1 for each x *)
sol1=Solve[eqs/.subst,rightvars] ;(* solution for odds in terms of empirical odds, trivial for saturated models *)
tsol=MapThread[Solve,{Thread[odds==ts],params}]//Flatten ;
sol2=tsol[[All,2]]//.tsol;
sol2/.sol1
]

solvelog[2]//Simplify