function E = sdp(J)

N = size(J,1);
Z = sdpvar(N,N);

% Objective
O = trace(J*Z);

% Constrains
C = set(Z>=0) + set(diag(Z) == 1);

% Solve
solvesdp(C, O);
X = double(Z);

[Cp,Op,free] = primalize(C,-O);
solvesdp(Cp,Op)
setsdpvar(free,dual(Cp(end)));
xp = double(free);

x = X(1,:);
E = x*J*x';