% Function to calculate the Hessian of lambda to the entries of projection
% matrix A.
%
% Input:
%    A = the projection matrix (n by n matrix)
%
% Function outputs result of Equation (24) (n^2 by n^2 matrix).
%

function H = Hlambda_A(A)

[lambda, w, v] = domeig(A);
n = length(v); 
In = eye(n);
et = ones(n, 1)';

dw = inv(lambda*In - A + w*et*A) * kron(w', In-w*et); % Eq (19)
dv = inv(lambda*In - A' + lambda*v*w') * (kron(In-v*w', v') - lambda*kron(v,v')*dw); % Eq (20)

B = kron(In, v)*dw + kron(w, In)*dv; % Eq (23)
H = 1/2*(B+B'); % Eq (24)

end