Bug when compiling Mex-files

Ones might encounters problems from compilers when trying to compile MEX files. Assume that you are in Linux, with installed GCC/G++ 4.x and MATLAB R2009x or R2010a. The machine can announces that some libs are absent, i.e. libstdcXXX but the fact is that you already install everything! Okay, the problem is just that there are small problems in the configuration file mexopts.sh, the bash file mentions everything about OS environment for MATLAB. By default, Mathworks encapsulates compilers with Matlab versions, and this compiler version soon to be outdated from the regular update of compiler, i.e. GCC. So, it might cause version conflict when MATLAB tries to find an exact compiler having version XYZ while your machine just have version ZYZ.01 which is more recent.

The solution for this is quite simple. You do not need to create any further symbolic links. Just go into your home directory by typing cd ~; then cd .matlab; then cd R20XXX; then open file mexopts.sh and start modifying.

Depending on your machine which is glnx86 or glnxa64, you must change the lines CC and CXX from gcc-x.xx.x to gcc, and g++-x.xx.x to g++. That’s it! Save it and enjoy Mex files.

Remark: every time you call mex -setup, this file is overrided, then you have to do this stuff again.

Constrained Quadratic Programming Solver – qpOASES

Xin giới thiệu cho những ai muốn làm optimization với constrained quadratic programming, toolbox free qpOASES. Thực sự số lượng toolbox tốt cho optimization không nhiều. Ví dụ nếu muốn solve một cái problem nào thì xài Matlab cracked có thể OK nhưng muốn tăng tốc thì hơi bị chua. Giải pháp là kiếm cái lib chạy trên C/C++. Theo thống kê trên wiki thì số toolbox xài free được rất hiếm.

May mắn là phát hiện ra chú này,

http://www.kuleuven.be/optec/software/qpOASES

giải rất kool. Biên dịch cũng cực kỳ đơn giản vì chả cần third party libs nào hết. Thằng này có interface cho cả matlab và trong option của nó có nhiều lựa tuỳ biến cho solver. Qua su dung thay thang nay kha robust!

Có điều cần lưu ý là nó không tolerable bằng Matlab opt. toolbox nhá! Nó dễ gặp các vấn đề về numerical, nếu data của mình ko được hiệu chỉnh kỹ lưỡng. Cái này thì problem dependent, nên ko thể có một trick chung được. Tuy nhiên, ko nên để magnitude quá chênh lệch giữa Hessian matrix H và f:

Thông thường, nhớ kiểm tra điều kiện semi-positive definite của H. Có thể kiểm tra nhanh bằng cách dùng hàm eigs của Matlab để tìm eigenvalues của H. Nếu tồn tại một negative eigenvalues thì coi như tiêu.Trường hợp H suy biến (singular/ill-conditioned) thì nên cộng vào một lượng bias nhỏ dọc đưòng chéo chính của H:

</pre>
H = H + diag(1e-4*ones(size(H,1)));

Dictionary Learning with CVX Toolbox

In this post we do an exercise that help you implement sparse coding and dictionary learning without using algorithms dedicated to optimizing a -regularized least squares problem. We will use a generic convex optimization CVX toolbox for sparse encoding and adopt simple updating rule for dictionary .

For simplicity, a single image “lena” is used to learn the dictionary. Firstly, patches are extracted and normalized:

clear all;
I=double(imread('lena.png'))/255;
% extract 8 x 8 patches
X=im2col(I,[8 8],'sliding');
X=X-repmat(mean(X),[size(X,1) 1]);
X=X ./ repmat(sqrt(sum(X.^2)),[size(X,1) 1]);
X=X(:,1:100:end);

Since the number of patches is so huge, we just sub-sample 1%. This is enough for us to check the correctness of our algorithm.

Next, the number of iterations, dictionary size, and regularization coefficient are set:

nIter = 25;
gamma = .15;
p = 200;

% init dictionary
sel = randperm(size(X,2));
sel = sel(1:p);
D = X(:,sel);
n = size(X,2);
%D = D ./ repmat( sqrt(sum(D.^2)), [w^2, 1] );
S = zeros(p,n);

We individually encode each of input sample by optimizing the following objective function:

cvx_begin quiet
variable s(p);
minimize (.5*norm(X(:,j)-D*s) + gamma*norm(s,1));
cvx_end

In order to update , we just need to solve a least squares problem:

This objective function can be solved easily by taking the first-order derivative respect to :

If the Hessian is positive semidefinite, the above quadratic form reachs its global minimum value with the minimizer:

There is the case that will grow big and shrinks small to fit . In order to prevent it from happenning, a -norm applied for dictionary atoms is neccessary:

The Matlab code is very simple:

D = X * pinv(S);
D=D ./ repmat(sqrt(sum(D.^2)),[size(D,1) 1]);

After iterations, we can visualize the content of dictionary atoms by calling from SPAMS Toolbox:

ImD=displayPatches(D);
figure, imagesc(ImD); colormap('gray');

The result after 6 hours training (!) is displayed below:

Learned dictionary from Lena image with 100 atoms

Visual Object Dictionary Learning based on Sparse Coding (part II)

Now, let’s take the traditional approach of spliting images into patches and learn them. Finishing this experiment, we can understand a little about the behaviour of dictionary learning in recognizing visual objects.

Image used for training

I took patches of pixels with sliding fashion just on one image. If you do with more images, the total number of patches is very huge! So, my training set is ~250k samples.

Dictionary Visualization

You can observe that the dictionary atoms look quite similar to the experiments done in papers. This is the good news since our code work correctly.

-norm magnitude over

Comments: It is quite strange when -norm monotically increases when changes from very small to bigger. Normally, is the regularization term on the magnitude of coefficients . When increases, the magnitude of penalized term should decrease. But the above figure does not show such the thing.

Reconstruction error of the fidelity term

Comments: This figure also mentions strange things. The regularization term is used to prevent overfitting, in which the result regression function fits perfectly with the data. Therefore, the residual should increase when the regularization term is put more weights. In contrast, residual decreases when decreases!

Percentage of non-zeroes coefficients over

Comments: The sparsity of the encoded vectors also decreases! What happening?

My correction: if you cross out line 38 and 48 in the code snippet below, it will produce correct figures with the  increase of cost-function and decrease of non-zeroes when grows up. Since I am lazy, you correct for yourself, in Part I and Part II too. In fact, I think the documentation of SPAMS Toolbox has some sort of confusion in describing running modes. Anyway, it works well if we does not explicitly declare any mode.

% load images and extract patches
imglist = dir('dat/smallawa/');
if length(imglist) < 3
disp 'empty directory';
end

% just take first 5 images, otherwise it's huge!
imglist = imglist(3:4);
% imgsize = [32 32];
bag = cell(1,length(imglist));

wndsize = [32 32];

for j=1:length(imglist)
I = imread(['dat/smallawa/' imglist(j).name]);
I = rgb2gray(I);
I = im2double(I);
%I = imresize(I,[32 32]);
%X(:,j) = reshape(I,1, []);
X=im2col(I,wndsize,'sliding');
bag{j} = X;
end

X = cell2mat(bag);

disp (['Total patches: ' num2str(size(X,2))]);
X = X - repmat(mean(X),[size(X,1) 1]);
X = X./ repmat(sqrt(sum(X.^2)),[size(X,1) 1]);

lambda = logspace(-2,1,10);
R = zeros(length(lambda),1);
NR = zeros(length(lambda),1);

for i=1:length(lambda)
param.K = 100;
param.lambda = lambda(i);
param.iter = 100;
param.mode = 2;
param.modeD = 0;

tic;

D = mexTrainDL_Memory(X,param);
t=toc;
fprintf('time of computation for Dictionary Learning: %f\n',t);

fprintf('Evaluating cost function...\n');
param.mode = 0;
alpha=mexLasso(X,D,param);
alpha(isnan(alpha)) = 0;
R(i)=mean(0.5*sum((X-D*alpha).^2)+param.lambda*sum(abs(alpha)));
ImD=displayPatches(D);
%   title (['lambda=' num2str(lambda(i))]); axis square; axis off;
imwrite(ImD,sprintf('%02i.png',i));
NR(i) = sum(sum(alpha>0))/size(X,2);
%disp (['Non-zeros ' num2str(NR(i)*100), '%']);
%disp (['Residual:' num2str(R(i))]);
end

figure,semilogx(lambda,R); xlabel('lambda'); ylabel('residual'); grid on;
figure,semilogx(lambda,NR); xlabel('lambda'); ylabel('non-zeroes'); grid on;

Visual Object Dictionary Learning based on Sparse Coding (Part I)

In sparse coding for image denoising or inpainting, dictionary is learned over the patch level, i.e. images are divided into rectangular patches and the unsupervised dictionary learning take those to induce a dictionary for sparse coding. The famous example of coding natural images shows that the learned dictionary contains “edges” patches, a higher visual level compared to pixel. I wonder if patches are the whole images, how the dictionary looks like.

I did a funny experiment with Corel dataset. A small set of training images were taken from “tiger” category. There are ~500 images and I used all for dictionary training.

Firstly, I prepared training data by doing some preprocessing: convert to grayscale colormap, rescale it to [0..1], center the data by their mean, normalize with L2 norm.

Secondly, I set parameters for the algorithm dictionary provided by SPAMS Toolbox and wait for minutes. Because the training set is small, it was not take much time to finish the training.

Thirdly, I recorded some useful informations such as the number of non-zeroes and the residuals (between original signal and the reconstructed one).

Note that I train the dictionary over couples of values of , says from . It is interesting to observe the correlations between , non-zeroes, and dictionary visualization.

Dictionary Visualization

Comments:

Residuals over lambda

Comments:

Average number of non-zeroes per vector

Comments:

Here is the Matlab code for my experiment

% load images and extract patches
imglist = dir('dat/smallawa/');
if length(imglist) < 3
disp 'empty directory';
end imgsize = [32 32];
X = zeros(prod(imgsize),length(imglist)-2);
for j=1:length(imglist)-2
I = imread(['dat/smallawa/' imglist(j+2).name]);
I = rgb2gray(I); I = im2double(I);
I = imresize(I,[32 32]);
X(:,j) = reshape(I,1, []);
end
disp (['Total images = ' num2str(size(X,2))]);
X = X - repmat(mean(X),[size(X,1) 1]);
X = X./ repmat(sqrt(sum(X.^2)),[size(X,1) 1]);
lambda = logspace(-2,1,20);
R = zeros(length(lambda),1);
NR = zeros(length(lambda),1);
for i=1:length(lambda)
param.K = 64;
param.lambda = lambda(i);
param.iter = 500;
param.mode = 2;
param.modeD = 0;
tic;
D = mexTrainDL_Memory(X,param);
t=toc;
fprintf('time of computation for Dictionary Learning: %f\n',t);
fprintf('Evaluating cost function...\n');
param.mode = 0;
alpha=mexLasso(X,D,param);
R(i)=mean(0.5*sum((X-D*alpha).^2)+param.lambda*sum(abs(alpha)));
ImD=displayPatches(D);
%   title (['lambda=' num2str(lambda(i))]); axis square; axis off;
imwrite(ImD,sprintf('%02i.png',i));
NR(i) = sum(sum(alpha>0))/param.K;
%disp (['Non-zeros ' num2str(NR(i)*100), '%']);
%disp (['Residual:' num2str(R(i))]);
end

figure,semilogx(lambda,R); xlabel('lambda'); ylabel('residual'); grid on;
figure,semilogx(lambda,NR); xlabel('lambda'); ylabel('non-zeroes'); grid on;

Partiview & Ndaona: A good way to study Machine Learning

Partiview is a piece of software for interactive 4D-viewer dataset. It is also useful for 3D visualization of dataset in Machine Learning algorithm, i.e. show manifolds or classification results.

Ndaona is a piece of Matlab functions for Partiview-compatible data conversion.

Downloads

• PartiView can be reached here (binaries + source codes)
• Ndaona is here

Screenshots

I borrow some screenshots from Ndaona’s website to illustrate here:

How well a linear svm (in LIBSVM) classifies Mandarin tones

How well does Principal Components Analysis separate faces

[Fastfood] Digital Image Processing + Matlab