invariances

 INVARIANCES computes the invariances of a quadratic form.
   [V,D] = INVARIANCES(H,F, XP) computes the invariances V of the
   quadratic form 0.5*x'*H*x + f'*x + c in the optimal stimulus XP.
   Return in V the direction of the invariances, sorted by increasing
   magnitude of the second derivative and in D the corresponding second
   derivative.

   Reference: Pietro Berkes and Laurenz Wiskott (2005) On the analysis
   and interpretation of quadratic forms as receptive fields.

0001 function [V,d] = invariances(H,f, xp),
0002 % INVARIANCES computes the invariances of a quadratic form.
0003 %   [V,D] = INVARIANCES(H,F, XP) computes the invariances V of the
0004 %   quadratic form 0.5*x'*H*x + f'*x + c in the optimal stimulus XP.
0005 %   Return in V the direction of the invariances, sorted by increasing
0006 %   magnitude of the second derivative and in D the corresponding second
0007 %   derivative.
0008 %
0009 %   Reference: Pietro Berkes and Laurenz Wiskott (2005) On the analysis
0010 %   and interpretation of quadratic forms as receptive fields.
0011 
0012   % input dimension
0013   dim = size(H,1);
0014 
0015   % radius of the sphere
0016   r = norm(xp);
0017   % constant term in the second derivative
0018   cte = -r^(-2)*( xp'*H*xp + f'*xp );
0019 
0020   %% get a basis of the tangential plane in xp
0021   one = eye(dim);
0022   basis = [xp, one(:,2:dim)];
0023   B = gram_schmidt(basis);
0024   B = B(:,2:dim);
0025   % restrict H to the tangential plane
0026   Ht = B'*H*B;
0027   [V,D] = eig(Ht);
0028   % second derivative in the direction of the eigenvectors
0029   d = diag(D)+cte;
0030   [tmp,idx] = sort(abs(d));
0031   % sort by increasing absolute value of the second derivative
0032   d = d(idx);
0033   V = V(:,idx);
0034   % project the eigenvectors back in R^N
0035   V = B*V;