From 3895221063a4701acda6a73539aefdee11659823 Mon Sep 17 00:00:00 2001 From: Edouard Date: Mon, 4 May 2026 16:41:25 -0400 Subject: [PATCH 1/7] Add Hesian computation for objective minimization --- best/main/be_main.m | 1 + best/mem/be_free_energy_simplified.m | 89 ++++++++++++++++++++++++++++ best/mem/be_memstruct.m | 1 + best/mem/be_minimize_free_energy.m | 30 +++++++--- best/mem/be_slice_obj.m | 25 ++++++-- 5 files changed, 132 insertions(+), 14 deletions(-) create mode 100644 best/mem/be_free_energy_simplified.m diff --git a/best/main/be_main.m b/best/main/be_main.m index 323b543..e01edc9 100644 --- a/best/main/be_main.m +++ b/best/main/be_main.m @@ -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; diff --git a/best/mem/be_free_energy_simplified.m b/best/mem/be_free_energy_simplified.m new file mode 100644 index 0000000..5be162f --- /dev/null +++ b/best/mem/be_free_energy_simplified.m @@ -0,0 +1,89 @@ +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 . +% ------------------------------------------------------------------------- + + 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' + nSensor = size(M,1); + curvature_term = reshape( reshape(G_active_var_Gt, nSensor * nSensor, []) * s1', nSensor, nSensor); + + % 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 diff --git a/best/mem/be_memstruct.m b/best/mem/be_memstruct.m index e5bc6d5..67f0d01 100644 --- a/best/mem/be_memstruct.m +++ b/best/mem/be_memstruct.m @@ -139,6 +139,7 @@ mem.nb_clusters = nb_clusters; mem.optim_algo = OPTIONS.solver.Optim_method; mem.optimoptions = OPTIONS.solver.optimoptions; +mem.useHessian = OPTIONS.solver.useHessian; mem.G_active_var_Gt = cat(3,obj.G_active_var_Gt{:}); % Lambda initialization diff --git a/best/mem/be_minimize_free_energy.m b/best/mem/be_minimize_free_energy.m index e07a2ad..a784b54 100644 --- a/best/mem/be_minimize_free_energy.m +++ b/best/mem/be_minimize_free_energy.m @@ -41,15 +41,27 @@ % 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 diff --git a/best/mem/be_slice_obj.m b/best/mem/be_slice_obj.m index c90c776..9a28ac6 100644 --- a/best/mem/be_slice_obj.m +++ b/best/mem/be_slice_obj.m @@ -160,12 +160,27 @@ if strcmpi(OPTIONS.solver.Optim_method, 'fminunc') && ... license('test', 'Optimization_Toolbox') && ... exist('fminunc', 'file') + + if OPTIONS.solver.useHessian + disp('Using Hessian') + OPTIONS.solver.optimoptions = optimoptions('fminunc',... + 'MaxIter', 1000, ... + 'MaxFunEvals', 1000, ... + 'algorithm', 'trust-region', ... + 'SpecifyObjectiveGradient',true, ... + 'GradObj', 'on', ... + 'HessianFcn', 'objective', ... + 'Display', 'off'); + - OPTIONS.solver.optimoptions = optimoptions('fminunc','GradObj', 'on', ... - 'MaxIter', MAX_ITER, ... - 'MaxFunEvals', MAX_ITER, ... - 'algorithm', 'quasi-newton',... % 'quasi-newton' trust-region' - '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' ); + end else if strcmpi(OPTIONS.solver.Optim_method, 'fminunc') From 873b7eb9b773638f37ae7b5a855b5db9bc0e8801 Mon Sep 17 00:00:00 2001 From: Edouard Date: Tue, 5 May 2026 21:09:15 -0400 Subject: [PATCH 2/7] Correctly fill the datatype default option based on the modality --- best/cmem/solver/be_cmem_pipelineoptions.m | 7 +++++++ best/main/be_main.m | 2 +- best/mne/solver/be_cmne_pipelineoptions.m | 7 +++++++ best/rmem/solver/be_rmem_pipelineoptions.m | 7 +++++++ best/wmem/solver/be_wmem_pipelineoptions.m | 7 +++++++ 5 files changed, 29 insertions(+), 1 deletion(-) diff --git a/best/cmem/solver/be_cmem_pipelineoptions.m b/best/cmem/solver/be_cmem_pipelineoptions.m index be0da7e..272578b 100644 --- a/best/cmem/solver/be_cmem_pipelineoptions.m +++ b/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'; diff --git a/best/main/be_main.m b/best/main/be_main.m index e01edc9..c12c73d 100644 --- a/best/main/be_main.m +++ b/best/main/be_main.m @@ -128,7 +128,7 @@ % ----------------------------------------------------------------------- % -function Def_OPTIONS= BEst_defaults() +function Def_OPTIONS = BEst_defaults() % ===== common I/O arguments ===== %% diff --git a/best/mne/solver/be_cmne_pipelineoptions.m b/best/mne/solver/be_cmne_pipelineoptions.m index 00028c5..7baee39 100644 --- a/best/mne/solver/be_cmne_pipelineoptions.m +++ b/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'; diff --git a/best/rmem/solver/be_rmem_pipelineoptions.m b/best/rmem/solver/be_rmem_pipelineoptions.m index e9558cc..56f5fce 100644 --- a/best/rmem/solver/be_rmem_pipelineoptions.m +++ b/best/rmem/solver/be_rmem_pipelineoptions.m @@ -2,6 +2,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 = 'blockwise'; diff --git a/best/wmem/solver/be_wmem_pipelineoptions.m b/best/wmem/solver/be_wmem_pipelineoptions.m index 1b2299f..378b39b 100644 --- a/best/wmem/solver/be_wmem_pipelineoptions.m +++ b/best/wmem/solver/be_wmem_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'; From a9c4856aa10da7fe32050d4226ff3f0575596f4e Mon Sep 17 00:00:00 2001 From: Edouard Date: Wed, 6 May 2026 11:43:22 -0400 Subject: [PATCH 3/7] Clarify code intent --- best/mem/be_free_energy_simplified.m | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/best/mem/be_free_energy_simplified.m b/best/mem/be_free_energy_simplified.m index 5be162f..e392c27 100644 --- a/best/mem/be_free_energy_simplified.m +++ b/best/mem/be_free_energy_simplified.m @@ -74,10 +74,9 @@ % Compute the Hessian if nargout >= 3 - % Curvature term: sum_i s1_i * G*Sigma_i*G' - nSensor = size(M,1); - curvature_term = reshape( reshape(G_active_var_Gt, nSensor * nSensor, []) * s1', nSensor, nSensor); - + % 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'; From b3f41d9432ae0a74fc49a54630442d894c19714e Mon Sep 17 00:00:00 2001 From: Edouard Date: Wed, 6 May 2026 13:29:24 -0400 Subject: [PATCH 4/7] small optim --- best/mem/be_minimize_free_energy.m | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/best/mem/be_minimize_free_energy.m b/best/mem/be_minimize_free_energy.m index a784b54..13b4e51 100644 --- a/best/mem/be_minimize_free_energy.m +++ b/best/mem/be_minimize_free_energy.m @@ -35,9 +35,7 @@ % along with BEst. If not, see . % ------------------------------------------------------------------------- - 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. From aea0b95acea6ae2b79f724e26986475f69bee9d0 Mon Sep 17 00:00:00 2001 From: Edouard Date: Wed, 6 May 2026 14:08:56 -0400 Subject: [PATCH 5/7] update the estimation of G_active_var_Gt - avoid tmp allocation --- best/mem/be_memstruct.m | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/best/mem/be_memstruct.m b/best/mem/be_memstruct.m index 67f0d01..734bbad 100644 --- a/best/mem/be_memstruct.m +++ b/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 @@ -140,7 +142,7 @@ mem.optim_algo = OPTIONS.solver.Optim_method; mem.optimoptions = OPTIONS.solver.optimoptions; mem.useHessian = OPTIONS.solver.useHessian; -mem.G_active_var_Gt = cat(3,obj.G_active_var_Gt{:}); +mem.G_active_var_Gt = G_active_var_Gt; % Lambda initialization switch OPTIONS.model.initial_lambda From 73ea5858ad59466b0b6f9559d413776b31af1093 Mon Sep 17 00:00:00 2001 From: Edouard Date: Wed, 6 May 2026 14:46:13 -0400 Subject: [PATCH 6/7] minor perf improvement --- best/mem/be_memstruct.m | 4 ++-- best/mem/be_slice_obj.m | 35 ++++++----------------------------- 2 files changed, 8 insertions(+), 31 deletions(-) diff --git a/best/mem/be_memstruct.m b/best/mem/be_memstruct.m index 734bbad..0e99262 100644 --- a/best/mem/be_memstruct.m +++ b/best/mem/be_memstruct.m @@ -148,10 +148,10 @@ 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 diff --git a/best/mem/be_slice_obj.m b/best/mem/be_slice_obj.m index 9a28ac6..c688a68 100644 --- a/best/mem/be_slice_obj.m +++ b/best/mem/be_slice_obj.m @@ -113,38 +113,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; From 7bf696653ee9d8f1cd856cd532fdc695788005f1 Mon Sep 17 00:00:00 2001 From: Edouard Date: Wed, 6 May 2026 15:30:50 -0400 Subject: [PATCH 7/7] avoid display when verbose is false --- best/main/be_main_data_preprocessing.m | 16 +++++++++++----- best/mem/be_slice_obj.m | 11 ++++++----- best/mne/solver/be_jmne_lcurve.m | 11 ++++++++--- 3 files changed, 25 insertions(+), 13 deletions(-) diff --git a/best/main/be_main_data_preprocessing.m b/best/main/be_main_data_preprocessing.m index dc35888..63c4aa6 100644 --- a/best/main/be_main_data_preprocessing.m +++ b/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 diff --git a/best/mem/be_slice_obj.m b/best/mem/be_slice_obj.m index c688a68..25852fc 100644 --- a/best/mem/be_slice_obj.m +++ b/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 @@ -139,7 +140,6 @@ exist('fminunc', 'file') if OPTIONS.solver.useHessian - disp('Using Hessian') OPTIONS.solver.optimoptions = optimoptions('fminunc',... 'MaxIter', 1000, ... 'MaxFunEvals', 1000, ... @@ -179,6 +179,7 @@ OPTIONS.solver.optimoptions = options; end - fprintf('done. \n'); - + if OPTIONS.optional.verbose + fprintf('done. \n'); + end end diff --git a/best/mne/solver/be_jmne_lcurve.m b/best/mne/solver/be_jmne_lcurve.m index 008b8cf..4a8ecb4 100644 --- a/best/mne/solver/be_jmne_lcurve.m +++ b/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();