daniele@169: function [Dhat cost W] = rotatematrix(D,Phi,method,param)
daniele@169: %
daniele@169: %
daniele@169: %
daniele@169: % REFERENCE
daniele@169: % M.D. Plumbley, Geometrical Methods for Non-Negative ICA: Manifolds, Lie
daniele@169: % Groups and Toral Subalgebra, Neurocomputing
daniele@170: 
daniele@170: %% Parse inputs and set defaults
daniele@169: if ~nargin, testrotatematrix; return, end
daniele@169: 
daniele@170: if ~exist('param','var')  || isempty(param),  param  = struct; end
daniele@170: if ~exist('method','var') || isempty(method), method = 'conjgradLie'; end
daniele@170: if ~isfield(param,'nIter'), param.nIter = 100; end	%number of iterations
daniele@170: if ~isfield(param,'eps'), param.eps = 1e-9; end		%tolerance level
daniele@170: if ~isfield(param,'step'), param.step = 0.01; end
daniele@169: 
daniele@170: J = @(W) 0.5*norm(D-W*Phi,'fro');		%cost function
daniele@169: 
daniele@170: % Initialise variables
daniele@170: cost = zeros(param.nIter,1);			%cost at each iteration
daniele@170: W	 = eye(size(Phi,1));				%rotation matrix
daniele@170: grad = ones(size(W));					%gradient
daniele@170: t	 = param.step;						%step size
daniele@170: Gprev = 0;								%previous gradient
daniele@170: Hprev = 0;								%previous Lie search direction
daniele@170: iIter = 1;								%iteration counter
daniele@169: 
daniele@170: %% Main algorithm
daniele@170: while iIter<=param.nIter && norm(grad,'fro')>eps
daniele@170: 	cost(iIter) = J(W);				%calculate cost
daniele@170: 	grad = (W*Phi-D)*Phi';			%calculate gradient
daniele@169: 	switch method
daniele@169: 		case 'unconstrained'		% gradient descent
daniele@169: 			eta = param.step;
daniele@169: 			W = W - eta*grad;		% update W by steepest descent
daniele@169: 		case 'tangent'				% self correcting tangent
daniele@169: 			eta = param.step;
daniele@170: 			W = W - 0.5*eta*(grad - W*grad'*W);
daniele@170: 			[U , ~, V] = svd(W);
daniele@170: 			W = U*V';
daniele@170: 		case 'steepestlie'			%steepest descent in Lie algebra
daniele@169: 			eta = param.step;
daniele@170: 			B = 2*skew(grad*W');	% calculate gradient in Lie algebra
daniele@169: 			W = expm(-eta*B)*W;		% update W by steepest descent
daniele@170: 		case 'linesearchlie'		% line search in Lie algebra
daniele@170: 			B = 2*skew(grad*W');	% calculate gradient in Lie algebra
daniele@169: 			H = -B;					% calculate direction as negative gradient
daniele@170: 			t = searchline(J,H,W,t);% line search in one-parameter Lie subalgebra
daniele@169: 			W = expm(t*H)*W;		% update W by line search
daniele@170: 		case 'conjgradlie'			% conjugate gradient in Lie algebra
daniele@170: 			G = 2*skew(grad*W');	% calculate gradient in Lie algebra
daniele@169: 			H = -G + polakRibiere(G,Gprev)*Hprev; %calculate conjugate gradient direction
daniele@170: 			t = searchline(J,H,W,t);% line search in one-parameter Lie subalgebra
daniele@169: 			W = expm(t*H)*W;		% update W by line search
daniele@170: 			Hprev = H;				% save search direction
daniele@170: 			Gprev = G;				% save gradient
daniele@169: 	end
daniele@170: 	iIter = iIter+1;				% update iteration counter
daniele@169: end
daniele@170: Dhat = W*Phi;						%rotate matrix
daniele@172: cost(iIter:end) = cost(iIter-1);	%zero-pad cost vector
daniele@169: end
daniele@169: 
daniele@170: %% Support functions
daniele@169: function gamma = polakRibiere(G,Gprev)
daniele@170: %Polak-Ribiere rule for conjugate direction calculation
daniele@169: gamma = G(:)'*(G(:)-Gprev(:))/(norm(Gprev(:))^2);
daniele@169: if isnan(gamma) || isinf(gamma)
daniele@169: 	gamma = 0;
daniele@169: end
daniele@169: end
daniele@169: 
daniele@169: function t = searchline(J,H,W,t)
daniele@170: %Line search in one-parameter Lie subalgebra
daniele@169: t = fminsearch(@(x) J(expm(x*H)*W),t);
daniele@169: end
daniele@169: 
daniele@169: function B = skew(A)
daniele@170: %Skew-symmetric matrix
daniele@169: B = 0.5*(A - A');
daniele@169: end
daniele@169: 
daniele@169: 
daniele@170: %% Test function
daniele@169: function testrotatematrix
daniele@169: clear, clc, close all
daniele@170: n = 256;							%ambient dimension
daniele@170: m = 512;							%number of atoms
daniele@170: param.nIter = 300;				%number of iterations
daniele@170: param.step  = 0.001;			%step size
daniele@170: param.mu    = 0.01;			%regularization factor (for tangent method)
daniele@170: methods = {'unconstrained','tangent','linesearchlie','conjgradlie'};
daniele@169: 
daniele@170: Phi = randn(n,m);				%initial dictionary
daniele@170: Qtrue = expm(skew(randn(n)));	%rotation matrix
daniele@170: D = Qtrue*Phi;					%target dictionary
daniele@170: 
daniele@170: cost = zeros(param.nIter,length(methods));
daniele@170: for iIter=1:length(methods)
daniele@170: 	tic
daniele@170: 	[~, cost(:,iIter)] = rotatematrix(D,Phi,methods{iIter},param);
daniele@170: 	time = toc;
daniele@170: 	sprintf('Method %s completed in %f seconds \n',methods{iIter},time)
daniele@170: end
daniele@169: 
daniele@169: figure, plot(cost)
daniele@169: set(gca,'XScale','log','Yscale','log')
daniele@170: legend(methods)
daniele@169: grid on
daniele@169: xlabel('number of iterations')
daniele@169: ylabel('J(W)')
daniele@169: end