diff DL/dl_ramirez.m @ 177:714fa7b8c1ad danieleb

added ramirez dl (to be completed) and MOCOD dictionary update

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DL/dl_ramirez.m	Thu Nov 17 11:18:25 2011 +0000
@@ -0,0 +1,189 @@
+function DL = dl_ramirez(Problem,DL)
+%% Dictionary learning with incoherent dictionary
+%
+% REFERENCE
+% I. Ramirez, F. Lecumberry and G. Sapiro, Sparse modeling with universal
+% priors and learned incoherent dictionaries.
+
+%%
+%   Centre for Digital Music, Queen Mary, University of London.
+%   This file copyright 2011 Daniele Barchiesi.
+%
+%   This program is free software; you can redistribute it and/or
+%   modify it under the terms of the GNU General Public License as
+%   published by the Free Software Foundation; either version 2 of the
+%   License, or (at your option) any later version.  See the file
+%   COPYING included with this distribution for more information.
+
+%% Test function
+if ~nargin, testdl_ramirez; return; end
+
+%% Parameters & Defaults
+X = Problem.b;					%matrix of observed signals
+
+% determine dictionary size %
+if (isfield(DL.param,'initdict')) %if the dictionary has been initialised
+	if (any(size(DL.param.initdict)==1) && all(iswhole(DL.param.initdict(:))))
+		dictSize = length(DL.param.initdict);
+	else
+		dictSize = size(DL.param.initdict,2);
+	end
+end
+if (isfield(DL.param,'dictsize'))
+	dictSize = DL.param.dictsize;
+end
+
+if (size(X,2) < dictSize)
+	error('Number of training signals is smaller than number of atoms to train');
+end
+
+
+% initialize the dictionary %
+if (isfield(DL.param,'initdict')) && ~isempty(DL.param.initdict);
+	if (any(size(DL.param.initdict)==1) && all(iswhole(DL.param.initdict(:))))
+		D = X(:,DL.param.initdict(1:dictSize));
+	else
+		if (size(DL.param.initdict,1)~=size(X,1) || size(DL.param.initdict,2)<dictSize)
+			error('Invalid initial dictionary');
+		end
+		D = DL.param.initdict(:,1:dictSize);
+	end
+else
+	data_ids = find(colnorms_squared(X) > 1e-6);   % ensure no zero data elements are chosen
+	perm = randperm(length(data_ids));
+	D = X(:,data_ids(perm(1:dictSize)));
+end
+
+
+% coherence penalty factor
+if isfield(DL.param,'zeta')
+	zeta = DL.param.zeta;
+else
+	zeta = 0.1; 
+end
+
+% atoms norm penalty factor
+if isfield(DL.param,'eta')
+	eta = DL.param.eta;
+else
+	eta = 0.1; 
+end
+
+% number of iterations (default is 40) %
+if isfield(DL.param,'iternum')
+    iternum = DL.param.iternum;
+else
+    iternum = 40;
+end
+
+% show dictonary every specified number of iterations
+if isfield(DL.param,'show_dict')
+    show_dictionary=1;
+    show_iter=DL.param.show_dict;
+else
+    show_dictionary=0;
+    show_iter=0;
+end
+
+tmpTraining = Problem.b1;
+Problem.b1 = X;
+if isfield(Problem,'reconstruct')
+    Problem = rmfield(Problem, 'reconstruct');
+end
+
+
+%% Main Algorithm
+Dprev = D;						%initial dictionary
+Aprev = D\X;				    %set initial solution as pseudoinverse
+for i = 1:iternum
+	%Sparse Coding by 
+	A = sparsecoding(X,D,Aprev);
+	%Dictionary Update
+	D = dictionaryupdate(X,A,Dprev,zeta,eta);
+	
+	Dprev = D;
+	Aprev = A;
+   if ((show_dictionary)&&(mod(i,show_iter)==0))
+       dictimg = SMALL_showdict(dico,[8 8],...
+            round(sqrt(size(dico,2))),round(sqrt(size(dico,2))),'lines','highcontrast');  
+       figure(2); imagesc(dictimg);colormap(gray);axis off; axis image;
+       pause(0.02);
+   end
+end
+
+Problem.b1 = tmpTraining;
+DL.D = D;
+
+end
+
+function A = sparsecoding(X,D,Aprev)
+%Sparse coding using a mixture of laplacians (MOL) as universal prior.
+
+%parameters
+K = size(D,2);	%number of atoms
+M = size(X,2);	%number of signals
+
+mu1 = mean(abs(Aprev(:)));			%first moment of distribution of Aprev
+mu2 = (norm(Aprev(:))^2)/numel(Aprev);%second moment of distribution of Aprev
+kappa = 2*(mu2-mu1^2)/(mu2-2*mu2^2); %parameter kappa of the MOL distribution
+beta = (kappa-1)*mu1;				%parameter beta of the MOL distribution
+
+E = X-D*Aprev;						%error term
+sigmasq = mean(var(E));				%error variance
+tau = 2*sigmasq*(kappa+1);			%sparsity factor
+
+%solve a succession of subproblems to approximate the non-convex cost
+%function
+nIter = 10;							%number of iterations of surrogate subproblem
+Psi   = zeros(K,M);					%initialise solution of subproblem
+for iIter=1:nIter
+	Reg = 1./(abs(Psi) + beta);
+	Psi = solvel1(X,D,tau,Reg);
+end
+A = Psi;
+end
+
+function Psi = solvel1(X,D,tau,A)
+	[K M] = size(A);
+	Psi = zeros(K,M);
+	for m=1:M
+		cvx_begin quiet
+			variable v(K)
+			minimise (norm(X(:,m)-D*v) + tau*norm(A(:,m).*v,1));
+		cvx_end
+		Psi(:,m) = v;
+	end
+end
+
+function D = dictionaryupdate(X,A,Dprev,zeta,eta)
+	D = (X*A' + 2*(zeta + eta)*Dprev)/(A*A' + 2*zeta*(Dprev'*Dprev) + 2*eta*diag(diag(Dprev'*Dprev)));
+end
+
+
+
+function Y = colnorms_squared(X)
+% compute in blocks to conserve memory
+Y = zeros(1,size(X,2));
+blocksize = 2000;
+for i = 1:blocksize:size(X,2)
+	blockids = i : min(i+blocksize-1,size(X,2));
+	Y(blockids) = sum(X(:,blockids).^2);
+end
+end
+
+function testdl_ramirez
+	clc
+	N = 10;		%ambient dimension
+	K = 20;		%number of atoms
+	M = 30;		%number of observed signals
+	X = randn(N,M);		%observed signals
+	D = normcol(randn(N,K)); %initial dictionary
+	Problem.b = X;		%sparse representation problem
+	Problem.b1 = X;
+	DL = SMALL_init_DL('dl_ramirez');
+	DL.param.initdict = D;
+	DL.param = struct('initdict',D,...
+					  'zeta',0.5,...
+					  'eta',0.5);
+	DL = SMALL_learn(Problem,DL);
+end