remnenmf_full.m 3.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160
  1. function [ T , RMSE ] = remnenmf_full( W , X , G , F , Omega_G, Omega_F, Phi_G, Phi_F , InnerMinIter , InnerMaxIter , Tmax , v, F_theo, delta_measure, nu)
  2. r = nu;
  3. X0 = X;
  4. Omega_G = (Omega_G == 1); % Logical mask is faster than indexing in matlab.
  5. Omega_F = (Omega_F == 1); % Logical mask is faster than indexing in matlab.
  6. nOmega_G = ~Omega_G; % Logical mask is faster than indexing in matlab.
  7. nOmega_F = ~Omega_F; % Logical mask is faster than indexing in matlab.
  8. [~, num_sensor] = size(F);
  9. num_sensor = num_sensor-1;
  10. em_iter_max = round(Tmax / delta_measure) ;
  11. T = nan(1,em_iter_max);
  12. RMSE = nan(2,em_iter_max);
  13. nW = (1-W);
  14. % X = G*F+W.*(X0-G*F);
  15. X = X0 + nW.*(G*F);
  16. GG = G' * G;
  17. GX = G' * X ;
  18. GradF = GG * F - GX;
  19. FF = F * F';
  20. XF = X * F' ;
  21. GradG = nOmega_G.*(G * FF - XF);
  22. d = Grad_P([GradG',GradF],[G',F]);
  23. StoppingCritF = 1.e-3*d;
  24. StoppingCritG = 1.e-3*d;
  25. T_E = [];
  26. T_M = [];
  27. tic
  28. i = 1;
  29. niter = 0;
  30. RMSE(:,i) = vecnorm(F(:,1:end-1)- F_theo(:,1:end-1),2,2)/sqrt(num_sensor);
  31. T(i) = toc;
  32. while toc<Tmax
  33. t_e = toc;
  34. % Estimation step
  35. X = X0 + nW.*(G*F);
  36. % if i>1
  37. % [L,R]=RSI_compression(X,r,L,R);
  38. % else
  39. % [L,R]=RSI_compression(X,r);
  40. % end
  41. [L,R]=RSI_compression(X,r);
  42. % Compress left and right
  43. X_L = L * X;
  44. X_R = X * R;
  45. T_E = cat(1,T_E,toc - t_e);
  46. % Maximization step
  47. for j =1:v
  48. t_m = toc;
  49. % F_R = F * R;
  50. % FF = F_R * F_R';
  51. FF = F * F';
  52. XF = X_R * (F * R)' - Phi_G * FF;
  53. % G(Omega_G) = 0; % Convert G to \Delta G
  54. [ GradG , iterG ] = MaJ_G_EM_NeNMF( FF , XF , GradG , InnerMinIter , InnerMaxIter , StoppingCritG , nOmega_G); % Update \Delta G
  55. % G(Omega_G) = Phi_G(Omega_G); % Convert \Delta G to G
  56. G = GradG + Phi_G;
  57. niter = niter + iterG;
  58. if(iterG<=InnerMinIter)
  59. StoppingCritG = 1.e-1*StoppingCritG;
  60. end
  61. % G_L = L * G;
  62. % GG = G_L' * G_L;
  63. GG = G' * G;
  64. GX = (L * G)' * X_L - GG * Phi_F;
  65. F(Omega_F) = 0; % Convert F to \Delta F
  66. % F = F - Phi_F;
  67. [ F , iterF ] = MaJ_F_EM_NeNMF( GG , GX , F , InnerMinIter , InnerMaxIter , StoppingCritF , nOmega_F); % Update \Delta F
  68. F(Omega_F) = Phi_F(Omega_F); % Convert \Delta F to F
  69. % F = F + Phi_F;
  70. niter = niter + iterF;
  71. if(iterF<=InnerMinIter)
  72. StoppingCritF = 1.e-1*StoppingCritF;
  73. end
  74. if toc - i*delta_measure >= delta_measure
  75. i = i+1;
  76. if i > em_iter_max
  77. break
  78. end
  79. T(i) = toc;
  80. RMSE(:,i) = vecnorm(F(:,1:end-1) - F_theo(:,1:end-1),2,2)/sqrt(num_sensor);
  81. end
  82. T_M = cat(1,T_M,toc - t_m);
  83. end
  84. end
  85. niter
  86. disp(['rem E step : ',num2str(median(T_E))])
  87. disp(['rem M step : ',num2str(median(T_M))])
  88. end
  89. function [ L,R ] = RSI_compression(X,r,varargin)
  90. % Tepper, M., & Sapiro, G. (2016). Compressed nonnegative
  91. % matrix factorization is fast and accurate. IEEE Transactions
  92. % on Signal Processing, 64(9), 2269-2283.
  93. % see: https://arxiv.org/pdf/1505.04650
  94. % The corresponding code is originally created by the authors
  95. % Then, it is modified by F. Yahaya.
  96. % Date: 13/04/2018
  97. %
  98. compressionLevel=2;
  99. [m,n]=size(X);
  100. l = min(min(n,m), max(compressionLevel, r ));
  101. switch nargin
  102. case 2
  103. OmegaL = randn(n,l);
  104. OmegaR = randn(l, m);
  105. q = 4;
  106. case 4
  107. OmegaL = varargin{2};
  108. OmegaR = varargin{1};
  109. q = 1;
  110. end
  111. Y = X * OmegaL;
  112. for i=1:q
  113. [Y,~]=qr(Y,0);
  114. S=X'*Y;
  115. [Z,~]=qr(S,0);
  116. Y=X* Z;
  117. end
  118. [L,~]=qr(Y,0);
  119. L=L';
  120. Y = OmegaR * X;
  121. for i=1:q
  122. [Y,~]=qr(Y',0);
  123. S=X*Y;
  124. [Z,~]=qr(S,0);
  125. Y=Z'*X;
  126. end
  127. Y=Y';
  128. [R,~] = qr(Y,0);
  129. end