# Weighted principal components analysis

After some searching, I find very little on the incorporation of observation weights/measurement errors into principal components analysis. What I do find tends to rely on iterative approaches to include weightings (e.g., here). My question is why is this approach necessary? Why can't we use the eigenvectors of the weighted covariance matrix?

• In addition to the answer(s) below, please see thread stats.stackexchange.com/q/141754/3277, where weighted PCA (with weights on columns and/or rows) is explained as primarily equivalent to weighted (generalized) svd/biplot. Oct 16, 2016 at 9:35

It depends on what exactly your weights apply to.

### Row weights

Let $$\mathbf{X}$$ be the data matrix with variables in columns and $$n$$ observations $$\mathbf x_i$$ in rows. If each observation has an associated weight $$w_i$$, then it is indeed straightforward to incorporate these weights into PCA.

First, one needs to compute the weighted mean $$\boldsymbol \mu = \frac{1}{\sum w_i}\sum w_i \mathbf x_i$$ and subtract it from the data in order to center it.

Then we compute the weighted covariance matrix $$\frac{1}{\sum w_i}\mathbf X^\top \mathbf W \mathbf X$$, where $$\mathbf W = \operatorname{diag}(w_i)$$ is the diagonal matrix of weights, and apply standard PCA to analyze it.

### Cell weights

The paper by Tamuz et al., 2013, that you found, considers a more complicated case when different weights $$w_{ij}$$ are applied to each element of the data matrix. Then indeed there is no analytical solution and one has to use an iterative method. Note that, as acknowledged by the authors, they reinvented the wheel, as such general weights have certainly been considered before, e.g. in Gabriel and Zamir, 1979, Lower Rank Approximation of Matrices by Least Squares With Any Choice of Weights. This was also discussed here.

As an additional remark: if the weights $$w_{ij}$$ vary with both variables and observations, but are symmetric, so that $$w_{ij}=w_{ji}$$, then analytic solution is possible again, see Koren and Carmel, 2004, Robust Linear Dimensionality Reduction.

• Thank you for the clarification. Can you explain why no analytic solution is possible with off-diagonal weights? I this is what I am missing from both Tamuz et al 2013 and Gabriel and Zamir 1979. Aug 27, 2014 at 21:00
• @noname: I am not aware of such a proof, and moreover I would not be surprised if it were not known. It is generally quite tricky to prove that something is not possible, in particular that something is not possible analytically. The impossibility of angle trisection famously waited for its proof for over 2000 years... (cont.) Aug 27, 2014 at 21:44
• @noname: (cont.) What you are asking is to show that the problem of minimizing $\sum_{i,j} w_{ij}(X_{ij} - A_{ij})^2$ with respect to $A$ constrained to have low rank $q$. is not reducible to an eigenvector problem. I am afraid you would need another forum for that (maybe mathoverflow?). But note that finding eigenvectors is also not exactly an analytical solution: it's just that the iterations are usually silently performed by a standard library function. Aug 27, 2014 at 21:44
• +1. The first section of the answer can be conceptualized also in terms of Weighted (Generalized) Biplot as described here. Keeping in mind how PCA is a "specific case of" Biplot (also concerned in the lined answer). Oct 6, 2015 at 10:07
• @ttnphns: After your comment and another thread being closed as a duplicate, I re-read my answer and expanded the explanation of how to deal with row weights. I think previously it was not completely correct or at least was not complete because I did not mention the centering with a weighted mean. I hope it makes more sense now! Oct 6, 2015 at 20:44

Thank you very much amoeba for the insight regarding row weights. I know that this is not stackoverflow, but I had some difficulties to find an implementation of row-weighted PCA with explanation and, since this is one of the first results when googling for weighted PCA, I thought it would be good to attach my solution, maybe it can help others in the same situation. In this Python2 code snippet, a PCA weighted with an RBF kernel as the one described above is used to calculate the tangents of a 2D dataset. I will be very happy to hear some feedback!

def weighted_pca_regression(x_vec, y_vec, weights):
"""
Given three real-valued vectors of same length, corresponding to the coordinates
and weight of a 2-dimensional dataset, this function outputs the angle in radians
of the line that aligns with the (weighted) average and main linear component of
the data. For that, first a weighted mean and covariance matrix are computed.
Then u,e,v=svd(cov) is performed, and u * f(x)=0 is solved.
"""
input_mat = np.stack([x_vec, y_vec])
weights_sum = weights.sum()
# Subtract (weighted) mean and compute (weighted) covariance matrix:
mean_x, mean_y =  weights.dot(x_vec)/weights_sum, weights.dot(y_vec)/weights_sum
centered_x, centered_y = x_vec-mean_x, y_vec-mean_y
matrix_centered = np.stack([centered_x, centered_y])
weighted_cov = matrix_centered.dot(np.diag(weights).dot(matrix_centered.T)) / weights_sum
# We know that v rotates the data's main component onto the y=0 axis, and
# that u rotates it back. Solving u.dot([x,0])=[x*u[0,0], x*u[1,0]] gives
# f(x)=(u[1,0]/u[0,0])x as the reconstructed function.
u,e,v = np.linalg.svd(weighted_cov)
return np.arctan2(u[1,0], u[0,0]) # arctan more stable than dividing

# USAGE EXAMPLE:
# Define the kernel and make an ellipse to perform regression on:
rbf = lambda vec, stddev: np.exp(-0.5*np.power(vec/stddev, 2))
x_span = np.linspace(0, 2*np.pi, 31)+0.1
data_x = np.cos(x_span)[:-1]*20-1000
data_y = np.sin(x_span)[:-1]*10+5000
data_xy = np.stack([data_x, data_y])
stddev = 1 # a stddev of 1 in this context is highly local
for center in data_xy.T:
# weight the  points based on their euclidean distance to the current center
euclidean_distances = np.linalg.norm(data_xy.T-center, axis=1)
weights = rbf(euclidean_distances, stddev)
# get the angle for the regression in radians
# plot for illustration purposes
line_x = np.linspace(-5,5,10)
plt.plot(line_x+center, line_y+center, c="r")
plt.scatter(*data_xy)
plt.show()


And a sample output (it does the same for every dot): Cheers,
Andres

• Now, since some time has past: in a modern Python 3 setup one would probably use @ instead of .dot(). Also, you should make this general for any kind of input matrix - in terms of number of features. Also, your function should not be called "regression" - as it is no regression. A regression's cost function relates to output error, while a PCA's cost function relates to error orthogonal to input. (This was not a rigorous discription. ;-)) So a PCA is basically a (linear) manifold learner instead of a regression. Jan 23, 2021 at 10:49
• Thanks @Make42 for the feedback! Indeed Python 2 is now very outdated, and a PCA is not a regression (which in this example fails to fit the tangents close to the "vertical" asymptote). Hopefully there are modern libraries by now that take proper care of this, do you know any? otherwise feel free to edit (I can give it a go in a few weeks) Jan 23, 2021 at 11:17
• Regarding the inverse, naybe the fact that eigh yields u=v, orthogonal, helps (inv=transposed) Jan 23, 2021 at 12:19
• My f(x) in the general case is basically the strongest eigenvector (first col of u). The arctan is just a robust division to get the gradient. In the general case, the gradient is a vector itself Jan 23, 2021 at 13:44
• Actually... maybe github.com/jakevdp/wpca is a valid alternative! Jan 23, 2021 at 15:03