Add Hesian computation for objective minimization

best/cmem/solver/be_cmem_pipelineoptions.m

@@ -3,6 +3,13 @@
         if nargin < 1 || isempty(DataTypes)
             DataTypes = {'EEG'};
         end
+
+        % Datatypes
+        if ischar(DataTypes)
+            DEF.mandatory.DataTypes = {DataTypes};
+        else
+            DEF.mandatory.DataTypes = DataTypes;
+        end
 
         % clustering
         DEF.clustering.clusters_type     	= 'static';

best/main/be_main.m

@@ -128,7 +128,7 @@
 
 % ----------------------------------------------------------------------- %
 
-function Def_OPTIONS=   BEst_defaults()
+function Def_OPTIONS =   BEst_defaults()
 
     % ===== common I/O arguments ===== %%
 
@@ -218,6 +218,7 @@
     Def_OPTIONS.solver.Optim_method                 = 'fminunc';
     Def_OPTIONS.solver.covariance_scale             = 1;
     Def_OPTIONS.solver.parallel_matlab              = be_canUseParallelPool();
+    Def_OPTIONS.solver.useHessian                   = 0;
 
     % output options
     Def_OPTIONS.output.save_factor                  =  1;

best/main/be_main_data_preprocessing.m

@@ -91,8 +91,10 @@
             % OPTIONS.mandatory.Data                            = OPTIONS.automatic.Modality(ii).emptyroom;
             OPTIONS.automatic.Modality(ii).data             	= OPTIONS.automatic.Modality(ii).emptyroom;
             OPTIONS.mandatory.DataTime                       	= OPTIONS.automatic.Emptyroom_time;  
-            fprintf('MEM : New noise covariance is computed using empty room data\n');
-
+
+            if OPTIONS.optional.verbose
+                fprintf('MEM : New noise covariance is computed using empty room data\n');
+            end
         elseif ~isempty( OPTIONS.automatic.Modality(ii).baseline )
             idNb = any(isnan(OPTIONS.automatic.Modality(ii).baseline));
             OPTIONS.automatic.Modality(ii).baseline(:,idNb)     = 0;
@@ -103,10 +105,14 @@
                 TLM                                             =   be_closest(OPTIONS.optional.BaselineSegment([1 end]), OPTIONS.mandatory.DataTime);
                 OPTIONS.mandatory.DataTime                      = OPTIONS.mandatory.DataTime( TLM(1):TLM(end) );      
             end
-            fprintf('MEM : New noise covariance is computed using baseline\n');
-
+
+            if OPTIONS.optional.verbose
+                fprintf('MEM : New noise covariance is computed using baseline\n');
+            end
         else
-        	fprintf('MEM : No baseline nor covariance matrix provided. Covariance set to identity\n');
+            if OPTIONS.optional.verbose
+        	    fprintf('MEM : No baseline nor covariance matrix provided. Covariance set to identity\n');
+            end
             OPTIONS.automatic.Modality(ii).baseline             =   eye(OPTIONS.automatic.Modality(ii).channels);
             OPTIONS.solver.NoiseCov_method                      =   0;
         end

best/mem/be_free_energy_simplified.m

@@ -0,0 +1,88 @@
+function [D, dD, H] = be_free_energy_simplified(lambda, M, noise_var, G_active_var_Gt, clusters)
+%CALCULATE_FREE_ENERGY gives the free energy of a system
+%   [D, dD] = CALCULATE_FREE_ENERGY(LAMBDA, M, NOISE_VAR, CLUSTERS, NB_CLUSTERS)
+%   calculates the free energy of the system for a given LAMBDA and 
+%   returns it in D. The function also returns the derivative of D with 
+%   respect to LAMBDA in dD.  This method should not be called by itself.
+%   
+%   Warning: This function assume that the active mean is 0, and inactive
+%   parcels are dirac at 0. For more options, use be_free_energy. 
+%   This limitation allows us to also compute the Hessian H. 
+% 
+%   The method uses the methodoly described in :
+%   Amblard, C., Lapalme, E., Lina, J. 2004. Biomagnetic Source Detection
+%       by Maximyum Entropy and Graphical  Models, IEEE Transactions on 
+%       Biomedical Engineering, vol. 51, no 3, p. 427-442.
+%
+%   The formulas are :
+%       D(lambda) = lambda' * M - 
+%                   (1/2)* noise_var * lambda' * lambda - 
+%                   sum(F*(G' * lambda))
+%
+%       F*(xi) = ln[(1- alpha) * exp(F0) + alpha * exp(F1)]
+%           where F0 is the inactive state and
+%                 F1 is the active state
+%       F0(xi) = 1/2 * xi' * omega * xi
+%       F1(xi) = 1/2 * xi' * sigma * xi + xi' * mu
+%
+%% ==============================================
+% Copyright (C) 2011 - LATIS Team
+%
+%  Authors: LATIS team, 2011
+%
+%% ==============================================
+% License 
+%
+% BEst is free software: you can redistribute it and/or modify
+%    it under the terms of the GNU General Public License as published by
+%    the Free Software Foundation, either version 3 of the License, or
+%    (at your option) any later version.
+%
+%    BEst is distributed in the hope that it will be useful,
+%    but WITHOUT ANY WARRANTY; without even the implied warranty of
+%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+%    GNU General Public License for more details.
+%
+%    You should have received a copy of the GNU General Public License
+%    along with BEst. If not, see <http://www.gnu.org/licenses/>.
+% -------------------------------------------------------------------------
+
+    lambda_trans = lambda';      
+    p            = [clusters.active_probability];
+
+    % Estimate dF1 and F1 (separating the contribution  of the mean and
+    % covariance for optimization purpose)
+    dF1     = squeeze(be_pagemtimes(G_active_var_Gt,lambda));
+    F1      =  1/2 * lambda_trans*dF1; 
+
+    % Finally estimate dF and F
+    coeffs_free_energy  = (1-p) .* exp(-F1)  +  p;
+    s1   = p ./ coeffs_free_energy;
+
+    F = F1 + log(coeffs_free_energy);
+    dF  = s1 .* dF1;
+
+
+    % The outcome of the equations produces a strictly convex function (with a maximum).
+    D   = -(lambda_trans * M - sum(F) - (1/2) * lambda_trans * noise_var * lambda);
+
+    % Compute the gradient
+    if nargout >= 2
+        dD  = -(               M - sum(dF,2)                     - noise_var * lambda);
+    end
+
+    % Compute the Hessian
+    if nargout >= 3
+
+        % Curvature term: sum_i s1_i * G*Sigma_i*G'               
+        curvature_term = squeeze(pagemtimes(permute(G_active_var_Gt, [2, 3, 1]), s1'));
+
+        % Rank-1 term: sum_i s1_i*(1-s1_i) * delta_i*delta_i'
+        w_delta      = dF1 .* sqrt(s1 .* (1 - s1));
+        rank1_term   = w_delta * w_delta'; 
+
+        H = curvature_term + rank1_term + noise_var;
+    end
+
+
+end

best/mem/be_memstruct.m

@@ -80,6 +80,8 @@
 % some basic parameters
 nb_clusters = max(obj.clusters);
 nb_sources  = obj.nb_sources;
+nb_sensors   = length(obj.data);
+
 
 % some definitions 
 clusters_struct(nb_clusters) = struct;
@@ -88,6 +90,7 @@
 active_mean = cell(1,nb_clusters);
 inactive_var= cell(1,nb_clusters);
 active_var_out = zeros(obj.nb_sources,1);
+G_active_var_Gt = zeros(nb_sensors, nb_sensors, nb_clusters);
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % the following loop goes though each of the clusters (non null clusters)
@@ -109,7 +112,7 @@
 
     % Multiply the active variance by the MNE energy
     obj.active_var{ii} = obj.active_var{ii}     * obj.mne_energy(ii);
-    obj.G_active_var_Gt{ii}     = obj.G_active_var_Gt{ii} * obj.mne_energy(ii);
+    G_active_var_Gt(:, :, ii)   = obj.G_active_var_Gt{ii} * obj.mne_energy(ii);
     active_var_out(idx_cluster) = diag( obj.active_var{ii} );
 
     % SIGMA0: variance of the inactive state (not relevant for the present version)
@@ -128,7 +131,6 @@
 [clusters_struct.active_mean]        = active_mean{:};
 [clusters_struct.active_var]         = obj.active_var{:};
 [clusters_struct.inactive_var]       = inactive_var{:};  
-[clusters_struct.G_active_var_Gt]    = obj.G_active_var_Gt{:};  
 [clusters_struct.indices]            = cID{:};
 
 % Creation of the structure to solve the MEM
@@ -139,16 +141,17 @@
 mem.nb_clusters                      = nb_clusters;
 mem.optim_algo                       = OPTIONS.solver.Optim_method;
 mem.optimoptions                     = OPTIONS.solver.optimoptions;
-mem.G_active_var_Gt                  = cat(3,obj.G_active_var_Gt{:});
+mem.useHessian                       = OPTIONS.solver.useHessian;
+mem.G_active_var_Gt                  = G_active_var_Gt;
 
 % Lambda initialization
 switch OPTIONS.model.initial_lambda 
 
     case 0
-        mem.lambda	=   zeros(size(mem.M));
+        mem.lambda	=   zeros(nb_sensors, 1);
 
     case 1
-        mem.lambda	=   randn(size(mem.M)) / mean( abs(mem.M) ); % initial value for threshold diff 0
+        mem.lambda	=   randn(nb_sensors, 1) / mean( abs(mem.M), 1 ); 
 
 end
 

best/mem/be_minimize_free_energy.m

@@ -35,21 +35,31 @@
 %    along with BEst. If not, see <http://www.gnu.org/licenses/>.
 % -------------------------------------------------------------------------
 
-    if strcmpi(mem_struct.optim_algo, 'fminunc')  && ...
-       license('test', 'Optimization_Toolbox') && ...
-       exist('fminunc', 'file')
+    if strcmpi(mem_struct.optim_algo, 'fminunc')
 
         % call the unconstrained minimization routine from MATLAB optimization
         %  toolbox, if available.
-
-        [lambda, entropy, ~, out]  = fminunc(...
-            @be_free_energy , ...
-            mem_struct.lambda, ...
-            mem_struct.optimoptions, ...
-            mem_struct.M, ...
-            mem_struct.noise_var, ...
-            mem_struct.G_active_var_Gt, ...
-            mem_struct.clusters);    
+        if mem_struct.useHessian
+
+            [lambda, entropy, ~, out]  = fminunc(...
+                @be_free_energy_simplified , ...
+                mem_struct.lambda, ...
+                mem_struct.optimoptions, ...
+                mem_struct.M, ...
+                mem_struct.noise_var, ...
+                mem_struct.G_active_var_Gt, ...
+                mem_struct.clusters);  
+        else
+
+            [lambda, entropy, ~, out]  = fminunc(...
+                @be_free_energy , ...
+                mem_struct.lambda, ...
+                mem_struct.optimoptions, ...
+                mem_struct.M, ...
+                mem_struct.noise_var, ...
+                mem_struct.G_active_var_Gt, ...
+                mem_struct.clusters);    
+        end
     else
 
         % Call an alternate unconstrained minimization routine 

best/mem/be_slice_obj.m

@@ -2,10 +2,11 @@
 
     nbSmp       = size(Data,2); 
     nb_sensors  = size(Data,1); 
-
     obj_slice(nbSmp)    = struct();
 
-    fprintf('%s, finalizing MEM prior ...', OPTIONS.mandatory.pipeline);
+    if OPTIONS.optional.verbose
+        fprintf('%s, finalizing MEM prior ...', OPTIONS.mandatory.pipeline);
+    end
 
     for i = 1:nbSmp
 
@@ -113,38 +114,15 @@
 
 
     % Estimate the active variance 
-    % Multiply Signa_s by 5% of the MNE solution
-    if strcmp(OPTIONS.clustering.clusters_type, 'static')
-        clusters = obj.CLS(:,1);
-
-        for i = 1:nbSmp
-
-            energy = zeros(max(clusters), 1);
-
-            Jmne = OPTIONS.automatic.Modality(1).MneKernel * obj_slice(i).data;
-            Jmne = Jmne  ./ max(abs(Jmne));
-
-            for ii = 1:max(clusters)
-                energy(ii)  = OPTIONS.solver.active_var_mult * mean(Jmne(clusters == ii).^2);
-            end
-
-            obj_slice(i).mne_energy = energy;            
-        end
-    else
-        for i = 1:nbSmp
-            clusters = obj.CLS(:,i);
-            energy = zeros(max(clusters),1);
+    for i = 1:nbSmp
+        clusters = obj.CLS(:,i);
 
-            Jmne = OPTIONS.automatic.Modality(1).MneKernel * obj_slice(i).data;
-            Jmne = Jmne  ./ max(abs(Jmne));
+        Jmne = OPTIONS.automatic.Modality(1).MneKernel * obj_slice(i).data;
+        Jmne = Jmne  ./ max(abs(Jmne));
 
-            for ii = 1:max(clusters)
-                energy(ii) = OPTIONS.solver.active_var_mult * mean(Jmne(clusters == ii).^2);        
-            end
-            obj_slice(i).mne_energy = energy;
-        end
+        energy      =  accumarray(clusters, Jmne .^ 2, [max(clusters), 1], @(x)mean(x,1));
+        obj_slice(i).mne_energy = OPTIONS.solver.active_var_mult *  energy; 
     end
-
 
     obj_const.gain          = obj.gain;
     obj_const.nb_sources    = obj.nb_sources;
@@ -160,12 +138,26 @@
     if strcmpi(OPTIONS.solver.Optim_method, 'fminunc') && ...
        license('test', 'Optimization_Toolbox') && ...
        exist('fminunc', 'file')
+
+        if OPTIONS.solver.useHessian
+            OPTIONS.solver.optimoptions =   optimoptions('fminunc',...
+                                                    'MaxIter', 1000, ...
+                                                    'MaxFunEvals', 1000, ...
+                                                    'algorithm', 'trust-region', ...
+                                                    'SpecifyObjectiveGradient',true, ...
+                                                    'GradObj', 'on', ...
+                                                    'HessianFcn', 'objective', ...
+                                                    'Display', 'off');  
+
+
+        else
 
-        OPTIONS.solver.optimoptions =   optimoptions('fminunc','GradObj', 'on', ...
-                                                    'MaxIter', MAX_ITER, ...
-                                                    'MaxFunEvals', MAX_ITER, ...
-                                                    'algorithm', 'quasi-newton',... % 'quasi-newton' trust-region'
-                                                    'Display', 'off' );
+            OPTIONS.solver.optimoptions =   optimoptions('fminunc','GradObj', 'on', ...
+                                                        'MaxIter', MAX_ITER, ...
+                                                        'MaxFunEvals', MAX_ITER, ...
+                                                        'algorithm', 'quasi-newton',... % 'quasi-newton' trust-region'
+                                                        'Display', 'off' );
+        end
     else
 
         if strcmpi(OPTIONS.solver.Optim_method, 'fminunc')
@@ -187,6 +179,7 @@
         OPTIONS.solver.optimoptions = options;
     end
 
-    fprintf('done. \n');
-
+    if OPTIONS.optional.verbose
+        fprintf('done. \n');    
+    end
 end

best/mne/solver/be_cmne_pipelineoptions.m

@@ -4,6 +4,13 @@
             DataTypes = {'EEG'};
         end
 
+        % Datatypes
+        if ischar(DataTypes)
+            DEF.mandatory.DataTypes = {DataTypes};
+        else
+            DEF.mandatory.DataTypes = DataTypes;
+        end
+
         % clustering
         DEF.clustering.clusters_type     	= 'static';
         DEF.solver.mne_method               = 'mne_lcurve';

best/mne/solver/be_jmne_lcurve.m

@@ -42,9 +42,11 @@
     end
 
     G = obj.gain; 
-
-    fprintf('%s, solving MNE by L-curve ...', OPTIONS.mandatory.pipeline);
 
+    if OPTIONS.optional.verbose
+        fprintf('%s, solving MNE by L-curve ...', OPTIONS.mandatory.pipeline);
+    end
+
     p       = OPTIONS.model.depth_weigth_MNE;
     % Compute covariance matrices
     if OPTIONS.solver.mne_use_noiseCov && ~isempty(OPTIONS.automatic.Modality(1).covariance) && size(OPTIONS.automatic.Modality(1).covariance,3) == 1
@@ -132,8 +134,11 @@
     if ~OPTIONS.automatic.stand_alone
         bst_progress('stop');
     end
-    fprintf('done. \n');
 
+    if OPTIONS.optional.verbose
+        fprintf('done. \n');
+    end
+
     if OPTIONS.optional.display
         if isempty(sfig.hfig)
             sfig.hfig =  figure();