first commit

functions/ACIN_Cal.m

@@ -0,0 +1,61 @@
+function [ G , F , R ] = ACIN_Cal( W , X , G , F , Omega_G , Omega_F , Phi_G , Phi_F , Mean_F , mu , F_theo , N_iter )
+
+W2 = W.^2;
+Omega_Bar_G = ones(size(G))-Omega_G;
+Omega_Bar_F = ones(size(F))-Omega_F;
+Delta_G = G.*Omega_Bar_G;
+Delta_F = F.*Omega_Bar_F;
+m = size(F,2)-1;
+R = zeros(N_iter+1,1);
+R(1) = norm(F(2,1:end-1)-F_theo(2,1:end-1),2)/sqrt(size(F_theo,2)-1);
+
+for i = 1 : N_iter
+    
+    % updating G
+    Delta_G = Updt_Delta_G( W2 , X , Delta_G , Phi_G , F , Omega_Bar_G );
+    G = Phi_G + Delta_G;
+    
+    % updating F
+    Delta_F = Updt_Delta_F( W2 , X , G , Delta_F , Phi_F , Omega_Bar_F , Mean_F , mu , m );
+    F = Phi_F + Delta_F;
+    
+    R(i+1) = norm(F(2,1:end-1)-F_theo(2,1:end-1),2)/sqrt(size(F_theo,2)-1);
+    
+end
+
+end
+
+
+
+function [ Delta_F ] = Updt_Delta_F( W , X , G , Delta_F , Phi_F , Omega_B_F , Mean_F , mu , m )
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Update rule for Delta_F %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Delta_F = Delta_F.*(Omega_B_F.*(G'*(W.*secu_plus(X-G*Phi_F)))+mu/m*diag(Mean_F)*Omega_B_F)./(Omega_B_F.*(G'*(W.*(G*Delta_F)))+mu/m^2*diag(Delta_F*ones(m+1,1))*Omega_B_F);
+Delta_F(isnan(Delta_F))=0;
+end
+
+
+
+function [ Delta_G ] = Updt_Delta_G( W , X , Delta_G , Phi_G , F , Omega_B_G )
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Update rule for Delta_G %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Delta_G = Delta_G.*(Omega_B_G).*((W.*secu_plus(X-Phi_G*F))*F')./((W.*(Delta_G*F)*F')); % mise à jour
+Delta_G(isnan(Delta_G))=0;
+end
+
+
+
+function [toto] = secu_plus(tutu,s)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Goal: Security in the NMF procedure which project the negative data to
+% epsilon (small user-defined threshold)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if(nargin<2)
+    s=1.e-12;
+end
+toto=max(tutu,s);
+toto(isnan(tutu)) = 0;
+end
+

functions/IN_Cal.m

@@ -0,0 +1,60 @@
+function [ G , F , R ] = IN_Cal( W , X , G , F , Omega_G , Omega_F , Phi_G , Phi_F , F_theo , N_iter )
+
+W2 = W.^2;
+Omega_Bar_G = ones(size(G))-Omega_G;
+Omega_Bar_F = ones(size(F))-Omega_F;
+Delta_G = G.*Omega_Bar_G;
+Delta_F = F.*Omega_Bar_F;
+
+R = zeros(N_iter+1,1);
+R(1) = norm(F(2,1:end-1)-F_theo(2,1:end-1),2)/sqrt(size(F_theo,2)-1);
+
+for i = 1 : N_iter
+    
+    % updating G
+    Delta_G = Updt_Delta_G( W2 , X , Delta_G , Phi_G , F , Omega_Bar_G );
+    G = Phi_G + Delta_G;
+    
+    % updating F
+    Delta_F = Updt_Delta_F( W2 , X , G , Delta_F , Phi_F , Omega_Bar_F );
+    F = Phi_F + Delta_F;
+    
+    R(i+1) = norm(F(2,1:end-1)-F_theo(2,1:end-1),2)/sqrt(size(F_theo,2)-1);
+    
+end
+
+end
+
+
+
+function [ Delta_F ] = Updt_Delta_F( W , X , G , Delta_F , Phi_F , Omega_B_F )
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Update rule for Delta_F %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Delta_F = Delta_F.*(Omega_B_F).*((G'*(W.*secu_plus(X-G*Phi_F)))./(G'*(W.*(G*Delta_F))));
+Delta_F(isnan(Delta_F))=0;
+end
+
+
+
+function [ Delta_G ] = Updt_Delta_G( W , X , Delta_G , Phi_G , F , Omega_B_G )
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Update rule for Delta_G %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Delta_G = Delta_G.*(Omega_B_G).*((W.*secu_plus(X-Phi_G*F))*F')./((W.*(Delta_G*F)*F')); % mise à jour
+Delta_G(isnan(Delta_G))=0;
+end
+
+
+
+function [toto] = secu_plus(tutu,s)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Goal: Security in the NMF procedure which project the negative data to
+% epsilon (small user-defined threshold
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if(nargin<2)
+    s=1.e-12;
+end
+toto=max(tutu,s);
+toto(isnan(tutu)) = 0;
+end

functions/OMP.m

@@ -0,0 +1,26 @@
+function [ x  ] = OMP( y , D  , k )
+
+residual = y; 
+
+Atoms_sel = zeros(k,1);
+D_select = zeros(size(D,1),k);
+x_OMP = zeros(k,1);
+
+for iter = 1 : k
+    
+   [ ~ , idx_max ] = max(abs(residual'*D));
+    
+    Atoms_sel(iter) = idx_max;
+    D_select(:,iter) = D(:,Atoms_sel(iter));
+    
+    for j = 1 : iter-1
+        D_select(:,iter) = D_select(:,iter)-real(D_select(:,iter)'*D_select(:,j))*D_select(:,j);
+    end
+    
+    D_select(:,iter) = D_select(:,iter)/norm(D_select(:,iter));
+    x_OMP(iter) = real(D_select(:,iter)'*residual);
+    residual = residual-D_select(:,iter)*x_OMP(iter);
+    
+end
+x = zeros(size(D,2),1);
+x(Atoms_sel) = real(D(:,Atoms_sel)\y);

functions/SpAIN_Cal.m

@@ -0,0 +1,68 @@
+function [ G , F , R ] = SpAIN_Cal( W , X , G , F , Omega_G , Omega_F , Phi_G , Phi_F , D , k , lambda , Mean_F , mu , F_theo , N_iter )
+
+W2 = W.^2;
+tilde_W2 = [W2,lambda*ones(size(W,1),1)];
+
+Omega_Bar_G = ones(size(G))-Omega_G;
+Omega_Bar_F = ones(size(F))-Omega_F;
+Delta_G = G.*Omega_Bar_G;
+Delta_F = F.*Omega_Bar_F;
+m = size(F,2)-1;
+
+R = zeros(N_iter+1,1);
+R(1) = norm(F(2,1:end-1)-F_theo(2,1:end-1),2)/sqrt(size(F_theo,2)-1);
+
+for i = 1 : N_iter
+    
+    x = OMP(G(:,2),D,k);
+    tilde_X = [X,D*x];
+    tilde_F = [F,[0;1]];
+    
+    % updating G
+    Delta_G = Updt_Delta_G( tilde_W2 , tilde_X , Delta_G , Phi_G , tilde_F , Omega_Bar_G );
+    G = Phi_G + Delta_G;
+    
+    % updating F
+    Delta_F = Updt_Delta_F( W2 , X , G , Delta_F , Phi_F , Omega_Bar_F , Mean_F , mu , m );
+    F = Phi_F + Delta_F;
+    
+    R(i+1) = norm(F(2,1:end-1)-F_theo(2,1:end-1),2)/sqrt(size(F_theo,2)-1);
+    
+end
+
+end
+
+
+
+function [ Delta_F ] = Updt_Delta_F( W , X , G , Delta_F , Phi_F , Omega_B_F , Mean_F , mu , m )
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Update rule for Delta_F %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Delta_F = Delta_F.*(Omega_B_F.*(G'*(W.*secu_plus(X-G*Phi_F)))+mu/m*diag(Mean_F)*Omega_B_F)./(Omega_B_F.*(G'*(W.*(G*Delta_F)))+mu/m^2*diag(Delta_F*ones(m+1,1))*Omega_B_F);
+Delta_F(isnan(Delta_F))=0;
+end
+
+
+
+function [ Delta_G ] = Updt_Delta_G( W , X , Delta_G , Phi_G , F , Omega_B_G )
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Update rule for Delta_G %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Delta_G = Delta_G.*(Omega_B_G).*((W.*secu_plus(X-Phi_G*F))*F')./((W.*(Delta_G*F)*F')); % mise à jour
+Delta_G(isnan(Delta_G))=0;
+end
+
+
+
+function [toto] = secu_plus(tutu,s)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Goal: Security in the NMF procedure which project the negative data to
+% epsilon (small user-defined threshold)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if(nargin<2)
+    s=1.e-12;
+end
+toto=max(tutu,s);
+toto(isnan(tutu)) = 0;
+end
+

functions/SpIN_Cal.m

@@ -0,0 +1,66 @@
+function [ G , F , R ] = SpIN_Cal( W , X , G , F , Omega_G , Omega_F , Phi_G , Phi_F , D , k , lambda , F_theo , N_iter )
+
+W2 = W.^2;
+tilde_W2 = [W2,lambda*ones(size(W,1),1)];
+
+Omega_Bar_G = ones(size(G))-Omega_G;
+Omega_Bar_F = ones(size(F))-Omega_F;
+Delta_G = G.*Omega_Bar_G;
+Delta_F = F.*Omega_Bar_F;
+
+R = zeros(N_iter+1,1);
+R(1) = norm(F(2,1:end-1)-F_theo(2,1:end-1),2)/sqrt(size(F_theo,2)-1);
+
+for i = 1 : N_iter
+    
+    x = OMP(G(:,2),D,k);
+    tilde_X = [X,D*x];
+    tilde_F = [F,[0;1]];
+    
+    % updating G
+    Delta_G = Updt_Delta_G( tilde_W2 , tilde_X , Delta_G , Phi_G , tilde_F , Omega_Bar_G );
+    G = Phi_G + Delta_G;
+    
+    % updating F
+    Delta_F = Updt_Delta_F( W2 , X , G , Delta_F , Phi_F , Omega_Bar_F );
+    F = Phi_F + Delta_F;
+    
+    R(i+1) = norm(F(2,1:end-1)-F_theo(2,1:end-1),2)/sqrt(size(F_theo,2)-1);
+    
+end
+
+end
+
+
+
+function [ Delta_F ] = Updt_Delta_F( W , X , G , Delta_F , Phi_F , Omega_B_F )
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Update rule for Delta_F %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Delta_F = Delta_F.*(Omega_B_F).*((G'*(W.*secu_plus(X-G*Phi_F)))./(G'*(W.*(G*Delta_F))));
+Delta_F(isnan(Delta_F))=0;
+end
+
+
+
+function [ Delta_G ] = Updt_Delta_G( W , X , Delta_G , Phi_G , F , Omega_B_G )
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Update rule for Delta_G %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Delta_G = Delta_G.*(Omega_B_G).*((W.*secu_plus(X-Phi_G*F))*F')./((W.*(Delta_G*F)*F')); % mise à jour
+Delta_G(isnan(Delta_G))=0;
+end
+
+
+
+function [toto] = secu_plus(tutu,s)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Goal: Security in the NMF procedure which project the negative data to
+% epsilon (small user-defined threshold
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if(nargin<2)
+    s=1.e-12;
+end
+toto=max(tutu,s);
+toto(isnan(tutu)) = 0;
+end

script_demo.m

@@ -0,0 +1,151 @@
+% This matlab script provides a demonstration of the informed and
+% constrained NMF algorithms proposed in the article:
+
+% [1] C. DORFFER, M. Puigt, G. Delmaire and G. Roussel, "Informed Nonnegative
+% Matrix Factorization Methods for Mobile Sensor Network Calibration", in
+% IEEE Transactions on Signal and Information Processing over Networks,
+% vol. 4, no. 4, pp. 667-682, Dec. 2018.
+
+
+
+%%
+
+clear
+close all
+clc
+
+addpath('functions/')
+rng(3)
+
+
+%% Simulation parameters
+
+N_Ref = 3; % Nb. of reference measurements
+N_Cpt = 25; % Nb. of mobile sensors
+Mu_beta = .9; % Mean sensors gain
+Mu_alpha = 5; % Mean sensors offset
+Bound_beta = [.01;1.5]; % Gain boundaries
+Bound_alpha = [3.5;6.5]; % Offset boundaries
+
+MV = .9; % Missing Value prop.
+RV = 0; % RendezVous prop.
+var_n = 0; % Noise variance
+
+
+%% Scene simulation
+
+[xx,yy] = meshgrid((-1:2/9:1),(-1:2/9:1));
+sig = [1;1];
+y=@(xx,yy) exp(-[xx;yy]'*diag(sig)*[xx;yy]);
+for i = 1:100
+    g(i) = y(xx(i),yy(i));
+end
+g = g-min(g);
+g = .5*(g/max(g))+1e-5;
+G_theo = [ones(100,1),g']; % Theoretical matrix G (see eq.(3) of [1])
+
+figure
+subplot(221)
+imagesc(reshape(g,10,10))
+
+
+%% Sensors simulation
+
+F_theo = [max(Bound_beta(1),min(Bound_beta(2),Mu_beta+.5*randn(1,N_Cpt)));...
+     max(Bound_alpha(1),min(Bound_alpha(2),Mu_alpha+.5*randn(1,N_Cpt)))];
+F_theo = [F_theo,[0;1]]; % Theoretical matrix F (see eq.(3) of [1])
+
+
+%% Data simulation
+
+X_theo = G_theo*F_theo; %  Theoretical matrix X (see eq.(2) of [1])
+
+W = zeros(100,N_Cpt+1);
+
+idx_Ref = randperm(100);
+idx_Ref = idx_Ref(1:N_Ref); % Reference measurement locations
+W(idx_Ref,end) = 1;
+
+N_RV = round(N_Cpt*RV); % Nb. of sensors having a RendezVous
+idx_CptRV = randperm(N_Cpt);
+idx_CptRV = idx_CptRV(1:N_RV); % Selection of sensors having a RendezVous
+idx_RefRV = randi(N_Ref,1,N_Cpt);
+idx_RefRV = idx_Ref(idx_RefRV(1:N_RV)); % Selection of the references for each RendezVous
+
+for i = 1 : N_RV
+    W(idx_RefRV(i),idx_CptRV(i)) = 1;
+end
+
+N_data = round((1-MV)*(N_Cpt)*(100-N_Ref)); % Nb. of measurements in data matrix X
+xCpt = 1 : 100;
+xCpt(idx_Ref) = []; % Reference free locations
+[xx,yy] = meshgrid(xCpt,1:N_Cpt); % Possibly sensed locations
+
+idx_data = randperm((100-N_Ref)*N_Cpt);
+for i = 1 : N_data
+    W(xx(idx_data(i)),yy(idx_data(i))) = 1; % Sensor measurement placement
+end
+
+N = var_n*randn(100,N_Cpt+1); % Noise simulation
+N(:,end) = 0;
+N = max(N,-X_theo);
+SNR = snr(X_theo(W~=0),N(W~=0));
+
+X = W.*(X_theo+N); % Data matrix X
+subplot(222)
+imagesc(X)
+
+
+%% Calibration parameters
+
+% % Common parameters
+N_iter = 1.e5; % Maximum nb. of iterations
+Omega_G = [ones(100,1),W(:,end)]; % Mask on known values in G (see eq.(14) of [1])
+Omega_F = [zeros(2,N_Cpt),[1;1]]; % Mask on known values in F (see eq.(15) of [1])
+Phi_G = [ones(100,1),X(:,end)]; % Known values in G (see eq.(14) of [1])
+Phi_F = [zeros(2,N_Cpt),[0;1]]; % Known values in F (see eq.(15) of [1])
+Ginit = abs(randn(100,2)); % Initial matrix G
+Ginit = (1-Omega_G).*Ginit+Phi_G;
+Finit = [max(Bound_beta(1),min(Bound_beta(2),Mu_beta+.5*randn(1,N_Cpt)));...
+     max(Bound_alpha(1),min(Bound_alpha(2),Mu_alpha+.5*randn(1,N_Cpt)))];
+Finit = [Finit,[0;1]]; % Initial matrix F
+Finit = (1-Omega_F).*Finit+Phi_F;
+
+% % Parameters for the "average constrained" approach (called ACIN-Cal in [1])
+Mean_F = mean(F_theo(:,1:N_Cpt),2);
+mu = 10; % Regularization weight
+
+% % Parameters for the Sparsity based regularization (called SpIN-Cal in [1])
+% Dictionary construction
+D = real(ifft(diag(ones(100,1))));
+D = [D(:,1:15),g',D(:,1:15),g'];
+D(51:end,1:16) = 0;
+D(1:50,17:end) = 0;
+D = D*diag(1./sqrt(sum(D.^2))); % Atoms normalisation
+
+k = 2; % Nb. of atoms to be choosen by the OMP
+lambda = 10; % Regularization weight
+
+
+%% Calibration
+
+% % IN-Cal
+[ G_IN_Cal , F_IN_Cal , RMSE_IN_Cal ] = IN_Cal( W , X , Ginit , Finit , Omega_G , Omega_F , Phi_G , Phi_F , F_theo , N_iter );
+
+% % ACIN-Cal
+[ G_ACIN_Cal , F_ACIN_Cal , RMSE_ACIN_Cal ] = ACIN_Cal( W , X , Ginit , Finit , Omega_G , Omega_F , Phi_G , Phi_F , Mean_F , mu , F_theo , N_iter );
+
+% % SpIN-Cal
+[ G_SpIN_Cal , F_SpIN_Cal , RMSE_SpIN_Cal ] = SpIN_Cal( W , X , Ginit , Finit , Omega_G , Omega_F , Phi_G , Phi_F , D , k , lambda , F_theo , N_iter );
+
+% % SpAIN-Cal
+[ G_SpAIN_Cal , F_SpAIN_Cal , RMSE_SpAIN_Cal ] = SpAIN_Cal( W , X , Ginit , Finit , Omega_G , Omega_F , Phi_G , Phi_F , D , k , lambda , Mean_F , mu , F_theo , N_iter );
+
+
+subplot(223)
+semilogy(RMSE_IN_Cal)
+axis([0 N_iter 1.e-16 1])
+hold all
+semilogy(RMSE_ACIN_Cal)
+semilogy(RMSE_SpIN_Cal)
+semilogy(RMSE_SpAIN_Cal)