function t = tauplus_sos_tensor(A) % t = TAUPLUS_SOS_TENSOR(A) % % Computes the quantity \tau_+^{sos}(A) for a 3-dimensional tensor 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 3-dimensional nonnegative array % % OUTPUT % t: a real number equal to \tau_+^{sos}(A) % % EXAMPLE % >> A = zeros(2,2,2); A(:,:,1) = [1 1; 1 0]; A(:,:,2) = [0 1; 1 1]; % >> tauplus_sos_tensor(A) % Outputs 3 % % NOTE % Requires Yalmip (http://users.isy.liu.se/johanl/yalmip/) % % AUTHOR % Hamza Fawzi (hfawzi@mit.edu) if any(A(:) < 0) error('A has to be nonnegative'); end if numel(size(A)) ~= 3, error('A has to be a 3-dimensional array'); end a = A(:); [n1,n2,n3] = size(A); N = n1*n1*n3; X = sdpvar(N,N,'symmetric'); t = sdpvar(1); p_cons = [[t a'; a X] >= 0; diag(X) <= a.^2]; % Setup the linear equalities corresponding to the rank-one conditions for % 3-dimensional tensors s2i = @(x,y,z) sub2ind(size(A),x,y,z); if n1 == 2 && n2 == 2 && n3 == 2, % Use hard-coded equations for the 2x2x2 case (9 equations in total) p_cons = [ p_cons; % 2x2 minor in slice (:,:,1) X(s2i(1,1,1),s2i(2,2,1)) == X(s2i(1,2,1),s2i(2,1,1)); % 2x2 minor in slice (:,:,2) X(s2i(1,1,2),s2i(2,2,2)) == X(s2i(1,2,2),s2i(2,1,2)); % 2x2 minor in slice (:,1,:) X(s2i(1,1,1),s2i(2,1,2)) == X(s2i(1,1,2),s2i(2,1,1)); % 2x2 minor in slice (:,2,:) X(s2i(1,2,1),s2i(2,2,2)) == X(s2i(1,2,2),s2i(2,2,1)); % 2x2 minor in slice (1,:,:) X(s2i(1,1,1),s2i(1,2,2)) == X(s2i(1,1,2),s2i(1,2,1)); % 2x2 minor in slice (2,:,:) X(s2i(2,1,1),s2i(2,2,2)) == X(s2i(2,1,2),s2i(2,2,1)); % 2x2 minors between (1,1,1) and (2,2,2) X(s2i(1,1,1),s2i(2,2,2)) == X(s2i(1,1,2),s2i(2,2,1)); X(s2i(1,1,1),s2i(2,2,2)) == X(s2i(1,2,1),s2i(2,1,2)); X(s2i(1,1,1),s2i(2,2,2)) == X(s2i(2,1,1),s2i(1,2,2)); ]; else % Generate equations in the general case % Note: this code is not optimized and there may be redundant equations lst = [0 0 0 0]; % Loop over multi-indices alpha and beta for ind_alpha=1:N, for ind_beta=(ind_alpha+1):N, alpha = cell(1,3); beta = cell(1,3); [alpha{:}] = ind2sub(size(A),ind_alpha); [beta{:}] = ind2sub(size(A),ind_beta); alpha = [alpha{:}]; beta = [beta{:}]; for k=1:3, if alpha(k) ~= beta(k) % Swap index k of alpha and beta alpha1 = alpha; alpha1(k) = beta(k); beta1 = beta; beta1(k) = alpha(k); ind_alpha1 = s2i(alpha1(1),alpha1(2),alpha1(3)); ind_beta1 = s2i(beta1(1),beta1(2),beta1(3)); sig = [ind_alpha ind_beta ... min(ind_alpha1,ind_beta1) ... max(ind_alpha1,ind_beta1)]; if ~ismember(sig,lst,'rows') lst = union(lst,sig,'rows'); % Add constraint p_cons = [p_cons; X(ind_alpha,ind_beta) == ... X(ind_alpha1,ind_beta1)]; end end end end end fprintf('# of equations added = %d\n', size(lst,1)-1); end p_obj = t; p_sol = solvesdp(p_cons,p_obj); t = double(t); end