Add files via upload

Uploading code accompanying ICASSP'20 submission.

Author: shekkizh

DistEuclideanPiotrDollar.m

@@ -0,0 +1,92 @@
+% Calculates the Euclidean distance between vectors [FAST].
+%
+% Assume X is an m-by-p matrix representing m points in p-dimensional space and Y is an
+% n-by-p matrix representing another set of points in the same space. This function
+% compute the m-by-n distance matrix D where D(i,j) is the SQUARED Euclidean distance
+% between X(i,:) and Y(j,:).  Running time is O(m*n*p).
+%
+% If x is a single data point, here is a faster, inline version to use:
+%   D = sum( (Y - ones(size(Y,1),1)*x).^2, 2 )';
+%
+% INPUTS
+%   X   - m-by-p matrix of m p dimensional vectors 
+%   Y   - n-by-p matrix of n p dimensional vectors 
+%
+% OUTPUTS
+%   D   - m-by-n distance matrix
+%
+% EXAMPLE
+%   X=[randn(100,5)]; Y=randn(40,5)+2;
+%   D = dist_euclidean( [X; Y], [X; Y] ); im(D)
+%
+% DATESTAMP
+%   29-Sep-2005  2:00pm
+%
+% See also DIST_CHISQUARED, DIST_EMD
+
+% Piotr's Image&Video Toolbox      Version 1.03   
+% Written and maintained by Piotr Dollar    pdollar-at-cs.ucsd.edu 
+% Please email me if you find bugs, or have suggestions or questions! 
+ 
+function D = DistEuclideanPiotrDollar( X, Y )
+    if( ~isa(X,'double') || ~isa(Y,'double'))
+%         disp( 'Inputs must be of type double - converting...');
+        X = double(X); Y=double(Y);
+    end
+    m = size(X,1); n = size(Y,1);  
+    Yt = Y';  
+    XX = sum(X.*X,2);        
+    YY = sum(Yt.*Yt,1);      
+    D = XX(:,ones(1,n)) + YY(ones(1,m),:) - 2*X*Yt;
+    
+    
+    
+
+%%%% code from Charles Elkan with variables renamed
+%    m = size(X,1); n = size(Y,1);
+%    D = sum(X.^2, 2) * ones(1,n) + ones(m,1) * sum(Y.^2, 2)' - 2.*X*Y';
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    
+    
+%%%% LOOP METHOD - SLOW
+%     [m p] = size(X);  
+%     [n p] = size(Y);
+%     
+%     D = zeros(m,n);
+%     ones_m_1 = ones(m,1);
+%     for i=1:n
+%         y = Y(i,:);
+%         d = X - y(ones_m_1,:);
+%         D(:,i) = sum( d.*d, 2 );  
+%     end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+%%% PARALLEL METHOD THAT IS SUPER SLOW (slower then loop)!
+% % Code taken from "MATLAB array manipulation tips and tricks" by Peter J. Acklam
+%     Xb = permute(X, [1 3 2]);  
+%     Yb = permute(Y, [3 1 2]);
+%     D = sum( (Xb(:,ones(1,n),:) - Yb(ones(1,m),:,:)).^2, 3);    
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
+
+
+%%%%% USELESS FOR EVEN VERY LARGE ARRAYS X=16000x1000!! and Y=100x1000
+%     % call recursively to save memory
+%     if( (m+n)*p > 10^5 && (m>1 || n>1))
+%         if( m>n )
+%             X1 = X(1:floor(end/2),:);
+%             X2 = X((floor(end/2)+1):end,:);
+%             D1 = dist_euclidean( X1, Y );
+%             D2 = dist_euclidean( X2, Y );
+%             D = cat( 1, D1, D2 );
+%         else
+%             Y1 = Y(1:floor(end/2),:);
+%             Y2 = Y((floor(end/2)+1):end,:);
+%             D1 = dist_euclidean( X, Y1 );
+%             D2 = dist_euclidean( X, Y2 );
+%             D = cat( 2, D1, D2 );
+%         end
+%         return;
+%     end 
+%        

GD_BuildDirectedKnnGraph.m

@@ -0,0 +1,67 @@
+function W = GD_BuildDirectedKnnGraph(M,k,which_matrix)
+% Usage: W = GD_BuildDirectedKnnGraph(M,k,which_matrix) 
+%
+% Input: 
+% M                = either the distance or the similarity matrix, needs to be square, symmetric, non-negative
+% k                = connectivity parameter of the kNN graph
+% which_matrix     = either 'sim' or 'dist' (similarity or distance matrix)
+% 
+% Output: 
+% W              = adjacency matrix of the directed knn graph
+%
+% For a similarity matrix S, returns the directed knn graph, edges are weighted by S. 
+% For a distance matrix D, returns the directed (unweighted!) knn graph. If you want to get 
+% a weighted graph in this case, you need to take care of transforming D to S yourself and then 
+% call the function with a similarity matrix.  
+% Self-edges are excluded in both graphs. 
+
+% implemented by brute force sorting
+
+% % testing whether matrix is square: 
+% if (size(M,1) ~= size(M,2))
+%   error('Matrix not square!')
+% end
+[n, m] = size(M); 
+
+% to exclude self-edges, set diagonal of sim/dissim matrix to Inf or 0
+for it=1:min(n,m)
+  if (strcmp(which_matrix,'sim'))
+    M(it,it) = 0; 
+  elseif (strcmp(which_matrix, 'dist'))
+    M(it,it) = Inf; 
+  else
+    error('Unknown matrix type')
+  end
+end
+
+
+% now do it: 
+W = M; 
+
+if (strcmp(which_matrix, 'sim'))
+  
+  for it = 1:n
+    % sort points according to similarities: 
+    [sorted_row,order] = sort(M(it,:), 'descend'); 
+    
+    % for all points which are not among the k nearest neighbors, set W to 0: 
+    W(it, order(k+1:end)) = 0;
+    %W(it,order(1:k) just stays the same
+  end
+  
+elseif (strcmp(which_matrix, 'dist'))
+  
+  for it = 1:n
+    % sort points according to distances: 
+    [sorted_row,order] = sort(M(it,:), 'ascend'); 
+    % for all points which are not among the k nearest neighbors, set W to 0: 
+    W(it, order(k+1:end)) = 0;
+    W(it, order(1:k)) = 1; %unweighted! 
+  end
+  
+else 
+  error('build_directed_knn_graph: unknown matrix type')
+end
+
+
+  

GD_BuildSymmetricKnnGraph.m

@@ -0,0 +1,31 @@
+function W = GD_BuildSymmetricKnnGraph(M,k,which_matrix)
+% Usage:  W = GD_BuildSymmetricKnnGraph(M,k,which_matrix)
+%
+% Input: 
+% M                = either the distance or the similarity matrix, needs to be square
+% k                = parameter of the knn graph
+% which_matrix     = either 'sim' or 'dist'
+% 
+% Output: 
+% W              = adjacency matrix of the symmetric knn graph without self-edges
+%                  
+%
+% For a similarity matrix S, returns the undirected knn graph, edges are weighted by S. 
+% For a distance matrix D, returns the undirected (unweighted!) knn graph. If you want to get 
+% a weighted graph in this case, you need to take care of transforming D to S yourself and then 
+% call the function with a similarity matrix 
+%
+% Excludes self-edges. 
+
+%implemented brute force 
+
+
+
+if (size(M,1) ~= size(M,2))
+  error('Matrix not square!')
+end
+
+% compute directed KNN graph: 
+W = GD_BuildDirectedKnnGraph(M,k+1,which_matrix); 
+% transform it to symmetric one: 
+W = max(W,W'); 

compute_error.m

@@ -0,0 +1,7 @@
+function err = compute_error(labels,labels_hat,samples)
+% Error ratio on unlabeled nodes
+n = length(labels);
+unlabeled = true(n,1);
+unlabeled(samples) = false; % true => node is unlabeled
+err = sum(labels(unlabeled) ~= labels_hat(unlabeled))/sum(unlabeled);
+

compute_laplacians.m

@@ -0,0 +1,31 @@
+function [L, Ln, Lr] = compute_laplacians(A)
+% Computes L, symm. norm. L and RW L
+
+N = length(A);
+% Remove self edges if any
+A = A - diag(diag(A));
+deg = sum(A,2);
+L = spdiags(deg,0,N,N) - A;
+
+d = deg;
+d = d.^(-1/2);
+d(deg==0) = 0;
+Dinv = spdiags(d,0,N,N);
+Ln = speye(N)-Dinv*A*Dinv;
+Ln = 0.5*(Ln+Ln');
+
+d = deg;
+d = d.^(-1);
+d(deg==0) = 0;
+Lr = speye(N)-spdiags(d,0,N,N)* A;
+
+
+% % normalization by number of links
+% % remove self loops
+% A =  A - diag(diag(A));
+% d_link = sum(A~=0, 2);
+% Dinv_link = spdiags(d_link.^(-1/2), 0, N, N);
+% Ln_link = Dinv_link*Lun*Dinv_link;
+% Ln_link = 0.5*(Ln_link+Ln_link');
+
+

gfhf_recon_multiclass.m

@@ -0,0 +1,22 @@
+function [err,labels_hat] = gfhf_recon_multiclass(L,labeled,y)
+% SSL using Zhu et al. 2003
+% L = Laplacian (or any other HP function of L e.g. L^k)
+% labeled = indices of labeled nodes
+% y = ground truth labels (one vs rest 1-0 membership)
+% (Assumes non-empty labeled set)
+
+[~,labels] = max(y,[],2); % true labels
+unlabeled = setdiff(1:1:size(y,1),labeled);
+y_hat = y;
+full_matrix = full(L(unlabeled,unlabeled));
+nn = any(isnan(full_matrix(:)));
+in = any(isinf(full_matrix(:)));
+if (nn || in)
+    disp('NaN or Inf detected')
+end
+
+y_hat(unlabeled,:) = -pseudo_inverse(L(unlabeled,unlabeled))*(L(unlabeled,labeled)*y(labeled,:));
+labels_hat = labels;
+[~,labels_hat(unlabeled)] = max(y_hat(unlabeled,:),[],2); % unknown predicted labels
+err = compute_error(labels,labels_hat,labeled);
+end

nnk_graph_demo.m

@@ -0,0 +1,50 @@
+clc;
+clear;
+close all;
+%% read data
+knn_param=10; k_choice=10; reg=1e-6;
+data = 'two_moons_1k';
+load([data,'.mat']); % data: X (dim x num) and labels: y
+
+N = size(X,2);
+results_folder = ['results/']; 
+dir_result = mkdir(results_folder);
+%% compute graphs
+tic
+D = DistEuclideanPiotrDollar(X',X'); % pairwise squared Euclidean distances
+directed_knn_mask = sparse(GD_BuildDirectedKnnGraph(D,knn_param,'dist'));
+distance_mask_time = toc;
+kD = sort(sqrt(D), 'ascend');
+sigma = full(mean(kD(k_choice, :)))/3;
+G = exp(-D./(2*sigma*sigma));
+similarity_time = toc;
+knn_mask = max(directed_knn_mask, directed_knn_mask');
+symmetrization_time = toc;
+%%
+fprintf('Computing the adj and L of %d-NN graph with sigma %0.4f...\n', knn_param, sigma)
+W_knn = G .* knn_mask;
+W_knn(W_knn<reg) = 0;
+knn_time = toc;
+
+%%
+fprintf('Computing the adj and L NNK...\n')
+tic
+W_nnk = nnk_inverse_kernel_graph(G, directed_knn_mask, knn_param, reg); % choose the min k-NN sim
+nnk_time = toc + similarity_time;
+
+%%
+time_values = {knn_time, nnk_time};
+sparsity_values = {length(find(W_knn))/2, length(find(W_nnk))/2};
+%%
+fname = [results_folder, data,'_k_',num2str(knn_param),'_sig_', num2str(round(sigma))]; %
+save([fname, '.mat'], 'knn_param', 'sigma', 'W_knn', ...%
+    'W_nnk', 'time_values', 'sparsity_values'); 
+%% Adjacency plots
+figure();
+subplot(1,2,1); axis off
+spy(W_knn);
+title(['KNN (t=' num2str(time_values{1}) ', edges=' num2str(sparsity_values{1}) ')'])
+
+subplot(1,2,2); axis off
+spy(W_nnk);
+title(['NNK (t=' num2str(time_values{2}) ', edges=' num2str(sparsity_values{2}) ')'])

nnk_inverse_kernel_graph.m

@@ -0,0 +1,38 @@
+function W = nnk_inverse_kernel_graph(G, mask, knn_param, reg)
+% Constructs NNK graph
+% G = similairty or kernel matrix
+% mask - Neighbors to consider at each node (each row corresponds to neighbors)
+% k_choice - Not used. For adaptive kernels
+% reg - for removing weights smaller than reg value
+
+if nargin < 4
+    reg=1e-6;
+end
+n = length(G);
+neighbor_indices = zeros(n, knn_param);
+weight_values = zeros(n, knn_param);
+error_values = zeros(n,knn_param);
+
+for i = 1:n
+    nodes_i = find(mask(i,:)); % only consider similar enough neighbors
+    nodes_i(nodes_i==i) = []; 
+    G_i = full(G(nodes_i, nodes_i));% removing i-th row and column from G
+    g_i = full(G(nodes_i,i)); % removing the i-th elem from 
+
+    qp_output = nonnegative_qp_solver(G_i, g_i, reg, g_i);
+    qpsol = qp_output.xopt;
+    
+    neighbor_indices(i,:) = nodes_i;
+    weight_values(i,:) = qpsol;
+    error_values(i,:) = (G(i,i) - 2*qpsol'*g_i + qpsol'*G_i*qpsol);
+end
+
+%%
+row_indices = repmat(1:n, knn_param,1)';
+W = sparse(row_indices(:), neighbor_indices(:), weight_values(:), n, n);
+error = sparse(row_indices(:), neighbor_indices(:), error_values(:), n, n);
+W(error > error') = 0;
+W = max(W, W');
+
+
+