How are models combined in locally weighted linear regression?

Question

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.

Answer 1

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

• Savitzky–Golay filtering
• kernel smoothening
• LOcally WEighted/Estimated Scatterplot Smoothening (LOESS/LOWESS)

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)}$.

---

Example using your question 5 from your lecture notes

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.

Answer 2

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:

Answer 3

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}