% Function to calculate the Hessian of r 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 (41) (s by s matrix).
%

function H = Hr_theta(A, DAt, HAt)

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

H = -1/lambda^2*DAt'*kron(w*w', v*v')*DAt + 1/lambda*Hlambda_theta(A, DAt, HAt); % Eq (41)

end