# Locally Weighted Linear Regression

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?

• Why are you using Minimize[]? Do you really need symbolic results? If not, use NMinimize[] instead. – J. M. will be back soon Jun 8 '16 at 11:33
• @J.M. This is essentially an answer. NMinimize gives immediate results. – Feyre Jun 8 '16 at 12:14
• @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! – Frederik Brinck Jensen Jun 8 '16 at 23:09

## 3 Answers

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"]


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}]


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}]