# These are the standard imports for CS 111
# This list may change as the quarter goes on
import sys
sys.path.append("/Users/youngmin/Documents/CS/Introduction to Computational Science - F21/CS111-2021-fall/Python")
import os
import time
import math
import numpy as np
import numpy.linalg as npla
from scipy import linalg as spla
import scipy.sparse
import scipy.sparse.linalg
from scipy import integrate
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import axes3d
%matplotlib inline
import cs111
# Loading a black-and-white face image for its visualization.
Face = plt.imread( "faces/face000.bmp", format="bmp" )
print( "Face image has shape:", Face.shape )
print( " and type:", Face.dtype ) # Face is a 64x64 matrix of gray scale
print( " Minimum value:", np.min( Face ) ) # values thay may range between 0 and 255.
print( " Maximum value:", np.max( Face ) )
# Let's show the image.
plt.figure( figsize=(3,3) )
plt.title( r"A $64\times 64$ black-and-white face image" )
plt.imshow( Face, cmap="gray" )
Face image has shape: (64, 64) and type: uint8 Minimum value: 0 Maximum value: 182
<matplotlib.image.AxesImage at 0x123c78d90>
We will create a data set $F$ of gray-scale images. Each image will be a row in $F$, and their pixels will be placed in the columns.
The accompanying folder faces/
has a data set of 178 black-and-white face images. Out of these 178 faces, there are a few of them (indicated in face_descriptions.py
) that should be omitted. Each image is of size $64\times 64$ pixels.
To build our face data set matrix $F$,
First, let's define some utility functions. These functions assume that images are gray-scale matrices of uint8
values between 0 and 255, and vectors are double precision so that we can perform arithmetics on them.
# Some global variables.
HEIGHT = WIDTH = 64
################# Useful functions to go from images to vectors and back #################
imageToVector = lambda Bitmap: Bitmap.reshape( -1 ).astype( np.float64 )
vectorToImage = lambda Vector: Vector.clip( 0, 255 ).reshape( HEIGHT, WIDTH ).astype( np.uint8 )
renderVector = lambda Vector: plt.imshow( vectorToImage( Vector ), cmap="gray" )
scaleImageIntensities = lambda A: np.round( (A - np.min( A )) / (np.max( A ) - np.min( A )) * 255 ).astype( np.uint8 )
##########################################################################################
By using these functions we can build the $n$-by-$p$ matrix $F$ with $n = 169$ faces stacked as row vectors with $p = 64^2$ columns. This way, the image vector $\mathbf{f}_i^T$ is at the $i^{th}$ row of $F$:
$$ F = \begin{array}{r} \textrm{face image 0} \\ \textrm{face image 1} \\ \vdots \\ \textrm{face image }n - 1 \\ \end{array}\begin{pmatrix} f_{0,0} & f_{0,1} & \cdots & f_{0,p-1} \\ f_{1,0} & f_{1,1} & \cdots & f_{1,p-1} \\ \vdots & \vdots & \ddots & \vdots \\ f_{n-1,0} & f_{n-1,1} & \cdots & f_{n-1,p-1} \\ \end{pmatrix} \in \mathbb{R}^{169\times 64^2}. $$Let's next load all images in the faces/
directory to build $F$.
import face_descriptions # Includes the file names we are interested in.
F = [] # We append to F first as a list; then we make it into a matrix.
for fileName in sorted( face_descriptions.face_features ):
if fileName not in face_descriptions.image_to_omit:
bitMap = plt.imread( "faces/" + fileName )
F.append( imageToVector( bitMap ) ) # Flatenning the image.
F = np.array( F ) # Now F is a matrix.
print( "Shape of F is:", F.shape )
plt.figure( figsize=(10, 3) )
plt.title( r"Entire matrix $F$ with $169\times 64^2$ pixels in it" )
plt.imshow( F, cmap="gray" )
Shape of F is: (169, 4096)
<matplotlib.image.AxesImage at 0x123d5bf70>
We can compute the average face as $\mathbf{m} = \frac{1}{n}\sum_{i}\mathbf{f}_i$.
m = np.mean( F, axis=0 ) # Notice the axis=0 parameter!
# Plotting resulting mean face.
plt.figure( figsize=(3,3) )
plt.title( r"The mean face" )
renderVector( m )
<matplotlib.image.AxesImage at 0x123ec8eb0>
Let $\mathbf{m}^T$ be the mean row vector of matrix $F$ as defined above. Suppose we stack $\mathbf{m}^T$, yet again, $n$ times in a matrix $M$. Then, the covariance matrix of $F$ is given by
$$C =\textrm{cov}(F) = \frac{1}{n-1}(F - M)^T (F - M)$$If we define $X = (F - M)$, then
$$\textrm{cov}(X) = \frac{1}{n-1}X^T X = \frac{1}{n-1}(F - M)^T (F - M) = \textrm{cov}(F) = C,$$and
$$C = USV^T = VSV^T,$$since $C$ is symmetric positive semidefinite and $U = V$. We call the columns in $V$ the principal components of $F$, and we can use them to express any face $\mathbf{f}$. Suppose $\mathbf{x} = \mathbf{f} - \mathbf{m}$ is a vectorized/flattened face image with all pixel values centered around the corresponding means. We can express $\mathbf{x}$ as
$$\mathbf{x} = V\mathbf{c} = c_0\mathbf{e}_0 + c_1\mathbf{e}_1 + \cdots + c_{p-1}\mathbf{e}_{p-1},$$where $\mathbf{c} = (c_0, c_1, ..., c_{p-1})^T$ is a vector of coefficients, and
$$ V = \left(\begin{array}{c|c|c|c} ~ & ~ & ~ & ~ \\ \mathbf{e}_0 & \mathbf{e}_1 & \cdots & \mathbf{e}_{p-1} \\ ~ & ~ & ~ & ~ \\ \end{array}\right) $$is the $p$-by-$p$ matrix of singular vectors of $\textrm{cov}(X)$.
To find the vector of coefficients/projections $\mathbf{c}$, we use
$$\mathbf{c} = V^T \mathbf{x},$$and so $\mathbf{x} = V\mathbf{c} = V\left(V^T \mathbf{x}\right) = \mathbf{x}$. Furthermore, since $\mathbf{x} = \mathbf{f} - \mathbf{m}$, then
$$\mathbf{x} = (\mathbf{f} - \mathbf{m}) = V\mathbf{c},$$which implies
$$\mathbf{f} = \mathbf{m} + V\mathbf{c}.$$This is our cookbook to fabricate faces.
With PCA, we can perform data compression or low-rank approximation if we use only the first $k$ singular vectors in $V$. If we do so with any face $\mathbf{f}$, then
$$\mathbf{f}^{(k)} = \mathbf{m} + V^{(k)}\mathbf{c} = \mathbf{m} + c_0\mathbf{e}_0 + c_1\mathbf{e}_1 + \cdots + c_{k-1}\mathbf{e}_{k-1}.$$The singular vectors in $V$ are called eigenfaces.
# We already know F and m. Let's compute X: the centered matrix F around the column means.
X = F - m # Notice that m will be broadcasted to fit the dimensions of X.
# Compute the covariance matrix of the the centered data.
C = np.cov( X.T )
print( "Shape of cov(X):", C.shape )
# Perform PCA on C using the SVD.
U, sigma, Vt = spla.svd( C )
Eigenfaces = Vt.T # The Eigenfaces are the eigenvectors (U = V).
# Let's look at the first largest k eigen/singular values.
k = 30
# Plotting eigenfaces and the explained variance for first k singular values.
fig = plt.figure( figsize=(10,7) )
plt.subplot( 231 )
plt.title( "Singular values" )
plt.plot( sigma[:k], "." )
plt.subplot( 232 )
plt.title( "First Eigenface" )
renderVector( scaleImageIntensities( Eigenfaces[:,0] ) )
plt.subplot( 233 )
plt.title( "Second Eigenface" )
renderVector( scaleImageIntensities( Eigenfaces[:,1] ) )
plt.subplot( 234 )
plt.title( "Third Eigenface" )
renderVector( scaleImageIntensities( Eigenfaces[:,2] ) )
plt.subplot( 235 )
plt.title( "Fourth Eigenface" )
renderVector( scaleImageIntensities( Eigenfaces[:,3] ) )
plt.subplot( 236 )
plt.title( "Fifth Eigenface" )
renderVector( scaleImageIntensities( Eigenfaces[:,4] ) )
plt.show()
Shape of cov(X): (4096, 4096)
By modifying the coefficient of an eigenface (i.e., $c_j$ for $\mathbf{e}_j$ above), we can alter the appearance of the face. Let's modify the coefficient for the fourth (i.e., $j = 3$) eigenface, which emphasizes the tip of the nose and the eyebrows. In general, the brighter the pixels in an eigenface, the more it emphasizes those pixels/features in a face.
# At this point, we know the mean face m and the Eigenfaces V.
f = plt.imread( "faces/face000.bmp" ) # Let's modify the first face in the data set.
f = imageToVector( f )
x = f - m
c = Vt @ x # Equivalently, Eigenfaces.T @ x: How many coefficients are in c?
factors = [1, 3, 6, 9]
plt.figure( figsize=(10,4) )
subplotCode = 140
for i in range( len( factors ) ):
factor = factors[i]
plt.subplot( subplotCode + i + 1 )
cNew = c.copy()
cNew[3] *= factor
renderVector( m + Eigenfaces @ cNew ) # Recall that Eigenfaces = V.
plt.xlabel( "factor = {}".format( factor ) )
if i == 1:
plt.title( "A face image with 4th Eigenface coefficient multiplied by a constant factor" )
By having our face cookbook ($V$ = eigenfaces) at hand, we can approximate a face $\mathbf{f}$ with just $k$ numbers!
$$\mathbf{f}^{(k)} = \mathbf{m} + V^{(k)}\mathbf{c} = \mathbf{m} + c_0\mathbf{e}_0 + c_1\mathbf{e}_1 + \cdots + c_{k-1}\mathbf{e}_{k-1}.$$############## The kth approximation of a face from the cookbook ##############
k = 50 # You must use k <= 169 because min(row,col)=169.
f = plt.imread( "faces/face000.bmp" ) # Let's take again the first face in the data set.
f = imageToVector( f )
x = f - m
c = Eigenfaces[:,:k].T @ x
# Plotting the resulting kth approximation and the original image.
plt.figure( figsize=(6,3))
plt.subplot( 121 )
renderVector( f )
plt.xlabel( "Original face" )
plt.subplot( 122 )
renderVector( m + Eigenfaces[:,:k] @ c )
plt.xlabel( "Approximated face" )
plt.suptitle( "Approximating a face using k = {} Eigenfaces".format( k ) )
Text(0.5, 0.98, 'Approximating a face using k = 50 Eigenfaces')
############## The kth approximation of a face *not* using Eigenfaces ##############
k = 20
f = plt.imread( "faces/face000.bmp" ) # Let's take again the first face in the data set.
U2, sigma2, V2t = spla.svd( f )
fk = np.zeros( f.shape )
for i in range( k ):
fk += sigma2[i] * np.outer( U2[:,i], V2t[i,:] )
plt.figure( figsize=(6,3))
plt.subplot( 121 )
plt.imshow( f, cmap="gray" )
plt.xlabel( "Original face" )
plt.subplot( 122 )
plt.imshow( fk, cmap="gray" )
plt.xlabel( "Approximated face" )
plt.suptitle( "Approximating a face using k = {} singular vector and no eigenfaces".format( k ) )
Text(0.5, 0.98, 'Approximating a face using k = 20 singular vector and no eigenfaces')
By knowing the mean face $\mathbf{m}$ and the eigenfaces in $V$, we can synthesize any face $\tilde{\mathbf{f}}$ by sampling a vector of coefficients $\tilde{\mathbf{c}}$ from some normal distribution:
$$\tilde{\mathbf{f}}^{(k)} = \mathbf{m} + V^{(k)}\tilde{\mathbf{c}} = \mathbf{m} + \tilde{c}_0\mathbf{e}_0 + \tilde{c}_1\mathbf{e}_1 + \cdots + \tilde{c}_{k-1}\mathbf{e}_{k-1}.$$############## Generating a random face from the cookbook ##############
k = 20 # Use only the top k Eigenfaces.
Coeffs = X @ Eigenfaces[:,:k] # An nxk matrix of face coefficients.
mu = np.mean( Coeffs, axis=0 ) # Mean row vector of coefficients.
std = np.std( Coeffs, ddof=1, axis=0 ) # Standard deviations along each column of coefficients.
print( "Variance from function:" ) # Checking that the variances on coefficients should be the same
print( std ** 2 ) # as the singular values!
print( "\nVariance from SVD:" )
print( sigma[:k] )
# To generate the new face, we must generate new coefficients based on the
# statistics from the population: the mean and std obtained above.
newFaceCoeffs = std * np.random.randn( k ) + mu # This creates coefficients with a normal distribution
newFace = m + Eigenfaces[:,:k] @ newFaceCoeffs # whose mean is mu and variance is sigma^2.
# Plotting new random face.
plt.figure( figsize=(3,3) )
plt.title( "Random face" )
renderVector( newFace )
Variance from function: [594316.78282257 388739.27150947 220110.71998256 215765.40290314 148366.11611318 136747.66077291 83176.09921457 70937.49110233 59014.92801102 55035.63320068 49892.09229926 45004.20836223 40254.16226681 35670.46406856 34394.65172355 32401.94970633 30830.37818509 28070.31809972 26477.49452898 24139.03654889] Variance from SVD: [594316.78282257 388739.27150947 220110.71998256 215765.40290314 148366.11611318 136747.66077291 83176.09921457 70937.49110233 59014.92801102 55035.63320068 49892.09229926 45004.20836223 40254.16226681 35670.46406856 34394.65172355 32401.94970633 30830.37818509 28070.31809972 26477.49452898 24139.03654889]
<matplotlib.image.AxesImage at 0x11064daf0>