Weighted Minimum Norm solution for the static inverse problem.
Solve the following minimization problem
min {E(J) = (B-GJ)'X(B-GJ) + lam*J'WJ}
with optimization of 'lam' by ABIC or GCV.
When 'option.lapacian' exists, the following minimization problem
min {E(J) = (B-GJ)'X(B-GJ) + lam*J'LWLJ}
is solved ('LORETA' solution).
[Syntax]
> [lambda, sigma2, ABIC, Jest, F, Info] = wmn_solution(B, G, Ceps, W, option)
[Example]
>> option.criteria = 'ABIC'
>> option.flag_current = 1;
>> [lambda, sigma2, ABIC, Jest, F] = wmn_solution(B, G, Ceps, W, option)
returns optimal parameters, estimates of J, and the inverse filter using ABIC
as the statisitcal criteria for optimization.
[Output]
optlambda : optimum regularization parameter
optsigma2 : optimum observation variance
optcriteria : the value of criteria (ABIC or GCV) at optlambda
Jest : estimated J
F : inverse filter (mapping from B to J)
Info : information of the optimization
1st column : lambda
2nd column : criteria (ABIC or GCV)
3rd column : observation noise variance (SIG2)
[Input]
B : EEG/MEG observation of a single trial (Nch*Nt matrix)
G : lead field matrix (Nch * Nj matrix)
X : observation weight matrix (inverse of observation noise)
w : the diagonal elements of prior weight matrix W (should be diagonal, because of memory problem)
[Optional Input]
option.criteria : criteria for the estimation of reguralization parameter
option.laplacian : if exist, W is replaced by LWL (L is lapalcian)
option.flag_plot : if 1, the values of 'criteria' (ABIC or GCV) are displayed with respect to 'lambda'
option.flag_current : if 1, the current is reconstrcted in the third output (default:1)
option.Npoints : the number of grid points in lambda axis to be explored (default:100)
option.lambda : fix the regularization parameter (default : [])
option.search_margin : parameter to define range of lambda to be searched (default:5).
2016-03-22 O.Yamashita
* allow empty X and w to speed up (avoid time-consuming sqrtm function)
2014-03-03 O.Yamashita
* add a new field 'search_margin' to 'option' struct
2014-02-10 O.Yamashita
* add a new field 'lambda' to 'option' struct
* modify comments
2012-11-07 O.Yamashita
* save information for the optimization in 'Info'
2012-04-26 O.Yamashita
* bug fix for the lines to compute square root of a matrix (sqrt --> sqrtm)
previous version causes complex values
* add inverse filter 'F' as an output argument (important)
2005-04-26 O.Yamashita beta-version,0001 function [optlambda, optsigma2, optcriteria, Jest, F, Info] = wmn_solution(B, G, X, w, option) 0002 % Weighted Minimum Norm solution for the static inverse problem. 0003 % 0004 % Solve the following minimization problem 0005 % min {E(J) = (B-GJ)'X(B-GJ) + lam*J'WJ} 0006 % with optimization of 'lam' by ABIC or GCV. 0007 % When 'option.lapacian' exists, the following minimization problem 0008 % min {E(J) = (B-GJ)'X(B-GJ) + lam*J'LWLJ} 0009 % is solved ('LORETA' solution). 0010 % 0011 % [Syntax] 0012 % > [lambda, sigma2, ABIC, Jest, F, Info] = wmn_solution(B, G, Ceps, W, option) 0013 % 0014 % [Example] 0015 % >> option.criteria = 'ABIC' 0016 % >> option.flag_current = 1; 0017 % >> [lambda, sigma2, ABIC, Jest, F] = wmn_solution(B, G, Ceps, W, option) 0018 % 0019 % returns optimal parameters, estimates of J, and the inverse filter using ABIC 0020 % as the statisitcal criteria for optimization. 0021 % 0022 % [Output] 0023 % optlambda : optimum regularization parameter 0024 % optsigma2 : optimum observation variance 0025 % optcriteria : the value of criteria (ABIC or GCV) at optlambda 0026 % Jest : estimated J 0027 % F : inverse filter (mapping from B to J) 0028 % Info : information of the optimization 0029 % 1st column : lambda 0030 % 2nd column : criteria (ABIC or GCV) 0031 % 3rd column : observation noise variance (SIG2) 0032 % 0033 % [Input] 0034 % B : EEG/MEG observation of a single trial (Nch*Nt matrix) 0035 % G : lead field matrix (Nch * Nj matrix) 0036 % X : observation weight matrix (inverse of observation noise) 0037 % w : the diagonal elements of prior weight matrix W (should be diagonal, because of memory problem) 0038 % 0039 % [Optional Input] 0040 % option.criteria : criteria for the estimation of reguralization parameter 0041 % option.laplacian : if exist, W is replaced by LWL (L is lapalcian) 0042 % option.flag_plot : if 1, the values of 'criteria' (ABIC or GCV) are displayed with respect to 'lambda' 0043 % option.flag_current : if 1, the current is reconstrcted in the third output (default:1) 0044 % option.Npoints : the number of grid points in lambda axis to be explored (default:100) 0045 % option.lambda : fix the regularization parameter (default : []) 0046 % option.search_margin : parameter to define range of lambda to be searched (default:5). 0047 % 0048 % 2016-03-22 O.Yamashita 0049 % * allow empty X and w to speed up (avoid time-consuming sqrtm function) 0050 % 2014-03-03 O.Yamashita 0051 % * add a new field 'search_margin' to 'option' struct 0052 % 2014-02-10 O.Yamashita 0053 % * add a new field 'lambda' to 'option' struct 0054 % * modify comments 0055 % 2012-11-07 O.Yamashita 0056 % * save information for the optimization in 'Info' 0057 % 2012-04-26 O.Yamashita 0058 % * bug fix for the lines to compute square root of a matrix (sqrt --> sqrtm) 0059 % previous version causes complex values 0060 % * add inverse filter 'F' as an output argument (important) 0061 % 2005-04-26 O.Yamashita beta-version, 0062 0063 %profile on 0064 0065 [Nch, Nt] = size(B); 0066 [Nj] = size(G,2); 0067 0068 0069 % default 0070 criteria = 'ABIC'; 0071 flag_plot = 1; 0072 flag_current = 1; 0073 L = speye(Nj); 0074 Npoints = 100; % Number of points on the curve. 0075 lambda = []; 0076 search_margin = 5; % in log scale 0077 0078 %% Error check 0079 if exist('option') == 1 & isfield(option, 'criteria') 0080 criteria = option.criteria; 0081 end 0082 if exist('option') == 1 & isfield(option, 'flag_plot') 0083 flag_plot = option.flag_plot ; 0084 end 0085 if exist('option') == 1 & isfield(option, 'flag_current') 0086 flag_current = option.flag_current; 0087 end 0088 if exist('option') == 1 & isfield(option, 'laplacian') 0089 L = option.laplacian; 0090 L = L / max(L(:)); 0091 end 0092 if exist('option') == 1 & isfield(option, 'Npoints') 0093 Npoints = option.Npoints; 0094 end 0095 if exist('option') == 1 & isfield(option, 'lambda') 0096 lambda = option.lambda; 0097 end 0098 if exist('option') == 1 & isfield(option, 'search_margin') 0099 search_margin = option.search_margin; 0100 end 0101 0102 % pre-whitening 0103 if isempty(X) 0104 Xsqrt = eye(Nch); 0105 else 0106 Xsqrt = sqrtm(X); 0107 end 0108 if isempty(w) 0109 wsqrt = ones(Nj,1); 0110 else 0111 wsqrt = sqrt(w); 0112 end 0113 0114 Btil = Xsqrt * B; 0115 Gtil = Xsqrt * G / L * spdiags(1./wsqrt, 0, Nj, Nj) ; 0116 0117 [U s V] = csvd1(Gtil); 0118 0119 % Sets grid points 0120 if isempty(lambda) 0121 minloglx = log(s(end)) - search_margin; 0122 maxloglx = log(s(1)) + search_margin; 0123 bin = (maxloglx-minloglx)/(Npoints-1); 0124 loglx_points = [minloglx : bin : maxloglx]; 0125 else 0126 loglx_points = log(lambda); % fixed lambda 0127 end 0128 0129 Npoints = length(loglx_points); 0130 SIG2 = zeros(Npoints,1); 0131 CRITERIA = zeros(Npoints,1); 0132 0133 UtB = U'*Btil; 0134 BB = UtB * UtB'; 0135 s2=s.*s; 0136 0137 % find optimum regularization parameter 0138 switch criteria 0139 case 'ABIC' 0140 %% calculate ABIC 0141 for i = 1 : Npoints 0142 lam = exp(loglx_points(i)); 0143 l2=lam; 0144 r = l2./(s2+l2); 0145 sig2 = sum(sum(diag(r).*BB))/Nch/Nt; %tr(diag(r)*BB) 0146 ABIC1=Nch*Nt*log(sig2); 0147 ABIC2= -Nt*sum(log(r)); 0148 CRITERIA(i)=ABIC1+ABIC2; 0149 SIG2(i) = sig2; 0150 end 0151 0152 case 'GCV' 0153 %% calculate GCV 0154 for i = 1 : Npoints 0155 lam = exp(loglx_points(i)); 0156 l2=lam; 0157 r = l2./(s2+l2); 0158 nume = sum(sum(diag(r.^2).*BB)); 0159 denom = (sum(r))^2; 0160 GCV = nume/denom; 0161 sig2 = sum(sum(diag(r).*BB))/Nch/Nt; %tr(diag(r)*BB) 0162 CRITERIA(i) = GCV; 0163 SIG2(i) = sig2; 0164 end 0165 0166 otherwise 0167 disp('No such a method'); 0168 end 0169 0170 [tmp, ind] = min(CRITERIA); 0171 optlambda = exp(loglx_points(ind)); 0172 optsigma2 = SIG2(ind); 0173 optcriteria = CRITERIA(ind); 0174 Info = [exp(loglx_points(:)) CRITERIA(:) SIG2(:)]; 0175 0176 %% current reconstruction 0177 0178 if flag_current 0179 l2= optlambda; 0180 r = s./(s2+l2); 0181 Jtil = V *diag(r)*UtB; 0182 Jest = L\ Jtil ./ wsqrt(:,ones(1,Nt)); %Jest = L\ spdiags(1./wsqrt, 0, Nj, Nj) * Jtil; 0183 F = L\spdiags(1./wsqrt, 0, Nj, Nj)*V*diag(r)*U'*Xsqrt; 0184 else 0185 Jest = []; 0186 F = []; 0187 end 0188 0189 if flag_plot 0190 figure 0191 [AX,H1,H2] = plotyy(loglx_points, CRITERIA, loglx_points, SIG2); 0192 xlabel('Log of lambda'); 0193 set(get(AX(1),'Ylabel'),'String','ABIC or GCV') 0194 set(get(AX(2),'Ylabel'),'String','Observation Noise Variance'); 0195 hold on; 0196 plot([loglx_points(ind) loglx_points(ind)], [min(get(AX(1),'ytick')) max(get(AX(1),'ytick'))], 'k--'); 0197 end 0198 0199 0200 %profile viewer 0201 0202