Mercurial > hg > smallbox
diff util/classes/dictionaryMatrices/rotatematrix.m @ 170:68fb71aa5339 danieleb
Added dictionary decorrelation functions and test script for Letters paper.
author | Daniele Barchiesi <daniele.barchiesi@eecs.qmul.ac.uk> |
---|---|
date | Thu, 06 Oct 2011 14:33:41 +0100 |
parents | 290cca7d3469 |
children | 9c41f87dead7 |
line wrap: on
line diff
--- a/util/classes/dictionaryMatrices/rotatematrix.m Thu Sep 29 09:46:52 2011 +0100 +++ b/util/classes/dictionaryMatrices/rotatematrix.m Thu Oct 06 14:33:41 2011 +0100 @@ -5,54 +5,66 @@ % REFERENCE % M.D. Plumbley, Geometrical Methods for Non-Negative ICA: Manifolds, Lie % Groups and Toral Subalgebra, Neurocomputing + +%% Parse inputs and set defaults if ~nargin, testrotatematrix; return, end +if ~exist('param','var') || isempty(param), param = struct; end +if ~exist('method','var') || isempty(method), method = 'conjgradLie'; end +if ~isfield(param,'nIter'), param.nIter = 100; end %number of iterations +if ~isfield(param,'eps'), param.eps = 1e-9; end %tolerance level +if ~isfield(param,'step'), param.step = 0.01; end -if ~exist('method','var') || isempty(method), method = 'unconstrained'; end +J = @(W) 0.5*norm(D-W*Phi,'fro'); %cost function -J = @(W) 0.5*norm(D-W*Phi,'fro'); -cost = zeros(param.nIter,1); +% Initialise variables +cost = zeros(param.nIter,1); %cost at each iteration +W = eye(size(Phi,1)); %rotation matrix +grad = ones(size(W)); %gradient +t = param.step; %step size +Gprev = 0; %previous gradient +Hprev = 0; %previous Lie search direction +iIter = 1; %iteration counter -W = eye(size(Phi,1)); -t = 0; -Gprev = 0; -Hprev = 0; -for i=1:param.nIter - cost(i) = J(W); - grad = (W*Phi-D)*Phi'; +%% Main algorithm +while iIter<=param.nIter && norm(grad,'fro')>eps + cost(iIter) = J(W); %calculate cost + grad = (W*Phi-D)*Phi'; %calculate gradient switch method case 'unconstrained' % gradient descent eta = param.step; W = W - eta*grad; % update W by steepest descent case 'tangent' % self correcting tangent eta = param.step; - mu = param.reg; - W = W - 0.5*eta*(grad - W*grad'*W + mu*W*(W'*W-eye(size(W)))); - case 'steepestlie' + W = W - 0.5*eta*(grad - W*grad'*W); + [U , ~, V] = svd(W); + W = U*V'; + case 'steepestlie' %steepest descent in Lie algebra eta = param.step; - B = 2*skew(grad*W'); % calculate gradient in lie algebra + B = 2*skew(grad*W'); % calculate gradient in Lie algebra W = expm(-eta*B)*W; % update W by steepest descent - case 'linesearchlie' - B = 2*skew(grad*W'); % calculate gradient in lie algebra + case 'linesearchlie' % line search in Lie algebra + B = 2*skew(grad*W'); % calculate gradient in Lie algebra H = -B; % calculate direction as negative gradient - t = searchline(J,H,W,t);% line search in one-parameter lie subalgebra + t = searchline(J,H,W,t);% line search in one-parameter Lie subalgebra W = expm(t*H)*W; % update W by line search - case 'conjgradlie' - G = 2*skew(grad*W'); % calculate gradient in lie algebra + case 'conjgradlie' % conjugate gradient in Lie algebra + G = 2*skew(grad*W'); % calculate gradient in Lie algebra H = -G + polakRibiere(G,Gprev)*Hprev; %calculate conjugate gradient direction - t = searchline(J,H,W,t);% line search in one-parameter lie subalgebra + t = searchline(J,H,W,t);% line search in one-parameter Lie subalgebra W = expm(t*H)*W; % update W by line search - Hprev = H; % % save search direction - Gprev = G; % % save gradient + Hprev = H; % save search direction + Gprev = G; % save gradient end + iIter = iIter+1; % update iteration counter end -Dhat = W*Phi; +Dhat = W*Phi; %rotate matrix +cost(iIter:end) = cost(iIter-1); %pad cost vector end -% function C = matcomm(A,B) -% %Matrix commutator -% C = A*B-B*A; +%% Support functions function gamma = polakRibiere(G,Gprev) +%Polak-Ribiere rule for conjugate direction calculation gamma = G(:)'*(G(:)-Gprev(:))/(norm(Gprev(:))^2); if isnan(gamma) || isinf(gamma) gamma = 0; @@ -60,36 +72,41 @@ end function t = searchline(J,H,W,t) +%Line search in one-parameter Lie subalgebra t = fminsearch(@(x) J(expm(x*H)*W),t); end function B = skew(A) +%Skew-symmetric matrix B = 0.5*(A - A'); end +%% Test function function testrotatematrix clear, clc, close all -n = 256; -m = 512; -disp('A random matrix...'); -Phi = randn(n,m); -disp('And its rotated mate...'); -Qtrue = expm(skew(randn(n))); -D = Qtrue*Phi; -disp('Now, lets try to find the right rotation...'); -param.nIter = 1000; -param.step = 0.001; +n = 256; %ambient dimension +m = 512; %number of atoms +param.nIter = 300; %number of iterations +param.step = 0.001; %step size +param.mu = 0.01; %regularization factor (for tangent method) +methods = {'unconstrained','tangent','linesearchlie','conjgradlie'}; -cost = zeros(param.nIter,4); -[~, cost(:,1)] = rotatematrix(D,Phi,'unconstrained',param); -[~, cost(:,2)] = rotatematrix(D,Phi,'steepestlie',param); -[~, cost(:,3)] = rotatematrix(D,Phi,'linesearchlie',param); -[~, cost(:,4)] = rotatematrix(D,Phi,'conjgradlie',param); +Phi = randn(n,m); %initial dictionary +Qtrue = expm(skew(randn(n))); %rotation matrix +D = Qtrue*Phi; %target dictionary + +cost = zeros(param.nIter,length(methods)); +for iIter=1:length(methods) + tic + [~, cost(:,iIter)] = rotatematrix(D,Phi,methods{iIter},param); + time = toc; + sprintf('Method %s completed in %f seconds \n',methods{iIter},time) +end figure, plot(cost) set(gca,'XScale','log','Yscale','log') -legend({'uncons','settpestlie','linesearchlie','conjgradlie'}) +legend(methods) grid on xlabel('number of iterations') ylabel('J(W)')