Saturday, April 3, 2010

Matrix and linear algebra in F#, Part III: Eigen decomposition and face recognition

The theory

Face recognition is one of the most import research problems in computer vision. It is also an important application for everyday use. A lot of security system has a face recognition component, as well as other parts, e.g., finger print recognition.

While this application has a lot of algorithms, the most famous one must be Eigen face proposed in

M. Turk and A. Pentland (1991). "Eigenfaces for recognition". Journal of Cognitive Neuroscience 3 (1): 71–86.

The theory of Eigen face is very simple if one has learned the concept of Principal Component Analysis (PCA). Given a set of faces, each face is a n-by-m matrix, if we reshape the matrix into a vector, then each face can be viewed as a point in the (n*m) dimensional space. The idea of Eigen face is to find some directions, which are orthogonal, in this space such that the points(faces) have a minimal representation, or have maximum variance along these directions. As you can see, this is an optimization problem. The magic part is that this optimization is solved by Singular Value Decomposition (SVD).

The following is a summary of the procedures to get the eigen faces:

1. for a set of faces, F = {f1,…,fn}, calculate its mean face m. F is a d-by-n matrix, where d is the number of pixels of an face and n is the number of faces in the dataset. fi is a column vector.

2. subtract the mean face from F, get A = {f1-m, … fn-m}. We want to get the eigenvalues and eigenvectors of Cov(A) = A * A’. The eigenvectors associated with big eigenvalues are eigenfaces. (In the program part, we will visualize these eigenfaces).

However, Cov(A) is usually too big. (If there are 10K pixels in each face, then the size of Cov(A) would be 10K-by-10K). A trick is to calculate the eigenvalues and eigenvectors of A’ * A (which is of small size if the number of faces n is small.) and calculate Cov(A)’s based on them. This trick is well documented in the program section.

Once we get the Eigen faces, we can do face recognition. As said before, each eigen face is actually an direction, we can map our new face and the face in the database into the space defined by the set of Eigen faces. For example, if we have 10 Eigen faces, and we know that they are orthogonal to each other, thus a 10-D real space is defined by the 10 Eigen faces. Once real faces are mapped into this 10-D space, we can measure the distance between the new face and the existing faces, and find the closet one as the label face for the new face.

Here are some reference for interested readers:

1. Eigenface wiki page. It has the mathematics and links to programs and papers.

2. A Matab Tutorial. and this one. These two are very good tutorials using Matlab.

3. Principal Component Analysis (PCA), Singular Value Decomposition and Eigendecompsition. If you want to know more theory of it, you can move on to have a good understanding of PCA.

The ORL face dataset

There are more than 10 face datasets available online, which range from different complexity. I use an quite old one the ORL face dataset. This one is well pre-processed, e.g. under the same light condition, of the same size, etc.

Using this dataset greatly eases our work in processing the images and enables us to focus on the Eigen face algorithm in F#.

This dataset contains 400 faces from 40 subjects, with each having 10 faces. The images are of .PGM format. For the I/O of the images, I wrote a separate blog post on reading, writing and visualizing this image format. That post also features on how to write a simple GUI in F#.

The data structures for the dataset:

type pgm = matrix
type person = matrix array
type dataset = person array
Each PGM image is an F# matrix, each person is an array of matrices, and the dataset is an array of persons.

Then we define functions to read the dataset for a given folder:

let readPerson folder =
seq {
for i=1 to 10 do
name = folder + string(i) + ".pgm"
yield readPgm name
} |> Seq.toArray

let readDataset folder =
seq {
for i=1 to 40 do
name = folder + "s" + string(i) + "\\"
yield readPerson name
} |> Seq.toArray
We also define functions to split the dataset into training data and testing data. Training data are used to get the eigenfaces, and then faces in the testing data are recognized using the eigenfaces. For each person, the first 7 faces are used for training, and the 3 left ones are used for testing.

let splitDataSet (ds:dataset) =
let train = seq {
for p in ds do
seq {
for i=0 to 9 do
i<7 then yield p.[i]
} |> Seq.toArray )

let test = seq {
for p in ds do
seq {
for i=0 to 9 do
i>=7 then yield p.[i]
} |> Seq.toArray )
train |> Seq.toArray, test |> Seq.toArray

let strip (data: 'a array array) =
let dat = seq {
for x in data do
y in x do
let label = seq {
let i = ref 0
for x in data do
y in x do
i := !i + 1
dat |> Seq.toArray, label |> Seq.toArray

Function strip is used to transform a data set (matrix array array) into a flat array (matrix array) and remember their labels.

The program

First we write some help functions:

let mat2vec (A:matrix) =
let n = A.NumRows
let m = A.NumCols
Vector.init (n*m) (fun i->A.[i%n, i/n])

let vec2mat (vec:vector) nrow ncol =
if vec.Length <> nrow * ncol then failwith "vec length != nrow * ncol "
Matrix.init nrow ncol (fun i j -> vec.[j*nrow+i])

let colSum (m:matrix) =
let nrow = m.NumRows
Matrix.foldByRow (fun a b -> a+b) (Vector.create nrow 0.0) m

let colMean (m:matrix) =
let ncol = m.NumCols
colSum m * (1.0 / float ncol)
mat2vec is concate the columns of a matrix into a long column vector, while vec2mat does the inverse. Their behaviors are similar to the reshape function in Matlab. colSum and colMean are similar to sum and mean in Matlab.

To use the eigen decomposition for a symmetric matrix, we call lapack via math-provider. Read Part II of this F# matrix series to know to how use them.

open Microsoft.FSharp.Math
open Microsoft.FSharp.Math.Experimental
let isSucc = Experimental.LinearAlgebra.Lapack.Start()
printfn "service = %A" (if isSucc then "Netlib lapack" else "Managed code")
let eig = LinearAlgebra.EigenSpectrumWhenSymmetric
eig is a shorthand for EigenSpectrumWhenSymmetric. Notice that eig gets the eigenvectors in a matrix, each row of which is an eigenvector. This may be a bug in the math-provider wrapper, but since it only deals with square and symmetric matrices, this incontinence does not hurt too much.

The main function:

let eigCov (faces:matrix) numVecs =
    // each column is a face

let nrow = faces.NumRows
let ncol = faces.NumCols
let colm = colMean faces
let A = Matrix.init nrow ncol (fun i j -> faces.[i,j] - colm.[i])
// A - each column is a difference face
// Cov(A) = A * A', a big matrix
// trick: A'*A*v(i) = lambda(i)*v(i)
// A*A'*A*u(i) = lambda(i)*A*u(i)
// (A*u(i)) is an eignvector of Cov(A)
// L(A) = A' * A
let L = A.Transpose * A
let val1, vec1 = eig L // each row of vec1 is an eigenvector
// get the eigen values and eigen vectors for Cov(A)
let v = val1 * (1.0 / float (ncol - 1))
let u = A * vec1.Transpose

// normalize eigen vectors
let mutable cnt = 0
for i=0 to ncol-1 do
abs v.[i] < 1e-8 then
u.[0..nrow-1,i..i] <- Matrix.create nrow 1 0.0
cnt <- cnt + 1
let norm = Vector.norm (u.[0..nrow-1,i..i].ToVector())
u.[0..nrow-1,i..i] <- u.[0..nrow-1,i..i] * ( 1.0 / float norm )

// sort
// let u' = u.PermuteColumns (fun j -> j ) //u.NumCols - j - 1)
// the PermuteColumns function is buggy
let u' = Matrix.init u.NumRows u.NumCols (fun i j -> u.[i, u.NumCols - j - 1])

let v' = v.Permute (fun j -> v.NumRows - j - 1)

let select = min cnt numVecs
//printfn "%A" select
if cnt < numVecs then
printfn "warning: numvec is %d; only %d exist." numVecs cnt

let v1 = v'.[0..(select-1)]
let u1 = u'.Columns(0, select)
let transformed = (u1.Transpose) * A
v1, u1, colm, transformed

let eigenface (p:matrix array) numVecs =
let B = p |> mat2vec
let nrow = B.[0].Length // nrow is big == # of elements in a face
let ncol = B.Length // # of faces for a person
// put a 2d face into a vector as a column in matrix A
let A = Matrix.init nrow ncol (fun i j -> B.[j].[i])

let v, u, mface, transformed = eigCov A numVecs
(new PGMViewer(vec2mat mface 112 92, "mean face")).Show()
for i=0 to 2 do
let eface = scale (u.Column(i))
let f = vec2mat eface 112 92
(new PGMViewer(f, "eigenface " + i.ToString())).Show()

v, u, mface, transformed
This function transform an image matrix into a vector and call the eigCov to get the eigenfaces. The commented lines in this function is used to output the mean face and the first 3 eigenfaces:


Use eigenfaces to do recognition for a new face:

let recognition (u:matrix) (mface:vector) (transformed:matrix) (face:vector) =
    let d = face - mface
let t = u.Transpose * d

let distance a b = Vector.norm (a - b)
|> (fun j -> distance t (transformed.Column(j)), j)
|> Seq.min
|> snd
To put the recognition into a learning and testing framework:

let learning data =
let train, test = splitDataSet data
let dat, lab = strip train
let dat0, lab0 = strip test
let v, d, mface, tran = eigenface dat 10

let correctCnt =
|> Seq.mapi (fun i f ->
idx = recognition d mface tran (mat2vec f)
let label = lab.[idx]
if label = lab0.[i] then 1 else 0
|> Seq.sum
float correctCnt / (float dat0.Length)
Using 10 most significant eigen faces from the 280(7 * 40) faces, we recognize the remaining 120(3 * 40) faces, and get 113 faces recognized to the correct person, thus a 94.17% accuracy, which is reasonable good considering that we have done any normalization on the faces, and the classifier is only 1 nearest neighbor.


In this post, we have seen that how F# could be used as an prototype tool like Matlab. In the matrix computation part, we need to write some utility functions ourselves, while in Matlab they are built in. However, the strength of F# comes from the language part, we could control the training and testing more freely than in Matlab. Moreover, once we have developed this face recognition component, it is just straightforward to use C# or VB.Net to call it.

Although we have a high accuracy in the final test, the result is only for a well prepared dataset. In real situations, things would be very hard. Maybe over 90% time is spent on tuning best parameters and preprocessing the face images.




  1. When someone writes an article he/she keeps the thought of a user in his/her brain that how a user can be aware of it. Therefore that’s why this post is outstanding.Thanks!


  2. I really enjoyed reading your article. I found this as an informative and interesting post, so i think it is very useful and knowledgeable. I would like to thank you for the effort you have made in writing this article.