SHOW_OPTIMAL shows the optimal excitatory stimuli of an SFA object.
SHOW_OPTIMAL(HDL, NRM, TOL) Show the optimal excitatory stimuli of SFA
object HDL, under a fixed norm constraint. The fixed norm is specified
by NRM. TOL gives the maximal norm difference tollerated by the
maximization algorithm.
Optional arguments can be specified by appending
'ArgumentName',ArgumentValue pairs to the argument list
(e.g. SHOW_OPTIMAL(HDL,NRM,TOL,'start',5,'show_xm',1) ).
Possible optional arguments:
'show_xm' (default:0) if set to 1 shows the optimal inhibitory stimuli
'h' (default:16) height of the input patch
'w' (default:16) width of the input patch
'sh' (default:7)
'sw' (default:7) the optimal stimuli are displayed on a SH times SW
grid
'start' (default:1) the first unit to consider0001 function show_optimal(hdl, nrm, tol, varargin), 0002 % SHOW_OPTIMAL shows the optimal excitatory stimuli of an SFA object. 0003 % 0004 % SHOW_OPTIMAL(HDL, NRM, TOL) Show the optimal excitatory stimuli of SFA 0005 % object HDL, under a fixed norm constraint. The fixed norm is specified 0006 % by NRM. TOL gives the maximal norm difference tollerated by the 0007 % maximization algorithm. 0008 % 0009 % Optional arguments can be specified by appending 0010 % 'ArgumentName',ArgumentValue pairs to the argument list 0011 % (e.g. SHOW_OPTIMAL(HDL,NRM,TOL,'start',5,'show_xm',1) ). 0012 % 0013 % Possible optional arguments: 0014 % 'show_xm' (default:0) if set to 1 shows the optimal inhibitory stimuli 0015 % 0016 % 'h' (default:16) height of the input patch 0017 % 'w' (default:16) width of the input patch 0018 % 0019 % 'sh' (default:7) 0020 % 'sw' (default:7) the optimal stimuli are displayed on a SH times SW 0021 % grid 0022 % 0023 % 'start' (default:1) the first unit to consider 0024 0025 0026 %%%% default values 0027 0028 % set to 1 if you want to see the optimal inhibitory stimuli 0029 ctxt.show_xm = 0; 0030 0031 % input patch height and width 0032 ctxt.h = 16; ctxt.w = 16; 0033 0034 % number of optimal stimuli to display (vertically and horizontally) 0035 ctxt.sh = 7; ctxt.sw = 7; 0036 % first unit to consider 0037 ctxt.start = 1; 0038 % default window position 0039 ctxt.Position = [360 372 597 562]; % NxN 0040 0041 % overwrite with user-defined list of settings 0042 for k = 1:2:length(varargin); 0043 % error check: the optional arguments must be defined as name-value pairs 0044 if ~ischar(varargin{k}), 0045 error 'Setting names must be strings'; 0046 end 0047 % set variable value 0048 ctxt=setfield(ctxt,varargin{k},varargin{k+1}); 0049 end 0050 0051 %%%%%%%%%%%%%%%%%%%%%%%%%%%% 0052 0053 % range of units to consider 0054 range = ctxt.start:(ctxt.start+ctxt.sh*ctxt.sw-1); 0055 0056 % copy some useful quantities 0057 global SFA_STRUCTS 0058 avg = SFA_STRUCTS{hdl}.avg0; 0059 DW = SFA_STRUCTS{hdl}.DW0; 0060 D = SFA_STRUCTS{hdl}.D0; 0061 PCinv = DW*diag(D.^(-0.5)); 0062 0063 % set the figure up 0064 clf; colormap(gray); 0065 set(gcf, 'Position', ctxt.Position); 0066 0067 % subplot counter 0068 k=1; 0069 % loop over all subunits 0070 for i=range, 0071 % get the quadratic form corresponding to the current subunit in the 0072 % principal components space (much faster then in the input space). 0073 % ! technical detail: this is equivalent to finding the optimal 0074 % stimuli in the input _meanfree_ space. 0075 % This is similar to physiological experiments where the optimal stimuli 0076 % are defined as intensity changes with respect to a constant mean 0077 % luminance (and thus have negative and positive intensity values). 0078 [H,f,c] = sfa_getHf(hdl, i, 1); 0079 0080 % get the maximum and the minimum of the quadratic form 0081 xp = maximize_qform(H,f, [], nrm, tol); 0082 vp = 0.5*xp'*H*xp + f'*xp + c; 0083 xm = maximize_qform(-H,-f, [], nrm, tol); 0084 vm = 0.5*xm'*H*xm + f'*xm + c; 0085 0086 % xp has the largest value (or xm if show_xm==1) 0087 if abs(vm)>vp & ~ctxt.show_xm, 0088 tmp = xp; xp = xm; xm = tmp; 0089 end; 0090 0091 % project xp back to the input space 0092 xp = PCinv*xp; 0093 xm = PCinv*xm; 0094 0095 % plot 0096 subplot(ctxt.sh,ctxt.sw,k); 0097 imagesc(reshape(xp, ctxt.h, ctxt.w)); 0098 axis off; axis image; drawnow, 0099 0100 % increase the subplot counter 0101 k = k+1; 0102 end 0103 0104 0105 function x = maximize_qform(H,f, x0, nrm, tol), 0106 % maximize the quadratic form 1/2*x'*H*x + f'*x + c 0107 0108 % center the quadratic form around x0 0109 if ~isempty(x0), 0110 f = H*x0+f; 0111 %c = 0.5*x0'*H*x0 + f'*x0 + c; 0112 end 0113 0114 % input dimension 0115 dim = size(H,1); 0116 % norm of f 0117 nrm_f = norm(f); 0118 0119 % get eigenvalues and eigenvectors 0120 [V,D] = eig(H); 0121 0122 mu = diag(D)'; 0123 % coefficients of the eigenvectors decomposition of f 0124 alpha = V'*f; 0125 % v_i = alpha_i * v_i 0126 V=V.*repmat(alpha',dim,1); 0127 0128 % left bound for lambda 0129 % added 'real' to avoid numerical problems if you maximize in input space 0130 ll = max(real(mu)); 0131 % right bound for lambda 0132 lr = norm(f)/nrm + ll; 0133 0134 % search by bisection until norm(x)^2 = nrm^2 0135 nrm_2 = nrm^2; 0136 norm_x_2 = 0; 0137 while abs(norm_x_2-nrm_2)>tol, 0138 % bisection of the lambda-interval 0139 lambda=(lr-ll)/2+ll; 0140 % eigenvalues of (lambda*Id - H)^-1 0141 beta = (lambda-mu).^(-1); 0142 0143 % solution to the second lagragian equation 0144 norm_x_2 = sum(alpha.^2.*beta'.^2); 0145 0146 %[ll,lr] 0147 if norm_x_2>nrm_2, ll=lambda; 0148 else lr=lambda; end 0149 %[ll, lr, norm_x_2] 0150 %pause(1) 0151 end 0152 0153 x = sum(V.*repmat(beta,dim,1),2); 0154 0155 if ~isempty(x0), 0156 x=x+x0; 0157 end