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?

  • 1
    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. – ttnphns Oct 16 '16 at 9:35
up vote 27 down vote accepted

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. – noname Aug 27 '14 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.) – amoeba Aug 27 '14 at 21:44
  • 3
    @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. – amoeba Aug 27 '14 at 21:44
  • 2
    +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). – ttnphns Oct 6 '15 at 10:07
  • 1
    @JánЯabčan Thanks, it was a mistake. I fixed it! – amoeba Feb 6 '17 at 19:52

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
    p_grad = weighted_pca_regression(data_x, data_y, weights)
    # plot for illustration purposes
    line_x = np.linspace(-5,5,10)
    line_y = np.tan(p_grad)*line_x
    plt.plot(line_x+center[0], line_y+center[1], c="r")
    plt.scatter(*data_xy)
    plt.show()

And a sample output (it does the same for every dot): enter image description here

Cheers,
Andres

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.