Description

The purpose of this homework is to compare matrix factorization techniques through a facial recognition study. For this purpose, the CBCL database, originally used in [3], which is made available in Moodle will be used. Note that the dataset is divided into train and test datasets. In order to perform your first task, you will use the train dataset. The techniques from [1] will be used in order to carry out the comparison. Proceed as described in the following steps and report your findings and results as required. Note that the required programming can be implemented in MATLAB or Python.
• Download the facial database from your Moodle account.
• Describe how the matrix X ∈ Rp×n (see [1]) should be formed.
1 Implementation and analysis
The following steps correspond to the comparison between Singular Value Decomposition (SVD) and Non-negative Matrix Factorization (NMF).
1. Singular Value Decomposition (SVD)
(a) Use SVD to factorize X as
X = UΣV T ,
where Σ is a 361 × 361 square matrix with diagonal entries σi2, 1 ≤ i ≤ 361 = 19 × 19.
(b) Plot the singular values as plot([1:361],diag(Σ)). Then compute the accumulated energy e of the singular values as: e , for 1 ≤ i ≤ 361. Afterwards, plot the normalized accumulated energy as plot([1:361],e/max(e)), where max(e) represents the maximum element of the vector e.
(c) Identify the indices I90, I95 and I99 of e that correspond to the indices when the normalized energy reaches to 0.9, 0.95 and 0.99, respectively.
(d) Check the corresponding singular faces which are the first I90 columns of U and display them as in Figure 1 of [1]. Comment on the obtained singular faces? Do they have localized or distributed features? Are they non-negative valued?
2. Non-negative Matrix Factorization (NMF)
(a) Use the Non-negative Matrix Factorization (NMF) technique called HALS in Section 3.1.5 of [1] to factorize X = WH, where both W ∈ R361×r and H ∈ Rr×2429 have non-negative entries, for some r ∈ N to be determined later. For this purpose do the following:
• Write a code that implements the HALS update given in Section 3.1.5 of [1].
• Initialize the algorithm by using the SVD-based initialization as described in Section 3.1.8 of [1]. Here you will use the decomposition that you found in the previous section in the following way. First, choose the value of r using SVD or by trial and error if the SVD choice does not seems satisfactory. Now let us define [x]+ = max(x,0) and [x]− = max(−x,0). Furthermore, for all 1 ≤ k ≤ r, we define [Mk]+ = [uk]+[vkT]+, where uk are columns of U and vk are columns of V . Similarly, we define [Mk]− = [uk]−[vkT]−. If k[M]+kF ≥ k[M]−kF, initialize the columns of W and rows of H as:
wk = [uk]+/k[uk]+k and h .
Otherwise, initialize them as:
wk = [uk]−/k[uk]−k and hTk = σk k[uk]−k[vk]−.
• As a Stopping criterion of the HALS iterations you can use Section 4.2.2 of [2].
• Plot the convergence of iterations as given in Figure 3 of [1]. Note that your plot will contain a single curve corresponding to the HALS results. Compare these results with the corresponding ones given in Figure 3 of [1].
(b) Check the eigenfaces in W and display them as images as in Figure 1 of [1]. Comment on the obtained eigenfaces? Do they have localized or distributed features? Are they non-negative valued?
2 Image recovery from noisy data
Your next task is to accurately reconstruct images from their noisy versions. For this, the test dataset will be used. Let us denote by Yk ∈ R19×19 the kth image from the test dataset, where 1 ≤ k ≤ 472. First, you need to add random noise to this image in order to obtain Yek, where the (i,j) entry of Yek is defined as
rand,255},
where η is a scaling parameter. You need to perform the following task for the values η = 1,10,25 and comment on the results.
Then, the noisy image should be vectorized. In order to do that, let us formally define the following function:
vec : Rm×n −→ Rmn
    v1
v1 … vn →7  …  ,
   
vn
where vi ∈ Rm, 1 ≤ i ≤ n are the columns of the input matrix. Now let us denote by yk the vectorized version of
Yk, that is, yk = vec(Yk). Similarly, we have yek = vec(Yek).
First, you should reconstruct the original image by using SVD. Recall the definition of I90 in the previous section. Let us denote it by rSV D from here on. You will use the first rSV D singular faces from the matrix U obtained above in order to reconstruct the original images. What we have so far are the noisy vectorized images yek, 1 ≤ k ≤ 472. For each yek, we denote by its reconstructed version. That is:
,
where is the dot product of yek and up and up are the singular faces of X (columns of U), for every
1 ≤ p ≤ rSV D.
Similarly, you should reconstruct the original image by using NMF. Denote by rNMF the value of r that you chose for its implementation. Therefore, we have the weight matrix W ∈ R361×rNMF . The reconstructed vector will in this case be:
,
where wp are the eigenfaces of X (columns of W) and β = (WTW)−1WTyek is the least squares solution of W corresponding to yek; note that βp is, as the pth element of β, a scalar.
Next, we transform ySV Dk and yNMFk back to a matrix via the following function:
mat : Rmn × R −→ Rm×n
     v1
 … , m 7→ v1 … vn .
     vn
Note that mat takes as input the vector and the the size m of the columns of the output matrix. We denote the SVDreconstructed image by Y kSV D, i.e. mat(ySV Dk ,19), and the NMF-reconstructed image by Y NMFk , i.e., mat since we want a 19 × 19 image. The following diagrams describe the whole procedure
for the SVD-based and NMF-based reconstruction respectively.
Yk −−−−−−−−+rand(19,19)→ Yek −−−−−vec( k→ ye −−−SV D→ ySV D mat ;
+rand(19,19) vecmat
Yk −−−−−−−−→ Yek.
SV D
For every reconstructed image Y k and Y k , you will compute the Frobenius norm of the corresponding error
EkSV D and EkNMF respectively, where the error is defined as Ek = Yk − Y k. So for each k, the norm and should be computed. The final step is to plot the average of the errors over the 472 test images versus different values of rSV D and rNMF. So the y-axis of the first graph contains the values of
error
corresponding to different values of rSV D in the x-axis. In the second graph you will show the values of error(rNMF) for different values of rNMF similarly.
Comment on the results and compare. What can you say about the change of the average error as rSV D and rNMF change? Is there a trend?
3 Image recovery from masked data
For the next task, you need to recover images from their masked versions. A masked image is the entry-wise product of the original image matrix with a shadow mask matrix S. For this homework, the (i,j) entry of the matrix S is defined as follows:
if j ≤ 10
S
if j > 10 Using the same notation as in the previous section we now have
,
where we denote by the entry-wise product, that is, Yek(i,j) = Yk(i,j) · S(i,j), for all 1 ≤ i,j ≤ 19. Next you
SV D NMF should proceed as in the previous section in order to obtain Y k and Y k corresponding to both methods. Note that in this case, the diagrams will look as follows:
vecmat
Yk;
vecmat
Yk.
Finally, you will again plot two graphs showing the change in error(rSV D) and error(rNMF) for different values of rSV D and rNMF respectively.
Comment on the results and compare. What can you say about the change of the average error as rSV D and rNMF change? Is there a trend?
References
[1] N Gillis. The why and how of nonnegative matrix factorization, 2014. arXiv Preprint arXiv:1401.5226v2, 2014.
[2] Nicolas Gillis et al. Nonnegative matrix factorization: Complexity, algorithms and applications. Unpublished doctoral dissertation, Université catholique de Louvain. Louvain-La-Neuve: CORE, 2011.
[3] Daniel D Lee and H Sebastian Seung. Learning the parts of objects by non-negative matrix factorization. Nature, 401(6755):788–791, 1999.