OPTIMAL_STIMULUS computes the optimal stimulus of a quadratic form. X = OPTIMAL_STIMULUS(H,F, R, EPS) computes the optimal excitatory stimulus of the quadratic form 0.5*x'*H*x + f'*x + c, i.e. the input vector X that maximizes the quadratic form given a fixed norm R. EPS tolerance of norm(X) from R. This function can be used to compute the optimal inhibitory stimulus by calling it with the negative of the quadratic form -H, -F . Reference: Pietro Berkes and Laurenz Wiskott (2005) On the analysis and interpretation of quadratic forms as receptive fields.
0001 function x = optimal_stimulus(H,f, r, eps), 0002 % OPTIMAL_STIMULUS computes the optimal stimulus of a quadratic form. 0003 % X = OPTIMAL_STIMULUS(H,F, R, EPS) computes the optimal excitatory stimulus 0004 % of the quadratic form 0.5*x'*H*x + f'*x + c, i.e. the input vector X 0005 % that maximizes the quadratic form given a fixed norm R. 0006 % 0007 % EPS tolerance of norm(X) from R. 0008 % 0009 % This function can be used to compute the optimal inhibitory stimulus 0010 % by calling it with the negative of the quadratic form -H, -F . 0011 % 0012 % Reference: Pietro Berkes and Laurenz Wiskott (2005) On the analysis 0013 % and interpretation of quadratic forms as receptive fields. 0014 0015 % input dimension 0016 dim = size(H,1); 0017 0018 % compute the eigenvalues mu and eigenvectors V 0019 [V,D] = eig(H); 0020 mu = diag(D)'; 0021 % compute the coefficients of the eigenvectors decomposition of f 0022 alpha = V'*f; 0023 0024 % compute the range of the parameter lambda 0025 % left bound for lambda 0026 % added 'real' to avoid numerical problems if you maximize a 0027 % ill-conditioned quadraitc form 0028 lambda_left = max(real(mu)); 0029 % right bound for lambda 0030 lambda_right = norm(f)/r + lambda_left; 0031 0032 % search by bisection until norm(x)^2 = r^2 0033 r_2 = r^2; 0034 % norm_x_2 holds the value of norm(x)^2 at the current lambda 0035 norm_x_2 = 0; 0036 while abs(sqrt(norm_x_2)-r) > eps, 0037 % bisect the lambda interval 0038 lambda = (lambda_right-lambda_left)/2 + lambda_left; 0039 % compute the eigenvalues of (lambda*Id - H)^-1 0040 beta = (lambda-mu).^(-1); 0041 0042 % compute norm(x)^2 at lambda 0043 norm_x_2 = sum(beta'.^2.*alpha.^2); 0044 0045 % update the interval limits 0046 if norm_x_2 > r_2, 0047 lambda_left = lambda; 0048 else 0049 lambda_right = lambda; 0050 end 0051 end 0052 0053 % lambda found, compute the solution 0054 x = sum(repmat(beta,dim,1).*V.*repmat(alpha',dim,1),2);