% Function to calculate the Hessian of lambda to the s by 1 parameter
% vector theta.
%
% Input:
%    A = the projection matrix (n by n matrix)
%    s = length of the parameter vector theta
%    DAt = D[vecA;theta] (n^2 by s matrix)
%    HAt = H[vecA;theta] (n^2*s by s matrix)
%
% Function outputs result of Equation (35) (s by s matrix).
%

function H = Hlambda_theta(A, DAt, HAt)

[lambda, w, v] = domeig(A);
[n2,s] = size(DAt);
Is = eye(s);

B = kron(kron(w', v'), Is)*HAt + DAt'*Hlambda_A(A)*DAt; % Eq (34)
H = 1/2*(B+B'); % Eq (35)

end