function t = tauplus_sos(A) % t = TAUPLUS_SOS(A) % % Computes the quantity \tau_+^{sos}(A) as defined in the paper % % "Self-scaled bounds for atomic cone ranks: % applications to nonnegative rank and cp-rank" % Hamza Fawzi, Pablo A. Parrilo % http://arxiv.org/abs/1404.3240 % % This is a lower bound on the nonnegative rank of A. % % INPUT % A: a nonnegative matrix % % OUTPUT % t: a real number equal to \tau_+^{sos}(A) % % EXAMPLE % >> A = [1 1 0 0; 1 0 1 0; 0 1 0 1; 0 0 1 1]; % >> tauplus_sos(A) % Outputs 4 % % NOTE % Requires Yalmip (http://users.isy.liu.se/johanl/yalmip/) % % AUTHOR % Hamza Fawzi (hfawzi@mit.edu) [m,n] = size(A); if any(A(:) < 0) error('A has to be a nonnegative matrix'); end % Extract strictly positive elements of A ix = find(A > 0); a = A(ix); N = numel(ix); IX = zeros(size(A)); IX(ix) = 1:N; % Setup SDP using Yalmip X = sdpvar(N,N,'symmetric'); t = sdpvar(1); p_cons = [[t a'; a X] >= 0; diag(X) <= a.^2]; % Setup the linear equalities X_{ij,kl} = X_{il,kj} corresponding % to the rank-one constraint (vanishing of 2x2 minors) for ii=1:m, for jj=1:n, for kk=(ii+1):m, for ll=(jj+1):n, ind_ij = IX(ii,jj); ind_kl = IX(kk,ll); ind_il = IX(ii,ll); ind_kj = IX(kk,jj); if A(ii,jj)*A(kk,ll) > 0 && A(ii,ll)*A(kk,jj) == 0 p_cons = [p_cons; X(ind_ij,ind_kl) == 0]; elseif A(ii,jj)*A(kk,ll) == 0 && A(ii,ll)*A(kk,jj) > 0 p_cons = [p_cons; X(ind_il,ind_kj) == 0]; elseif A(ii,jj)*A(kk,ll) > 0 && A(ii,ll)*A(kk,jj) > 0 p_cons = [p_cons; X(ind_ij,ind_kl) == X(ind_il,ind_kj)]; end end end end end p_obj = t; p_sol = solvesdp(p_cons,p_obj); t = double(t); end