annotate util/classes/dictionaryMatrices/rotatematrix.m @ 169:290cca7d3469 danieleb

Added dictionary decorrelation functions and test script for ICASSP paper.

function [Dhat cost W] = rotatematrix(D,Phi,method,param)
%
%
%
% REFERENCE
% M.D. Plumbley, Geometrical Methods for Non-Negative ICA: Manifolds, Lie
% Groups and Toral Subalgebra, Neurocomputing
if ~nargin, testrotatematrix; return, end


if ~exist('method','var') || isempty(method), method = 'unconstrained'; end

J = @(W) 0.5*norm(D-W*Phi,'fro');
cost = zeros(param.nIter,1);

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';
	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'
			eta = param.step;
			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
			H = -B;					% calculate direction as negative gradient
			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
			H = -G + polakRibiere(G,Gprev)*Hprev; %calculate conjugate gradient direction
			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
	end
end
Dhat = W*Phi;
end
% function C = matcomm(A,B)
% %Matrix commutator
% C = A*B-B*A;

function gamma = polakRibiere(G,Gprev)
gamma = G(:)'*(G(:)-Gprev(:))/(norm(Gprev(:))^2);
if isnan(gamma) || isinf(gamma)
	gamma = 0;
end
end

function t = searchline(J,H,W,t)
t = fminsearch(@(x) J(expm(x*H)*W),t);
end

function B = skew(A)
B = 0.5*(A - A');
end


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;

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);

figure, plot(cost)
set(gca,'XScale','log','Yscale','log')
legend({'uncons','settpestlie','linesearchlie','conjgradlie'})
grid on
xlabel('number of iterations')
ylabel('J(W)')
end