# How are models combined in locally weighted linear regression?

The cost function in locally weighted linear regression is $$J(\theta)= \frac{1}{2}\sum_{i=1}^mw^{(i)}(y^{(i)} - \theta^Tx^{(i)}) =(X\theta - y)^TW(X\theta-y),$$

where $x^{(i)}$ is the $i$-th instance, $y^{(i)}$ is its corresponding class label, $\theta$ are the model parameters, and $w^{(i)}$ is the weight of the $i$-th instance, given by $$w^{(i)} = \exp\Big(-\frac{(x - x^{(i)})^2}{2\tau^2}\Big),$$

where $\tau$ is the bandwidth, and $x$ is the query point, which is fixed for a given regression model and is typically one of the instances.

I am trying to implement locally weighted linear regression, where my dataset is such that $y\in\mathbb{R}^{450\times 1}$, $X\in\mathbb{R}^{450\times 2}$ (including the intercepts), and $W\in\mathbb{R}^{450\times 450}$. The point of weighted linear regression is to choose a query point, $x$, meaning that the instances close to $x$ are weighted more heavily in the regression model. Following this question: https://datascience.stackexchange.com/questions/16850/understanding-locally-weighted-linear-regression, I believe a model is supposed to be generated for each query point (i.e. a query point corresponding to each instance), meaning that 450 models will be generated in this case. The models are then supposed to be combined to give a final model.

I take the following steps:

thetas = []
for instance in X:
Set current instance as the query point
Compute weights for all instances using the equation above
Compute optimal parameters using the equation for theta above
Append these parameters to thetas


And this gives us 450 linear regression models for the data, with each model being weighted around a specific query point.

(If this is an incorrect implementation, please correct me.)

Main problem

The final stage is to combine these 450 models to give the final weighted linear regression model.

How do I do this?

Am I supposed to take some kind of average of the parameters across each model? My Google searches so far haven't been fruitful.

• Im not following what you're doing here. What's a "query point"? Why do you have 450 different regressions? I think there is some domain specific knowledge here we do not have. – Matthew Drury May 23 '18 at 23:49
• also what are "the weight of the i-th instance", and $x$, $x^{(i)}$, $w^{(i)}$? – Sextus Empiricus May 24 '18 at 0:13
• It seems to me that you only have one regression; $X^T$ is 2x450, $W$ is 450x450, $X$ is 450x2, so $(X^TWX)^{-1}$ will be 2x2, just as you want; $X^TWy$ is 2x1, so $\hat{\theta}$ is a vector of length 2, as it should be, and is computable using the obvious matrix algebra . – jbowman May 24 '18 at 3:43
• @MatthewDrury The query point is the point about which the linear regression is weighted. So when you do a weighted linear regression, you do it with respect to a query point, $x$. With the equation for $w$ I have given, this corresponds to allocating more importance to the instances in the dataset close to the query point. I believe you repeat this for each query point, and then combine the models: datascience.stackexchange.com/questions/16850/… – quanty May 24 '18 at 8:07
• @MartijnWeterings it is from Q5b(ii) in this sheet: cs229.stanford.edu/ps/ps1/ps1.pdf I believe the next stage is smoothing the data, yes. – quanty May 24 '18 at 11:18

What you are doing sounds like filtering or smoothening data. It resembles something like

for which there are many great resources.

If you do something like such smoothening or filtering then you do the algorithm step (fitting a regression line) once for each point (or in your case 450 times).

Every prediction/estimate/filtered-value $\hat{y}^{(i)}$ that you make, requires a separate regression. Thus 450 regressions to predict 450 $\hat{y}^{(i)}$.

In the image below you have 450 measurements (gray points) with 450 different locally weighted (using $\tau=5$) linear curves (if you like other orders could be used) and points fitted to it (blue/red).

For clarity only 11 of the fitted curves and points (the ones in red) are made salient.

Note that each predicted point is associated to a single regression curve. So the "model combining in locally weighted linear regression" is done by using a single model for each single point.

Note: effectively this is sort of linear smoothening like the trapezium rule or Simpson's rule (which can still be derived easily) only is it more flexible with uneven spaced $x$ values and the kernel may adapt from place to place.

### Combining the regressions

The combination of regressions is eventually done in a different way (in part ii of the question in your lecture notes). The 200 smoothened training samples are used to make predictions of the left side of the spectrum based on the right side of the spectrum (this is useful when you train the model on a selection of quasars for which you observe both sides very well and then wish to apply the trained-and tested model on more difficult quasars which have the left side obstructed).

So the "combination" is done by predicting the left side of the spectrum as a weighted function of the 200 training samples using weights defined by a distance between samples on the right side of the spectrum.

What is not explicitly done in this example is 'training' of some fitting parameters (I am not sure what the fitting parameter is maybe $h$, but it is not a linear regression which is performed). But, anyway this is how the 'models' are 'combined' in this example.

Martijn Weterings posts a valid answer to this question, but I am going to add a bit more to this. With respect to the pseudo-algorithm written in the question, it isn't quite complete. It should be like this:

models = []
for instance in X:
Set current instance as the query point
Compute weights for all instances using the equation above
Compute optimal parameters using the equation for theta above
Compute y_pred by taking the dot product of the optimal thetas with the instance
Append y_pred to models


The values of $\hat{y}$ inside the models list can then be plot against $X$ to see the smoothing.

I have done this on the first instance of the training data set that can be found on question 5 at this link: http://cs229.stanford.edu/ps/ps1/ps1.pdf

This is the result of a simple linear regression:

After the smoothing using weighted linear regression described above, the output is:

• Could you provide data of your sample? – Sextus Empiricus May 24 '18 at 16:06
• @MartijnWeterings edited to include the link to the data. If you are planning on re-implementing, note that this is for only the first instance of the dataset. (If that doesn't make sense, it will when you look at the data.) – quanty May 24 '18 at 16:10

As I am also working on the same problem as you (CS229!) I have been looking into implementing locally weighted regression in python. The following implementation is heavily based on the version provided by Alexandre Gramfort, an Sklearn developper, on his github page, but uses a bell shaped kernel as you have described in your question.

def lowess_bell_shape_kern(x, y, tau = .005):
"""lowess_bell_shape_kern(x, y, tau = .005) -> yest
Locally weighted regression: fits a nonparametric regression curve to a
scatterplot.
The arrays x and y contain an equal number of elements; each pair
(x[i], y[i]) defines a data point in the scatterplot. The function returns
the estimated (smooth) values of y.
The kernel function is the bell shaped function with parameter tau. Larger tau
will result in a
smoother curve.
"""
m = len(x)
yest = np.zeros(n)

#Initializing all weights from the bell shape kernel function
w = np.array([np.exp(- (x - x[i])**2/(2*tau)) for i in range(m)])

#Looping through all x-points
for i in range(n):
weights = w[:, i]
b = np.array([np.sum(weights * y), np.sum(weights * y * x)])
A = np.array([[np.sum(weights), np.sum(weights * x)],
[np.sum(weights * x), np.sum(weights * x * x)]])
theta = linalg.solve(A, b)
yest[i] = theta[0] + theta[1] * x[i]

return yest


The trick is to express the problem in matrix form and to solve the system of equations using np.linalg, for each point in the data set.

### Vectorized implementation

For those who might be interested here is my attempt at describing the vectorized implementation mathematically.

Consider the 1D case where $\Theta = [\theta_0, \theta_1]$ and $x$ and $y$ are vectors of size $m$. The cost function $J(\theta)$ is a weighted version of the OLS regression, where the weights $w$ are defined by some kernel function

\begin{aligned} J(\theta) &= \sum_{i=1}^m w^{(i)} \left( y^{(i)} - (\theta_0 + \theta_1 x^{(i)}) \right)^2 \\ \frac{\partial J}{\partial \theta_0} &= -2 \sum_{i=1}^m w^{(i)} \left( y^{(i)} - (\theta_0 + \theta_1 x^{(i)}) \right) \\ \frac{\partial J}{\partial \theta_1} &= -2 \sum_{i=1}^m w^{(i)} \left( y^{(i)} - (\theta_0 + \theta_1 x^{(i)}) \right) x^{(i)} \end{aligned}

Cancelling the $-2$ terms, equating to zero, expanding and re-arranging the terms: \begin{aligned} & \frac{\partial J}{\partial \theta_0} = \sum_{i=1}^m w^{(i)} \left( y^{(i)} - (\theta_0 + \theta_1 x^{(i)}) \right) = 0 \\ & \sum_{i=1}^m w^{(i)} \theta_0 + \sum_{i=1}^m w^{(i)} \theta_1 x^{(i)} = \sum_{i=1}^m w^{(i)} y^{(i)} &\text{Eq. (1)} \\ \\ & \frac{\partial J}{\partial \theta_1} = \sum_{i=1}^m w^{(i)} \left( y^{(i)} - (\theta_0 + \theta_1 x^{(i)}) \right) x^{(i)} = 0 \\ & \sum_{i=1}^m w^{(i)} \theta_0 + \sum_{i=1}^m w^{(i)} \theta_1 x^{(i)} x^{(i)} = \sum_{i=1}^m w^{(i)} y^{(i)} x^{(i)} &\text{Eq. (2)} \end{aligned}

Writing Eq. (1) and Eq. (2) in matrix form $\mathbf{A \Theta = b}$ allows us to solve for $\Theta$ \begin{aligned} & \sum_{i=1}^m w^{(i)} \theta_0 + \sum_{i=1}^m w^{(i)} \theta_1 x^{(i)} = \sum_{i=1}^m w^{(i)} y^{(i)} \\ & \sum_{i=1}^m w^{(i)} \theta_0 + \sum_{i=1}^m w^{(i)} \theta_1 x^{(i)} x^{(i)} = \sum_{i=1}^m w^{(i)} y^{(i)} x^{(i)} \\ & \begin{bmatrix} \sum w^{(i)} & \sum w^{(i)} x^{(i)} \\ \sum w^{(i)} x^{(i)} & \sum w^{(i)} x^{(i)} x^{(i)} \end{bmatrix} \begin{bmatrix} \theta_0 \\ \theta_1 \end{bmatrix} = \begin{bmatrix} \sum w^{(i)} y^{(i)} \\ \sum w^{(i)} y^{(i)} x^{(i)} \end{bmatrix} \\ & \mathbf{A} \Theta = \mathbf{b} \\ & \Theta = \mathbf{A}^{-1} \mathbf{b} \end{aligned}