# Softmax Regression Classification Algorithm

I'm following the Stanford Machine Learning course while simultaneously developing my Mathematica abilities. I'm having issues with my Softmax Classifier implementation (for ref, have a look at these notes, the last four pages).

Using my own dataset, when performing softmax classification on it, I get the following plot. As seen, I believe the hand-drawn red lines is a better fit, so how come the algorithm doesn't converge to that? I have gone over my code several times, but I can't find what I am doing wrong.

I define the hypothesis function

Hyp[x_, θ_, l_] := Exp[θ[[l]].x] / Sum[Exp[θ[[j]].x], {j, 1, Length[θ]}]


and the log likelihood function

l[θ_] :=
Sum[Log[Product[Hyp[multiX[[i]], θ, l]^Boole[multiY[[i]] == l],
{l, Length[Tally[multiY]]}]], {i, Length[multiX]}]


where multiX and multiY constitute my data in the format

multiX = {{10.5, 3}, {11, 5}, ...}
multiY = {1, 3, ...}


with labels 1 through 3.

I then define the Hessian matrix and the first derivate, such that I can perform Newton's Method

$$\theta:=\theta - H^{-1}\nabla_{\theta}l(\theta)$$

(* Hessian *)
Hessian[θ_] :=
Partition[
Derivative[##][l]@θ & /@
Partition[
Partition[
Flatten[(Table[#, {6}] + IdentityMatrix & /@
IdentityMatrix)], 2], 3], 6]

(* First Order Derivative *)
d[θ_] :=
Derivative[##][l]@θ & /@ Partition[Partition[Flatten[IdentityMatrix], 2], 3]


followed by the update rule

Step[θ_] := Partition[(Flatten[θ] - PseudoInverse[Hessian[θ]].d[θ]), 2]


Lastly, to obtain the optimal parameters, I step through the function

optimal = Nest[Step, {{1, 1}, {1, 1}, {1, 1}}, 100]


whose parameters I use for the hypothesis function, when I subsequently plot it.

Any clue as to why I don't get the optimal fit? I looked into this duplicate, but I couldn't make sense of it.

The data I used is the following:

multiX = {{10.5, 3}, {11, 5}, {8.4, 1}, {7.9, 3}, {9.1, 4}, {9.8, 2}, {8.7, 4.5}, {9.5, 8}, {10.5, 6}, {13.3, 11.2}, {12.7, 9.5}, {11.51, 6.8}, {12, 3.2}, {11.31, 3.24}, {14, 2.1}, {12.51, 4.41}, {13.6, 6.4}}
multiY = {1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3}

• You can use {min, sol} = NMaximize[l[{{x1, x2}, {y1, y2}, {z1, z2}}], {x1, x2, y1, y2, z1, z2}] and optimal = {{x1, x2}, {y1, y2}, {z1, z2}} /. sol to confirm that the problem is in your Newton method, not anywhere else. I recommend that you ask a new question about how to implement this version of Newton's method, how to compute the Hessian and so on. This more direct question will probably gain more interest than this one which is disguised as a more general problem many may not want to attempt. – C. E. Jun 22 '16 at 18:16