Pattern recognition and machine learning - Bernoulli mixture model

I have been reading the book Pattern Recognition and Machine Learning (Bishop) for a while, and recently I came across this figure, which was created using Bernoulli mixture model on the MNIST dataset: I figured it would be fun to code this, so I basically followed their algorithm:

• Suppose you have N images of handwritten digits from 2 to 4. Let D be the number of pixels of the image (each MNIST image is 28 x 28, so D = 28 * 28)

• You may want to model this image using a Bernoulli distribution for each pixel (remember a Bernoulli distribution with parameter 'mu' is just like flipping a damaged coin that has the probability mu of landing heads, and 1 - mu for landing tails). So you might toss a coin for each pixel, each with a different mu, if it lands head you set the pixel to 1, otherwise you set it to 0. But with this model, you can easily see that each pixel is independent of each other and there is no way you can model handwritten digits with this.

• Here is where the Bernoulli mixture model comes into play. Instead of using a Bernoulli for each pixel, we use a mixture of Bernoullis (that is, a weighted sum of Bernoullis), and this can be solved by using an algorithm called Expectation - Maximization. With this model, modelling digits suddenly becomes possible.

Of course this is just an informal treatment of the method, and for further readings you should definitely check out Chapter 9, Pattern recognition and machine learning, Section 9.2: Mixture of Bernoulli distribution

Let me summarize my results: I was able to recreate the figures in the book, and when I used 600 training images, I got 75% accuracy on the test set. But when I used all the training images that were from 2 to 4, the results were 90% (not bad for my first project, actually). Also the algorithm apparently worked very well with classifying the number 4, but it made many errors at number 2 (overfitting, anyone?).

However the code is not vectorized, and I fear that my coding style is not good. So any feedback would be great. Here is my code:

function [Labels] = ReadLabelsMNIST(filename)
% Read the labels of the MNIST dataset
% Written by Dang Manh Truong

fp = fopen(filename, 'rb');
assert(fp ~= -1, ['Could not open', filename, '']);

magic = fread(fp, 1, 'int32', 0, 'ieee-be');
assert(magic == 2049, ['Bad magic number in ', filename, '']);

numLabels = fread(fp, 1, 'int32', 0, 'ieee-be');
assert(size(Labels,1) == numLabels, 'Mismatch in label count');

fclose(fp);
end


function [images, Labels, numRows, numCols] = LoadMNIST(SelectedNumbers, type, numImages)
% Preprocessing code for Bernoulli mixture model
% Load a number of random images from the MNIST dataset
% All the digits are within the SelectedNumbers row array
% Written by Dang Manh Truong

% type = 1 : Train set. type = 2 : Test set
assert( (type == 1) | (type == 2), 'Type = 1 or Type = 2');

if type == 1
filename = 'train-images.idx3-ubyte';
else
filename = 't10k-images.idx3-ubyte';
end

fp = fopen(filename, 'rb');
assert(fp ~= -1, ['Could not open ', filename, '']);

magic = fread(fp, 1, 'int32', 0, 'ieee-be');
assert(magic == 2051, ['Bad magic number in ', filename, '']);

[~] = fread(fp, 1, 'int32', 0, 'ieee-be');
numRows = fread(fp, 1, 'int32', 0, 'ieee-be');
numCols = fread(fp, 1, 'int32', 0, 'ieee-be');

% Find 'numImages' random images that are from 2 to 4
% TODO: Replace 'find' with logical indexing
% Index = find(2 <= Labels & Labels <= 4);
Index = find(sum(bsxfun(@eq,Labels, SelectedNumbers),2) > 0);

s = RandStream('mt19937ar','Seed',0);
Permuted = randperm(s,size(Index,1));
% Permuted = randperm(size(Index,1));

Permuted = Permuted(1: numImages);
Index = sort(Index(Permuted));
Labels = Labels(Index);

images = zeros(numRows, numCols, numImages);
prev = 0;
ImageSize = numCols * numRows;
for i = 1 : numImages
% Ignore unneeded images
fread(fp, (Index(i) - prev - 1) * ImageSize , 'unsigned char');
Temp = reshape(Temp, numRows, numCols);
images(:,:,i) = Temp;
prev = Index(i);
end

fclose(fp);

images = permute(images,[2 1 3]);
% Reshape to #pixels x #examples
images = reshape(images, size(images, 1) * size(images, 2), size(images, 3));
% Convert to double and rescale to [0,1], then binarize the images
images = double(images) / 255;
images(images < 0.5) = 0;
images(images >= 0.5) = 1;

end


Function to help find cluster for each image. Used in testing phase:

function [Cluster] = GetClusterBMM(images,mu, K)
% Input : images : The input images (N x D)
%       : mu : Parameters for Bernoulli distribution for each pixel (K x D)
%       : phi: Mixing coefficients (K x 1)
%       : K: Number of mixtures
% Output: Cluster (N x 1) The clusters most likely to be associated with
% each image

N = size(images,1);
ClusterSum = zeros(N,K);
Temp1 = mu;
Temp2 = 1 - Temp1;
for n = 1 : N
for k = 1 : K
% In:  http://blog.manfredas.com/expectation-maximization-tutorial/
% They used the sum of mu's, but to be honest I don't know why
% Anyway it only gives about 80% accuracy, while mine gives 90%
%ClusterSum(n,k) = sum(Temp1(k,images(n,:) == 1)) + sum(Temp2(k,images(n,:) == 0));
ClusterSum(n,k) = prod(Temp1(k,images(n,:) == 1)) * prod(Temp2(k,images(n,:) == 0));
end
end
[~,Cluster] = max(ClusterSum,[],2);


Training the model:

function [Correct, MisClassified] = TestBMM(X, TestX, mu, Labels, TestLabels)
% Testing phase for Bernoulli mixture model
% Written by Dang Manh Truong
% The parameters here closely resemble those in the book Pattern
% recognition and machine learning (Bishop), chapter 9
% N : Number of data points
% N': Number of test points
% K : Number of mixtures
% D : Dimension of each data points
% Input: X (N x D) Train data
%      : TestX (N' x D) Test data
%      : mu (K x D) Bernoulli parameters learned from training phase
%      : numTestImages Number of test data needed
%      : Labels (N x 1) Train labels
%      : TestLabels (N' x 1) Test labels
% Output: Correct: The number of times the algorithm get it right
%       : MisClassified(10,10) : The misclassification matrix
% MisClassified(i,j) : The number of times that the digit 'i' is
% misclassified as digit 'j'. Of course the diagonal is zero

K = size(mu,1);
N = size(X,1);
numTestImages = size(TestX,1);
Correct = 0;
MisClassified = zeros(10,10);
digitsInTheSameCluster = zeros(10,numTestImages);

TrainClusters = GetClusterBMM(X,mu,K); % N x 1
TestClusters = GetClusterBMM(TestX,mu,K);
for i = 1 : numTestImages
for n = 1 : N
if TestClusters(i) == TrainClusters(n)
digitsInTheSameCluster(Labels(n),i) = digitsInTheSameCluster(Labels(n),i) + 1;
end
end
[~, AssignedLabel] = max(digitsInTheSameCluster(:,i));
if AssignedLabel == TestLabels(i)
Correct = Correct + 1;
else
MisClassified(TestLabels(i),AssignedLabel) = MisClassified(TestLabels(i), AssignedLabel) + 1;
end
end


Testing the model:

function [mu, phi, Res, effNum] = TrainBMM(X, mu, phi, Res, effNum)
% Training phase for Bernoulli mixture model using Expecation-Maximization
% Written by Dang Manh Truong
% The parameters here closely resemble those in the book Pattern
% recognition and machine learning (Bishop), chapter 9
% D: The dimension of each data points
% N: The number of data points
% K: The number of Bernoulli mixtures
% Input: X (N x D) Data points to be processed (each row - a data point)
%      : mu(K x D) Bernoulli parameters for each mixture
%      : phi(K x 1) Mixing coefficients for each mixture
%      : Res(N x K) Responsibilities of each component (1 - K) given a data
%      point (1 - n)
%      : effNum(K x 1) Effective number of observations for each mixture
% Output: The new values of mu, phi, Res and effNum
% Most of the time only mu will be used

% Size of each image. I don't want to pass these to the function
% because their only purpose is to show the images
numRows = 28; numCols = 28;

N = size(X,1); K = size(phi,1);
iterNum = 0;
uniform = 1 / K;
fprintf('E-M algorithm in progress. This may take a while.....\n');
while 1
% E-step

% Equivalent unvectorized code:
%     for n = 1 : N
%         for k = 1 : K
%             Res(n,k) = 1;
%             for i = 1 : D % D = size(X,2)
%                 if X(n,i) == 1
%                     Res(n,k) = Res(n,k) * mu(k,i);
%                 else
%                     Res(n,k) = Res(n,k) * (1 - mu(k,i));
%                 end
%             end
%         end
%    end
% TODO: Vectorize this part ASAP!!!!
for n = 1 : N
for k = 1 : K
Temp1 = mu(k,:);
Temp2 = 1 - mu(k,:);
Res(n,k) = prod(Temp1(X(n,:) == 1)) * prod(Temp2(X(n,:) == 0));
end
end
Res = bsxfun(@times, Res, phi');
% Divide by the denominator
Sum = sum(Res,2);
Sum(Sum == 0) = uniform;
Res = bsxfun(@rdivide, Res, Sum);

% M-step
effNum = sum(Res,1);
mu = Res' * X;
mu = bsxfun(@rdivide, mu, effNum');
% Check for convergence
iterNum = iterNum + 1;

for k = 1 : K
subplot(1,K,k)
Result = reshape(mu(k,:), numRows, numCols);
subimage(Result)
end
hold on
pause(1)
fprintf('Iteration %d \n',iterNum);

if iterNum >= 10
break;
end
end
fprintf('Press any key to continue\n\n\n');
pause
close(gcf)

end


The main function:

function [] = BMM()
% Bernoulli mixture model for classification of MNIST dataset
% Based on Figure 9.10 in the book Pattern recognition and machine learning
% Partly inspired by: http://blog.manfredas.com/expectation-maximization-tutorial/
% Written by Dang Manh Truong
% I was able to reproduce the 3 pictures of digits 2,3 and 4, and the
% results for classfication on the MNIST test set (2 to 4) were satisfying
% with about 90% correct, not too bad for my first project.
% But it's not a good score compared with state-of-the-art methods

pause(1);
fprintf('Bernoulli mixture model using Expectation - Maximization\n');
fprintf('to recreate Figure 9.10, chapter 9, Pattern recognition and machine learning\n');
% Change these lines if you wish
SelectedNumbers = [2 3 4 ]; % The numbers that we care about in the dataset
numTrainImages = 600; % For the numbers from 2-4: <= 17391
numTestImages = 3024; % For those from 2-4: <= 3024
rng(0,'twister');

% Step 1: Initialization

% images: #pixels * #examples
[images, Labels, numRows, numCols] = LoadMNIST(SelectedNumbers, 1, numTrainImages);

N = numTrainImages; % N : The number of train images
K = size(SelectedNumbers,2); % K : The number of mixtures
D = numRows * numCols; % Dimension of each image
phi = ones(K,1) * 1/K; % Mixing coefficients
mu = (0.75-0.25) * rand(K,D) + 0.25  ; % Means of each components
mu = mu ./ repmat(sum(mu,2),1,D);
Res = zeros(N,K); % Res(k,n): Responsibility of component 'k' given data point X(n,:)
effNum = zeros(K,1); % Effective number of data points associated with each component
X = images'; % (N x D) Each row is an image

% Step 2: Expectation - Maximization
[mu, ~, ~, ~] = TrainBMM(X, mu, phi, Res, effNum);

% Step 3: Testing

[TestImages, TestLabels ,~, ~] = LoadMNIST(SelectedNumbers, 2,numTestImages);
TestX = TestImages'; % Each image is in one row

[Correct, MisClassified] = TestBMM(X, TestX, mu, Labels, TestLabels);

fprintf('Correct: %f percents \n',100 * Correct / numTestImages);
fprintf('The misclassification matrix: \n');
MisClassfied = MisClassified(SelectedNumbers, SelectedNumbers)

% The MisClassified matrix when I used all 17931 train images (that are
% from 2 to 4) to test all 3024 test images (again, from 2 to 4):

%  0   130    48
% 84     0    21
% 17     1     0
% The algorithm appears to correctly labels all the digits 4, but fails
% spectacularly at the digit 2.


To run this program you need the MNIST dataset, available here: http://yann.lecun.com/exdb/mnist/ .Just type "BMM()" in the Matlab command line, and the program will run.

Here are the 3 resulting clusters: • Good job! @reviewers in the close queue: these close votes were made against a prior revision. – Mathieu Guindon Jun 18 '16 at 16:27

Your code is, as far as I can tell, very well written! I can't find much I would do differently.

Capitalization of variables is a matter of choice in MATLAB (there is no convention), but it's not common to use a capital first letter in variable names. Variables with capital letters in the beginning are usually used for methods, properties or the like. The exception is single character variable names (such as A).

Have you used the Profiler when executing your code? Is the bottleneck in the part of the code you have identifies as "Vectorize this ASAP"? The reason why I'm asking is because for loops aren't that slow in MATLAB anymore. The new and improved Execution Engine (JIT). In fact, the performance of loops have improved by 40% on average, making it competative against vectorized approaches.

The performance benefit of JIT compilation is greatest when MATLAB code is executed additional times and can re-use the compiled code. This happens in common cases such as for-loops or when applications are run additional times in a MATLAB session

Have a look at this answer to see some benchmarking results comparing vectorized appraches with for loops. A well written loop can easily outperform a badly written vectorized approach.

You might want to check out varargout. It can be used in your TrainBMM function, as you're only using the two first output arguments (in one of the function calls at least). I don't recommend that you change it, but it can be useful in other scenarios.

I didn't manage to find a way to vectorize the code. Neither did the users in the Matlab/Octave chat room on SO. A final tip: Pre-allocate memory for Res, it might help a bit.

• Indeed the "Vectorize this ASAP" takes the most time in my code. In fact my original version used 3 nested for loop, which was much slower when I checked with Profiler. But still the current version is quite slow, you can use Profiler to verify that, and it would be great to find a way to make it faster, whether it be using a loop or not. Definitely going to check out "varargout" by the way :D – Dang Manh Truong Jun 19 '16 at 12:18
• @DangManhTruong, I have asked in the Matlab/Octave chat room on SO, and no one there managed to vectorize your code. The mixed indexing makes it hard. Sorry I couldn't be of more help. I tried! :) – Stewie Griffin Jun 29 '16 at 8:38
• Tks anyway :) ! – Dang Manh Truong Jun 29 '16 at 9:52