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.

  • 3
    $\begingroup$ 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. $\endgroup$ May 23, 2018 at 23:49
  • $\begingroup$ also what are "the weight of the i-th instance", and $x$, $x^{(i)}$, $w^{(i)}$? $\endgroup$ May 24, 2018 at 0:13
  • $\begingroup$ 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 . $\endgroup$
    – jbowman
    May 24, 2018 at 3:43
  • $\begingroup$ @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/… $\endgroup$
    – quanty
    May 24, 2018 at 8:07
  • 1
    $\begingroup$ @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. $\endgroup$
    – quanty
    May 24, 2018 at 11:18

3 Answers 3


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

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.

example for locally weighted linear regression

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:

enter image description here

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

enter image description here

  • $\begingroup$ Could you provide data of your sample? $\endgroup$ May 24, 2018 at 16:06
  • $\begingroup$ @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.) $\endgroup$
    – quanty
    May 24, 2018 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 
    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}


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.