L2ノルム

センサー数をM個、格子点数をN個とし、格子点1つにつき電流双極子の自由度が２つあるとします。 そのときの以下の式が成り立つとします。

行列で書きます。

2N>Mであれば、この式を満たすqは幾らでも存在します。そこで以下のような制約を設けます。

pはpノルムの意味です。L1およびL2ノルムの場合は

となります。L2の場合に限って話を進めます。
Lagrangeの未定乗数法を用いて、

とおき、q_iで偏微分下値がゼロになるとします。

まとめます。

ここでλからなるベクトルを

とします。行列でまとめると

最後にラムダを求めた後、qを求めます。

但し、

が常に安定して逆行列が求められるわけではありません。 Helsinkiでは特異値分解が好きなようで、よく使われます。

とおいて、今までの式を書き換えます。 尚、nは任意に決めたランクの次元です。

右正中神経電気刺激による体性感覚誘発磁場のデータで検討してみます。
脳内の格子点の設定は、MCEを行った際に得られる
xxx-bem_pointsetXXX.mat
を用いることにします。
stacksize(50000000);
loadmatfile('.../xxx.mat');
loadmatfile('.../xxx-bem_pointsetXX.mat');
getf('.../CalcLFM.sce');
getf('.../CalcSSP2.sce');

Ch=ConversionCoordination(MEG2HEAD,ChInform);
Center=[0,0,40]/1000;
R=CalcTangentialVectors(pointset.points,Center);
wm=1;// gradiometer : magnetometer
LFM=CalcVectorViewLFM(R,Ch,Center,wm);//unstable to provide complex number?
[Usvd,Ssvd,Vsvd]=svd(real(LFM));
n=120;//default
Usvd=Usvd(:,1:n);Ln=Ssvd(1:n,:)*Vsvd';
clear Vsvd;//recommended to save memory
Ssvd=diag(Ssvd);
sum(Ssvd(1:n))/sum(Ssvd)
MN=Ln'*inv(Ln*Ln')*Usvd';
// MEG data, noise data, variance.
meg=CalcSSP2(ave_0001(:,1:306));
meg=meg';time=ave_0001(:,$);
Q=MN*meg;
これで各格子点の電流が求まりました。念のため元波形、計算波形、差分波形を見てみます。
scf();xset('use color',0);
subplot(311);plot2d(time,meg');
set(gca(),'data_bounds',[-50,-700;500,700]);
set(gca(),'tight_limits','on');
subplot(312);plot2d(time,Q'*LFM');
set(gca(),'data_bounds',[-50,-700;500,700]);
set(gca(),'tight_limits','on');
subplot(313);plot2d(time,meg'-Q'*LFM');
set(gca(),'data_bounds',[-50,-700;500,700]);
set(gca(),'tight_limits','on');

上段が元波形、中段が計算波形、下段が差分波形です。 元波形と計算波形が似ていることがわかります。
次に22msecの電流分布をみてみます。
t=min(find(time>=22));
nq=size(Q,1)/2;
q=(Q(1:nq,t).^2+Q(nq+1:$,t).^2).^0.5;//unstable to cause complex number?
q=abs(q)'*1e+9;//nA
scf();set(gcf(),'visible','off');
ncolor=32;
set(gcf(),'color_map',jetcolormap(ncolor));
subplot(121);
P=R.P';
param3d(P(1,:),P(2,:),P(3,:));
h=gca();h.box='off';
h.children.mark_mode='on';h.children.mark_size=1;
h.children.foreground=1;
for nn=1:ncolor-1;...
PP=P;n=1:nq;n(q < nn/ncolor*5)=[];PP=PP(:,n);...
param3d(PP(1,:),PP(2,:),PP(3,:));...
h.children(1).mark_mode='on';h.children(1).mark_size=round(nn/32*5+3);...
h.children(1).foreground=nn+1;...
end;
h.rotation_angles=[70,200];
h.axes_visible='off';h.isoview='on';
subplot(122);Matplot((ncolor:-1:1)');
hh=gca();hh.axes_visible='off';hh.box='off';
hh.tight_limits='on';hh.axes_bounds=[0.75,0.2,0.05,0.6];
h.axes_bounds=[0,0,0.9,1];h.margins=[0,0,0,0];
set(gcf(),'visible','on');

青が0nA以下、赤が5nAです。なんかうまくいきませんでした。

86msecの時の電流源推定の結果です。
小脳を除去して格子点を作成しているので、後頭部のセンサーの磁場を最小電流量で説明しようと、後頭葉底面の格子点に大きな電流源が推定されているようです。

Minimum Current Estimationでは次数30を推奨しています。 30にしてL2ノルムを解くことにします。
[Usvd,Ssvd,Vsvd]=svd(real(LFM));
n=30;//default
Usvd=Usvd(:,1:n);Ln=Ssvd(1:n,:)*Vsvd';
clear Vsvd;//recommended to save memory
Ssvd=diag(Ssvd);
sum(Ssvd(1:n))/sum(Ssvd)
MN=Ln'*inv(Ln*Ln')*Usvd';
Q=MN*meg;
元波形、計算波形、差分波形を見てみます。
scf();xset('use color',0);
subplot(311);plot2d(time,meg');
set(gca(),'data_bounds',[-50,-700;500,700]);
set(gca(),'tight_limits','on');
subplot(312);plot2d(time,Q'*LFM');
set(gca(),'data_bounds',[-50,-700;500,700]);
set(gca(),'tight_limits','on');
subplot(313);plot2d(time,meg'-Q'*LFM');
set(gca(),'data_bounds',[-50,-700;500,700]);
set(gca(),'tight_limits','on');

上段が元波形、中段が計算波形、下段が差分波形です。 次数120と次数30で差分波形はあまり変わりなさそうです。
次に22msecと86msecの電流分布をみてみます。
t=min(find(time>=22));
//t=min(find(time>=86));
nq=size(Q,1)/2;
q=(Q(1:nq,t).^2+Q(nq+1:$,t).^2).^0.5;//unstable to cause complex number?
q=abs(q)'*1e+9;//nA
scf();set(gcf(),'visible','off');
ncolor=32;
set(gcf(),'color_map',jetcolormap(ncolor));
subplot(121);
P=R.P';
param3d(P(1,:),P(2,:),P(3,:));
h=gca();h.box='off';
h.children.mark_mode='on';h.children.mark_size=1;
h.children.foreground=1;
for nn=1:ncolor-1;...
PP=P;n=1:nq;n(q < nn/ncolor*0.2)=[];PP=PP(:,n);...
param3d(PP(1,:),PP(2,:),PP(3,:));...
h.children(1).mark_mode='on';h.children(1).mark_size=round(nn/32*5+3);...
h.children(1).foreground=nn+1;...
end;
h.rotation_angles=[70,200];
h.axes_visible='off';h.isoview='on';
subplot(122);Matplot((ncolor:-1:1)');
hh=gca();hh.axes_visible='off';hh.box='off';
hh.tight_limits='on';hh.axes_bounds=[0.75,0.2,0.05,0.6];
h.axes_bounds=[0,0,0.9,1];h.margins=[0,0,0,0];
set(gcf(),'visible','on');

青が0nA以下、赤が0.2nAです。なんかいい感じです。

86msecの時の電流源推定の結果です。
次数をどのくらいに設定するかは経験則によるのだと思います。