% Hamza Fawzi
% October 25th, 2012
% Last updated: November 25th, 2012 (thanks to Johan Lofberg, who suggested
% using Yalmip's Automatic Dualization to make the case k=0 run faster)
% -------------------------------------------------------------------------
% This script computes a lower bound to the nonnegative rank of a matrix
% using the method described in the paper:
%
%   New lower bounds on nonnegative rank using conic programming
%   Hamza Fawzi, Pablo A. Parrilo
%   http://arxiv.org/abs/1210.6970
%
% The lower bound relies on the quantity (the input matrix is A)
%   \nu_+(A) = max_{W}  trace(A*W)   subject to   [I -W; -W^T I] copositive
% The lower bound on the nonnegative rank is then
%   \rank_+(A) >= (\nu_+(A) / ||A||_F)^2
% where ||A||_F is the Frobenius norm of A.
%
% This script computes approximations from below to the quantity \nu_+(A)
% using semidefinite-programming based relaxations of the copositive cone
% (see the paper for more details).
% The desired hierarchy level is controlled by the parameter k below.
% 
% This script uses the package YALMIP (and its Sum-Of-Squares module), 
% which can be downloaded from: http://users.isy.liu.se/johanl/yalmip/
% -------------------------------------------------------------------------

 A = [1 1 0 0;
      1 0 1 0;
      0 1 0 1;
      0 0 1 1];

k = 0; % SOS hierarchy level for copositive cone

m = size(A,1);
n = size(A,2);

if k == 0,
    % Special case for relaxation of order k=0
    W = sdpvar(m,n,'full');
    N = sdpvar(m+n,m+n,'symmetric');
    P = sdpvar(m+n,m+n,'symmetric');
    Objective = -trace(A'*W);
    Constraints = [[eye(m,m) -W; -W' eye(n,n)] == P+N, P >= 0, N(:) >= 0];
    solvesdp(Constraints,Objective,sdpsettings('verbose',1,'dualize',1));
else
    % Solve k'th relaxation of copositive program
    W = sdpvar(m,n,'full');
    x = sdpvar(m+n,1);
    pol = (x'*x)^k*(x.*x)'*[eye(m,m) -W; -W' eye(n,n)]*(x.*x);
    Objective = -trace(A'*W);
    Constraints = set(sos(pol));
    solvesos(Constraints,Objective,sdpsettings('verbose',1),W(:));
end
fprintf('Finished solving, optval = %.4f\n', double(-Objective));

fprintf('----------------------------\n');
fprintf('lower bound = %.2f\n', (double(-Objective) / norm(A,'fro'))^2);
fprintf('trivial lower bound from rank(A) = %d\n', rank(A));