Weighted principal components analysis

Question

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?

Answer 1

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.

Answer 2

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):

Cheers,
Andres