lcd_fitl2reg_holocalAR

 MAR parameter estimate with l2 reguralization for higher order local AR
 model (holocalAR model)

 * higher order local AR model 
  j(v,t) = sum_p a(v,p)*j(v,t-p) + sum_u b(v,u)*j(u,t-Delta(v,u)) + epsilon(t)
 
 [Input] 
 Zact :  current waveforms [Nv Nt]
 Delta : delay matrix (sparse) with unit of time points 
         (i.e. all the elements must be integer representing dicrete time
         points)
   p    : AR order for local dynamics 
 lambda : a regularization parameter common to all the vertices
          if set [], optimal lambda is computed for each vertex  
 localfactor : if specified, local regularization parameter is lambda*localfactor 
 
 [Output]
  A  : diagonal parts of MAR matrix (local dynamics)  [Nv,p]
  B  : non-digonal parts of MAR matrix (distant interactions) [Nv Nv]

 2016/05/09 O.Yamashita
 * avoid increasing memory allocation
 * allow for multiple lambdas
 2016/04/11 O.Yamashita 
 2016/01/29 O.Yamashita

0001 function [A,B,Lam,Sig]=lcd_fitl2reg_holocalAR(Zact,Delta,p,lambda, localfactor)
0002 % MAR parameter estimate with l2 reguralization for higher order local AR
0003 % model (holocalAR model)
0004 %
0005 % * higher order local AR model
0006 %  j(v,t) = sum_p a(v,p)*j(v,t-p) + sum_u b(v,u)*j(u,t-Delta(v,u)) + epsilon(t)
0007 %
0008 % [Input]
0009 % Zact :  current waveforms [Nv Nt]
0010 % Delta : delay matrix (sparse) with unit of time points
0011 %         (i.e. all the elements must be integer representing dicrete time
0012 %         points)
0013 %   p    : AR order for local dynamics
0014 % lambda : a regularization parameter common to all the vertices
0015 %          if set [], optimal lambda is computed for each vertex
0016 % localfactor : if specified, local regularization parameter is lambda*localfactor
0017 %
0018 % [Output]
0019 %  A  : diagonal parts of MAR matrix (local dynamics)  [Nv,p]
0020 %  B  : non-digonal parts of MAR matrix (distant interactions) [Nv Nv]
0021 %
0022 % 2016/05/09 O.Yamashita
0023 % * avoid increasing memory allocation
0024 % * allow for multiple lambdas
0025 % 2016/04/11 O.Yamashita
0026 % 2016/01/29 O.Yamashita
0027 
0028 if nargin < 5
0029     localfactor = [];
0030 end
0031 
0032 
0033 [Nv,Nt] = size(Zact);
0034 Dmax = max(Delta(:));
0035 Dmax = max([Dmax,p]);
0036 
0037 if isempty(lambda) | length(lambda) == 1
0038     B = zeros(Nv,Nv);
0039     A = zeros(Nv,p);
0040 else
0041     B = zeros(Nv,Nv,length(lambda));
0042     A = zeros(Nv,p,length(lambda));
0043 end
0044 
0045 
0046 for vv = 1 : Nv
0047     
0048     fprintf('process %04d ...\n',vv);
0049     
0050     y = [Zact(vv,Dmax+1:Nt)]';
0051     
0052     dvec = round(Delta(vv,:));
0053     iv = setdiff(find(dvec ~= 0), vv); % other than self loop
0054     
0055     % define X
0056     X = zeros(Nt-Dmax, length(iv)+p);
0057     for aa = 1 : p
0058     X(:,aa) = Zact(vv,Dmax+1-aa:Nt-aa)';
0059     end
0060     nn = p+1;
0061     for uu = iv
0062         X(:,nn) = Zact(uu,Dmax+1-dvec(uu):Nt-dvec(uu))';
0063         nn = nn + 1;
0064     end
0065     
0066     % weighted reguralization for local model
0067     w = ones(size(X,2),1);
0068     if ~isempty(localfactor)
0069         w(1:p) = localfactor;
0070     end
0071     
0072     
0073     if isempty(lambda) | length(lambda) == 1
0074         [lam,sig,val,a] = wmn_solution(y,X,[],[],struct('flag_plot',0,'lambda',lambda));
0075         A(vv,:) = a(1:p);    
0076         B(vv,iv) = a(p+1:end);
0077         Lam(vv) = lam;
0078         Sig(vv) = sig;
0079     else
0080         [aa] = wmn_solution_fixreg(y, X, [], [], lambda, struct('flag_plot',0));
0081         A(vv,:,:) = aa(1:p,1,:);
0082         B(vv,iv,:) = aa(p+1:end,1,:);
0083         Lam = lambda;
0084         Sig = [];
0085     end
0086     
0087     
0088     
0089     
0090     
0091 end
0092 
0093  
0094 
0095 
0096 
0097 
0098 
0099