% 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));