How to implement ZCA Whitening? Python

Question

Im trying to implement ZCA whitening and found some articles to do it, but they are a bit confusing.. can someone shine a light for me?

Any tip or help is appreciated!

Here is the articles i read :

http://courses.media.mit.edu/2010fall/mas622j/whiten.pdf http://bbabenko.tumblr.com/post/86756017649/learning-low-level-vision-feautres-in-10-lines-of

I tried several things but most of them i didnt understand and i got locked at some step. Right now i have this as base to start again :

dtype = np.float32
data = np.loadtxt("../inputData/train.csv", dtype=dtype, delimiter=',', skiprows=1)
img = ((data[1,1:]).reshape((28,28)).astype('uint8')*255)

Answer 1

Here is a python function for generating the ZCA whitening matrix:

def zca_whitening_matrix(X):
    """
    Function to compute ZCA whitening matrix (aka Mahalanobis whitening).
    INPUT:  X: [M x N] matrix.
        Rows: Variables
        Columns: Observations
    OUTPUT: ZCAMatrix: [M x M] matrix
    """
    # Covariance matrix [column-wise variables]: Sigma = (X-mu)' * (X-mu) / N
    sigma = np.cov(X, rowvar=True) # [M x M]
    # Singular Value Decomposition. X = U * np.diag(S) * V
    U,S,V = np.linalg.svd(sigma)
        # U: [M x M] eigenvectors of sigma.
        # S: [M x 1] eigenvalues of sigma.
        # V: [M x M] transpose of U
    # Whitening constant: prevents division by zero
    epsilon = 1e-5
    # ZCA Whitening matrix: U * Lambda * U'
    ZCAMatrix = np.dot(U, np.dot(np.diag(1.0/np.sqrt(S + epsilon)), U.T)) # [M x M]
    return ZCAMatrix

And an example of the usage:

X = np.array([[0, 2, 2], [1, 1, 0], [2, 0, 1], [1, 3, 5], [10, 10, 10] ]) # Input: X [5 x 3] matrix
ZCAMatrix = zca_whitening_matrix(X) # get ZCAMatrix
ZCAMatrix # [5 x 5] matrix
xZCAMatrix = np.dot(ZCAMatrix, X) # project X onto the ZCAMatrix
xZCAMatrix # [5 x 3] matrix

Hope it helps!

Details for why Edgar Andrés Margffoy Tuay's answer is not correct: As pointed out in R.M's comment, Edgar Andrés Margffoy Tuay's ZCA whitening function contains a small, but crucial mistake: the np.diag(S) should be removed. Numpy returns S as a m x 1 vector and not a m x m matrix (as is common to other svd implementations, e.g. Matlab). Hence the ZCAMatrix variable becomes a m x 1 vector and not a m x m matrix as it should be (when the input is m x n). (Also, the covariance matrix in Andfoy's answer is only valid if X is pre-centered, i.e mean 0).

Other references for ZCA: You can see the full answer, in Python, to the Stanford UFLDL ZCA Whitening exercise here.

Answer 2

Is your data stored in an mxn matrix? Where m is the dimension of the data and n are the total number of cases? If that's not the case, you should resize your data. For instance if your images are of size 28x28 and you have only one image, you should have a 1x784 vector. You could use this function:

import numpy as np

def flatten_matrix(matrix):
    vector = matrix.flatten(1)
    vector = vector.reshape(1, len(vector))
    return vector

Then you apply ZCA Whitening to your training set using:

def zca_whitening(inputs):
    sigma = np.dot(inputs, inputs.T)/inputs.shape[1] #Correlation matrix
    U,S,V = np.linalg.svd(sigma) #Singular Value Decomposition
    epsilon = 0.1                #Whitening constant, it prevents division by zero
    ZCAMatrix = np.dot(np.dot(U, np.diag(1.0/np.sqrt(np.diag(S) + epsilon))), U.T)                     #ZCA Whitening matrix
    return np.dot(ZCAMatrix, inputs)   #Data whitening

It is important to save the ZCAMatrix matrix, you should multiply your test cases if you want to predict after training the Neural Net.

Finally, I invite you to take the Stanford UFLDL Tutorials at http://ufldl.stanford.edu/wiki/index.php/UFLDL_Tutorial or http://ufldl.stanford.edu/tutorial/ . They have pretty good explanations and also some programming exercises on MATLAB, however, almost all the functions found on MATLAB are on Numpy by the same name. I hope this may give an insight.

Answer 3

I may be a little late to the discussion, but I found this thread recently as I struggled to implement ZCA in TensorFlow because my poor PC processor was too slow to process large volume of data.

If anyone is interested, I have made a gist of my implementation of the ZCA in TensorFlow:

import tensorflow as tf

from keras.datasets import mnist

import numpy as np

tf.enable_eager_execution()

assert tf.executing_eagerly()


class ZCA(object):
    """
    Simple ZCA aka Mahalanobis transformation class made in TensorFlow.
    The code was largely ported from Keras ImageDataGenerator
    """

    def __init__(self, epsilon=1e-5, dtype='float64'):

        """epsilon is the normalization constant, dtype refers to the data type used in the computation.
         WARNING: the default precision is set to float64 as i have found that when computing the mean tensorflow'
         and numpy results can differ by a substantial amount.
         Usage: fit method computes the principal components and should be called first,
                compute method returns the actual transformed tensor
         NOTE : The input to both methods must be a 4D tensor.
        """

        assert dtype is 'float32' or 'float64', "precision must be float32 or float64"

        self.epsilon = epsilon
        self.dtype = dtype
        self.princ_comp = None
        self.mean = None

    def _featurewise_center(self, images_tensor):

        if self.mean is None:
            self.mean, _ = tf.nn.moments(images_tensor, axes=(0, 1, 2))
            broadcast_shape = [1, 1, 1]
            broadcast_shape[2] = images_tensor.shape[3]
            self.mean = tf.reshape(self.mean, broadcast_shape)

        norm_images = tf.subtract(images_tensor, self.mean)

        return norm_images

    def fit(self, images_tensor):

        assert images_tensor.shape[3], "The input should be a 4D tensor"

        if images_tensor.dtype is not self.dtype:  # numerical error for float32

            images_tensor = tf.cast(images_tensor, self.dtype)

        images_tensor = self._featurewise_center(images_tensor)

        flat = tf.reshape(images_tensor, (-1, np.prod(images_tensor.shape[1:].as_list())))
        sigma = tf.div(tf.matmul(tf.transpose(flat), flat), tf.cast(flat.shape[0], self.dtype))
        s, u, _ = tf.svd(sigma)
        s_inv = tf.div(tf.cast(1, self.dtype), (tf.sqrt(tf.add(s[tf.newaxis], self.epsilon))))
        self.princ_comp = tf.matmul(tf.multiply(u, s_inv), tf.transpose(u))

    def compute(self, images_tensor):

        assert images_tensor.shape[3], "The input should be a 4D tensor"

        assert self.princ_comp is not None, "Fit method should be called first"

        if images_tensor.dtype is not self.dtype:
            images_tensor = tf.cast(images_tensor, self.dtype)

        images_tensors = self._featurewise_center(images_tensor)

        flatx = tf.cast(tf.reshape(images_tensors, (-1, np.prod(images_tensors.shape[1:]))), self.dtype)
        whitex = tf.matmul(flatx, self.princ_comp)
        x = tf.reshape(whitex, images_tensors.shape)

        return x


def main():
    import matplotlib.pyplot as plt

    train_set, test_set = mnist.load_data()
    x_train, y_train = train_set

    zca1 = ZCA(epsilon=1e-5, dtype='float64')

    # input should be a 4D tensor

    x_train = x_train.reshape(*x_train.shape, 1)
    zca1.fit(x_train)
    x_train_transf = zca1.compute(x_train)

    # reshaping to 28*28 and casting to uint8 for plotting

    x_train_transf = tf.reshape(x_train_transf, x_train_transf.shape[0:3])


    fig, axes = plt.subplots(3, 3)

    for i, ax in enumerate(axes.flat):
        # Plot image.
        ax.imshow(x_train_transf[i],
                  cmap='binary'
                  )

        xlabel = "True: %d" % y_train[i]
        ax.set_xlabel(xlabel)
        ax.set_xticks([])
        ax.set_yticks([])

    plt.show()


if __name__ == '__main__':
main()

I know this isn't a proper answer to the original question, but still it may be useful to anyone who is looking for a GPU implementation of ZCA but couldn't find one.

Answer 4

Although both answers refer to the UFLDL tutorial, none of them seems to use the steps as described in it.

Therefore, I thought it might not be bad idea to just provide an answer that simply implements PCA/ZCA-whitening according to the tutorial:

import numpy as np

# generate some random, 2D data
x = np.random.randn(1000, 2)
# and center it
x_c = x - np.mean(x, 0)

# compute the 2x2 covariance matrix
# (remember that covariance matrix is symmetric)
sigma = np.cov(x, rowvar=False)
# and extract eigenvalues and eigenvectors
# using the algorithm for symmetric matrices
l,u = np.linalg.eigh(sigma)
# NOTE that for symmetric matrices,
# eigenvalues and singular values are the same.
# u, l, _ = np.linalg.svd(sigma) should thus give equivalent results

# rotate the (centered) data to decorrelate it
x_rot = np.dot(x_c, u)
# check that the covariance is diagonal (indicating decorrelation)
np.allclose(np.cov(x_rot.T), np.diag(np.diag(np.cov(x_rot.T))))

# scale the data by eigenvalues to get unit variance
x_white = x_rot / np.sqrt(l)
# have the whitened data be closer to the original data
x_zca = np.dot(x_white, u.T)

I assume you can wrap this in a function by yourself...

For completeness, different implementation flavours and their runtime (evaluated on a centred version of CIFAR10):

x = np.random.randn(10_000, 3, 32, 32)
x_ = np.reshape(x, (len(x), -1))
x_c = x_ - np.mean(x_, axis=0)

def zca1(x):
    s, u = np.linalg.eigh(x.T @ x)
    scale = np.sqrt(len(x) / s)
    return (u * scale) @ u.T

def zca2(x):
    u, s, _ = np.linalg.svd(x.T @ x, hermitian=True)
    scale = np.sqrt(len(x) / s)
    return (u * scale) @ u.T

def zca3(x):
    _, s, v = np.linalg.svd(x, full_matrices=False)
    scale = np.sqrt(len(x)) / s
    return (v.T * scale) @ v


%timeit zca1(x_c)  
# 4.57 s ± 14.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit zca2(x_c)
# 4.62 s ± 22.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit zca3(x_c)
# 20.2 s ± 1.2 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

For the mathematics behind this, I refer to this excellent answer from cross validated.

Answer 5

This works with an array of 48x48:

def flatten_matrix(matrix):
    vector = matrix.flatten(order='F')
    vector = vector.reshape(1, len(vector))
    return vector

def zca_whitening(inputs): 
    sigma = np.dot(inputs, inputs.T)/inputs.shape[1] #Correlation matrix
    U,S,V = np.linalg.svd(sigma) #Singular Value Decomposition
    epsilon = 0.1                #Whitening constant, it prevents division by zero
    ZCAMatrix = np.dot(np.dot(U, np.diag(1.0/np.sqrt(np.diag(S) + epsilon))), U.T)  #ZCA Whitening matrix
    return np.dot(ZCAMatrix, inputs)   #Data whitening

def global_contrast_normalize(X, scale=1., subtract_mean=True, use_std=True,
                              sqrt_bias=10, min_divisor=1e-8):

    """
    __author__ = "David Warde-Farley"
    __copyright__ = "Copyright 2012, Universite de Montreal"
    __credits__ = ["David Warde-Farley"]
    __license__ = "3-clause BSD"
    __email__ = "wardefar@iro"
    __maintainer__ = "David Warde-Farley"
    .. [1] A. Coates, H. Lee and A. Ng. "An Analysis of Single-Layer
       Networks in Unsupervised Feature Learning". AISTATS 14, 2011.
       http://www.stanford.edu/~acoates/papers/coatesleeng_aistats_2011.pdf
    """
    assert X.ndim == 2, "X.ndim must be 2"
    scale = float(scale)
    assert scale >= min_divisor

    mean = X.mean(axis=1)
    if subtract_mean:
        X = X - mean[:, np.newaxis]  
    else:
        X = X.copy()
    if use_std:
        ddof = 1
        if X.shape[1] == 1:
            ddof = 0
        normalizers = np.sqrt(sqrt_bias + X.var(axis=1, ddof=ddof)) / scale
    else:
        normalizers = np.sqrt(sqrt_bias + (X ** 2).sum(axis=1)) / scale
    normalizers[normalizers < min_divisor] = 1.
    X /= normalizers[:, np.newaxis]  # Does not make a copy.
    return X

def ZeroCenter(data):
    data = data - np.mean(data,axis=0)
    return data

def Zerocenter_ZCA_whitening_Global_Contrast_Normalize(data):
    numpy_data = np.array(data).reshape(48,48)
    data2 = ZeroCenter(numpy_data)
    data3 = zca_whitening(flatten_matrix(data2)).reshape(48,48)
    data4 = global_contrast_normalize(data3)
    data5 = np.rot90(data4,3)
    return data5

for example from this image:

returns:

Here is the code:

https://gist.github.com/m-alcu/45f4a083cb5e388d2ed26ace4392ed66, needs to put fer2013.csv file in the same directory (https://www.kaggle.com/c/challenges-in-representation-learning-facial-expression-recognition-challenge/data)