def compute_St(data):
    """
    Given a dataset, computes the variance matrix of its features.
    """
    n_datapoints, n_features = data.shape

    # Computing the 'mean image'. A pixel at position (x,y) in this image is the
    # mean of all the pixels at position (x,y) of the images in the dataset.
    # This corresponds to the 'mu' we have seen in the lectures.
    mu = data.mean(axis=0) # apply along the rows for each columns.
    centered_data = data - mu

    # Computing the covariance matrix
    St = (1. / n_datapoints) * (centered_data.T * centered_data)
    return St