Initial commit

Author: adavoudi

DPLM.m

@@ -0,0 +1,91 @@
+% This is the supervised version of the DPLM algorithm which
+% needs the labels. Please cite the following paper if you use 
+% this code:
+%		Davoudi, Alireza, Saeed Shiry Ghidary, and Khadijeh Sadatnejad. 
+%		"Dimensionality reduction based on distance preservation to local 
+%		mean for symmetric positive definite matrices and its application 
+%		in brain–computer interfaces." Journal of Neural Engineering 14.3 
+%		(2017): 036019.
+% 
+% PARAMS:
+% 	data: An m*m*n dimensional matrix which contains n SPD matrix of size m*m
+% 	dim : The dimensionality of SPD matrices after DM
+%	labels : A verctor of size n which contains the label of each point
+%	k	: Number of neighbours
+%	Adj	: The adjacency matrix calculated by `DPLM_adjmat` function
+%	M	: The means matrix calculated by `DPLM_adjmat` function
+%
+% RETUTNS:
+%	U	: The calculated transformation matrix which can be used as 
+%		  below to transform an m*m SPD matrix `x` to an dim*dim 
+%		  SPD matrix y:
+%				y = U'*x*U
+%	obj	: Final value of the objective function
+%	Adj	: The adjacency matrix calculated by `DPLM_adjmat` function
+%	M	: The means matrix calculated by `DPLM_adjmat` function
+%
+function [U, obj, Adj, M] = DPLM( data, labels, dim, k, varargin )
+
+if(~isempty(varargin))
+    Adj = varargin{1};
+    M   = varargin{2};
+    if(length(varargin) > 2)
+        verbose = varargin{3};
+    else
+        verbose = 1;
+    end
+else
+    [ Adj, M ] = DPLM_adjmat( data, labels, k );
+    verbose = 1;
+end
+
+u0 = randn(size(data, 1), dim);
+u0 = orth(u0);
+
+opts.record = verbose;
+opts.mxitr  = 1000;
+opts.xtol = .1;
+opts.gtol = .1;
+opts.ftol = .1;
+%obj.tau = 1e-3;
+%opts.nt = 1;
+
+t = tic;
+[U, obj]= OptStiefelGBB(u0, @objfunc, opts);
+tsolve = toc(t);
+
+if(verbose)
+    [f,~] = objfunc(u0);
+    [flast,~] = objfunc(U);
+    
+    fprintf('Elapsed time: %f\n', tsolve);
+    fprintf('Stoped: %s , Init f: %f, Last f: %f \n',obj.msg, f, flast);
+end
+
+
+    function [F, G] = objfunc(u)
+        
+        F = 0;
+        G = 0;
+        for i = 1:size(data, 3)
+            for j = 1:size(data, 3)
+                if(Adj(i,j) == 0)
+                    continue;
+                end
+                
+                Ft = distance_ld(data(:,:,j),M(:,:,i))...
+                    -distance_ld(u' * data(:,:,j) * u, u' * M(:,:,i) * u);
+                F = F + abs(Ft);
+                
+                G = G -sign(Ft) * (2 * (data(:,:,j) + M(:,:,i)) * u ...
+                    * inv( u' * (data(:,:,j) + M(:,:,i)) * u) ...
+                    - data(:,:,j) * u * inv(u' * data(:,:,j) * u) ...
+                    - M(:,:,i) * u * inv(u' * M(:,:,i) * u));
+                
+            end
+        end
+        
+    end
+
+end
+

DPLM_adjmat.m

@@ -0,0 +1,38 @@
+% This function calculates the adjacency matrix (Adj) and 
+% the means (M) of every point's neighbourhood (according to the 
+% adjacency matrix).
+% Note that this is an within class adjacancy matrix (i.e. the
+% value of Adj(i,j) for the two points i and j of two different 
+% classes is zero.
+% 
+% 
+% Please cite the following paper if you use 
+% this code:
+%		Davoudi, Alireza, Saeed Shiry Ghidary, and Khadijeh Sadatnejad. 
+%		"Dimensionality reduction based on distance preservation to local 
+%		mean for symmetric positive definite matrices and its application 
+%		in brain–computer interfaces." Journal of Neural Engineering 14.3 
+%		(2017): 036019.
+% 
+%
+function [ Adj, M ] = DPLM_adjmat( data, labels, k )
+
+dist = inf * ones(numel(labels));
+
+for i = 1:numel(labels)
+    for j = i+1:numel(labels)
+        if(labels(i) == labels(j))
+            dist(i,j) = distance_riemann(data(:,:,i), data(:,:,j));
+            dist(j,i) = dist(i,j);
+        end
+    end  
+end
+
+M = zeros(size(data,1),size(data,1),numel(labels));
+Adj = zeros(numel(labels));
+for i = 1:numel(labels)
+    [~,idx] = sort(dist(i,:), 'ascend');
+    Adj(i,idx(1:k)) = 1;
+    M(:,:,i) = mean_covariances(data(:,:,idx),'riemann');
+end
+

Dependencies/FOptM/MGramSchmidt.m

@@ -0,0 +1,16 @@
+function V = MGramSchmidt(V)
+[n,k] = size(V);
+
+for dj = 1:k
+    for di = 1:dj-1
+        V(:,dj) = V(:,dj) - proj(V(:,di), V(:,dj));
+    end
+    V(:,dj) = V(:,dj)/norm(V(:,dj));
+end
+end
+
+
+%project v onto u
+function v = proj(u,v)
+v = (dot(v,u)/dot(u,u))*u;
+end