8
$\begingroup$

I am following the online Stanford Course in Machine Learning and I am trying to implement the algorithms on the fly in Mathematica.

Currently, I have issues with minimizing the Locally Weighted Linear Regression -- it hangs and just runs for more than 3 hours. I was hoping you might have some clues as to why.

The code I have is the following

trainX = {100, 320, 213, 512, 58, 84, 113, 142, 93, 121, 421, 432, 
   249, 254};
trainY = {140000, 400000, 241000, 489000, 78000, 123000, 139000, 
   143000, 97000, 134000, 392000, 458000, 311000, 378000};
trainComposed = 
  Transpose@{trainY, Table[1, {Length[trainX]}], trainX};

I define my hypothesis as

Hyp[x_, T_] := T.x;

and the cost function as

Cost[T_, case_, t_] := 
     1/2 Total[Exp[-(({##2} - case).({##2} - case)/(2t^2))]
               (Hyp[{##2}, T] - #1)^2 & @@@ trainComposed]

where T is a list of the parameters we are trying to find, case is the particular values on x for which we want to fit the model, and t is an arbitrary parameter that defines the weight of neighbouring data points.

I can calculate the cost without any issues, but when I try to minimise, it hangs. The code I use to run it with is

Clear[v]
T = Array[v, Length[trainComposed[[1]]] - 1];
Minimize[Cost[T, {1, 121}, 50], T]

I am essentially trying to minimise the following function

$$\sum_iw^{(i)}(y^{(i)}-\theta^Tx^{(i)})^2$$ where $$w^{(i)}=\exp\left(-\frac{(x^{(i)}-x)^2}{2\tau^2}\right)$$

and $y^{(i)}$ corresponds to the target value for the $y^{th}$ training example, $x^{(i)}$ is a feature vector of the $i^\text{th}$ example, and $x$ is the position vector for which we are currently doing the fit.

What am I doing wrong?

$\endgroup$
  • 3
    $\begingroup$ Why are you using Minimize[]? Do you really need symbolic results? If not, use NMinimize[] instead. $\endgroup$ – J. M. will be back soon Jun 8 '16 at 11:33
  • $\begingroup$ @J.M. This is essentially an answer. NMinimize gives immediate results. $\endgroup$ – Feyre Jun 8 '16 at 12:14
  • $\begingroup$ @J.M. it works, thanks. Can you please submit an answer? Moreover, I'd be happy if you could elaborate a bit on why Minimize wouldn't work in this example; I have used it earlier without any issues to a similar problem, and thus I didn't think of this distinction. Thanks! $\endgroup$ – Frederik Brinck Jensen Jun 8 '16 at 23:09
9
$\begingroup$

In part, your trouble stemmed from an attempt to determine exact/symbolic results from your data. Even for this modestly-sized set, the objective function you were feeding to Minimize[] is already sufficiently complicated, which is why it takes long.

As I previously noted, using NMinimize[] instead will give you the numerical values you need. Less expensively, tho, one can use the built-in functions specialized for the purpose. For instance, you can use LeastSquares[] like so:

With[{case = {1, 121}, t = 50}, 
     wd = DiagonalMatrix[Exp[-Map[Composition[#.# &, {1, #} - case &], trainX]/(4 t^2)]]; 
     LeastSquares[wd.DesignMatrix[Transpose[{trainX, trainY}], x, x], wd.(N @ trainY)]]

The sufficiently observant reader will notice two features in this weighted least squares problem: 1. the diagonal matrix contains the square roots of the original weights (think about why this must be so); and 2. we need to use N[] in this case as well (without it, we end up doing the same thing that doomed the Minimize[] approach).

Of course, since weighted regression is a relatively common operation, Mathematica provides for a function called, appropriately enough, LinearModelFit[]. Here's how to use it for your problem:

With[{case = {1, 121}, t = 50}, 
     LinearModelFit[Transpose[{trainX, trainY}], x, x, 
                    Weights -> Exp[-Map[Composition[#.# &, {1, #} - case &],
                                        trainX]/(2 t^2)]] @ "BestFitParameters"]
$\endgroup$
3
$\begingroup$

Just a follow-up to @J.M.'s answer to show the effect of the value of t:

trainX = {100, 320, 213, 512, 58, 84, 113, 142, 93, 121, 421, 432, 249, 254};
trainY = {140000, 400000, 241000, 489000, 78000, 123000, 139000, 
   143000, 97000, 134000, 392000, 458000, 311000, 378000};
rX = MinMax[trainX];
rY = MinMax[trainY];

Manipulate[

 (* Make a table of predictions across the range of the predictor variable *)
 predicted = Table[LinearModelFit[Transpose[{trainX, trainY}], x, x,
     Weights -> Exp[-Map[Composition[#.# &, {1, #} - {1, z} &], trainX]/(2 t^2)]]@z,
   {z, rX[[1]], rX[[2]], (rX[[2]] - rX[[1]])/100}];

 (* Plot results *)
 ListPlot[{
   Transpose[{trainX, trainY}],
   Transpose[{Range[rX[[1]], rX[[2]], (rX[[2]] - rX[[1]])/100], 
     predicted}]},
  PlotRange -> {rX, rY}, Joined -> {False, True}],

 {{t, 50}, 5, 300, Appearance -> "Labeled"},
 TrackedSymbols :> {t}]

Data and fit Data and fit with different tuning parameter

$\endgroup$
3
$\begingroup$

Another follow up, specifically to @Jim Baldwin's comment. Instead of plotting the full curve, we can also plot a certain number of points and see how they vary with the parameter t, i.e.:

Clear[v]
th = Array[v, Length[trainComposed[[1]]] - 1];
lwrPlot = 
 Manipulate[
  Show[listPlot, 
   With[{x = #}, 
      With[{bounds = 
         Values@NMinimize[
           Cost[th, {1, x}, 
             t], th][[2]]}, {Plot[
         Hyp[{1, y}, bounds], {y, x - span, x + span}, 
         PlotStyle -> {Orange, Thin}], 
        ListPlot[{{x, Hyp[{1, x}, bounds]}}, 
         PlotStyle -> Orange]}]] & /@ 
    Range[Min[trainX], Max[trainX], (Max[trainX] - Min[trainX])/(
     points - 1)]], {{t, 20}, 5, 80, 5}, {{span, 40}, 10, 200, 
   10}, {{points, 5}, 2, 10, 1}]

Example of executed code

$\endgroup$

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.