-
Notifications
You must be signed in to change notification settings - Fork 0
/
hcp_old_obj.m
83 lines (75 loc) · 1.78 KB
/
hcp_old_obj.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
%hcp_old_obj old method of computing HCP objective and derivatives
%
% Objective:
% f = -det(I-P+alpha*e*e')
%
% Gradient:
% G1 = I-P+alpha*e*e'
% G2 = -f*inv(G)'
% df/dp_{ij} = G2(i,j)
%
% Hessian:
% Note that df/dp_{ij} = -(-1)^(i+j)*det(G(r,c))
% G(r,c) is G with row i and col j removed
% Thus, grad(df/dp_{ij}) = -(-1)^(i_j)*grad(det(G(r,c))
% This code uses a recursive call to compute these gradients
%
% Note: there are different ways to compute all of these terms. This is
% not the most efficient way.
%
% Inputs:
% x = variable vector
% A = adjacency matrix
% alpha = parameter
% rows = logical index for included rows (optional)
% cols = logical index for included cols (optional)
%
% Outputs:
% f = objective function value
% g = gradient
% H = Hessian
%
% References:
% Matrix Reference Manual
% http://www.ee.ic.ac.uk/hp/staff/dmb/matrix/intro.html
% Matrix Cookbook
% http://matrixcookbook.com/ (down at the moment)
%
% 2010-03-26 (nwh) first working version
function [f g H] = hcp_old_obj(x,A,alpha,rows,cols)
n = size(A,1);
if nargin < 4
rows = true(n,1);
cols = true(n,1);
end
inds = (A == 1);
indr = (A(rows,cols) == 1);
P = double(A);
P(inds) = x;
m = length(x);
G = eye(n)-P+alpha*ones(n,n);
G = G(rows,cols);
% compute the objective function
f = -det(G);
if nargout > 1
% compute the gradient
B = -(f*inv(G))';
g = B(indr);
end
if nargout > 2
% compute the Hessian with recursive call
[i j] = find(A);
r = true(n,1);
c = true(n,1);
H = zeros(n,n);
for k = 1:m
r(i(k)) = false;
c(j(k)) = false;
v = (i~=i(k) & j~=j(k));
[f2 g2] = hcp_old_obj(x,A,alpha,r,c);
H(k,v) = -(-1)^(i(k)+j(k))*g2;
r(i(k)) = true;
c(j(k)) = true;
end
end
end