Spotless isos

spotopt/@spotprog/spotprog.m

@@ -623,23 +623,57 @@
             if ~spotprog.isScalarDimension(dim)
                 error('dim must be a row of non-negative scalars.');
             end
+
             ndim = spotprog.psdDimToNo(dim);
 
             [pr,Qvec] = pr.newFree(ndim);
             Q = mss_v2s(Qvec);
-
-
-           % Constrain Q to be in dual of DD (non-vectorized version)
-           n = length(Q);
-           pr = pr.withPos(diag(Q));
-           for i = 1:n
-               for j = [1:i-1, i+1:n] % j \neq i
-                   [pr,tau(i,j)] = pr.newFree(1);
-                   pr = pr.withPos(tau(i,j) - Q(i,j));
-                   pr = pr.withPos(tau(i,j) + Q(i,j));
-                   pr = pr.withPos(Q(i,i) + Q(j,j) - tau(i,j));
-               end
-           end
+
+            n = dim;
+            % Vectors with at most two non-zero entries with +/-1.
+            V = zeros(n,n^2); % preallocate for speed
+            V(end-n+1:end,end-n+1:end)=eye(n);
+            cnt = 0;
+            for i=1:n-1
+                for j=i+1:n
+                    cnt = cnt+1;
+                    v1=zeros(n,1); v1(i)=1; v1(j)=1;
+                    v2=zeros(n,1); v2(i)=1; v2(j)=-1;
+                    V(:,2*cnt-1) = v1;
+                    V(:,2*cnt) = v2;
+                end
+            end
+
+            % Setup the constraints in chunks (if we do it all at once it will eat up
+            % a huge amount of RAM).
+            chunkSize = 1000;
+            numChunks = ceil((n^2)/chunkSize);
+            chunkInds = 1 + [0:chunkSize:(numChunks-1)*chunkSize];
+            chunkInds = [chunkInds, n^2+1];
+
+            for chunk = 1:length(chunkInds)-1
+                startInd = chunkInds(chunk);
+                endInd = chunkInds(chunk+1)-1;
+
+                % Get chunk of V
+                V_chunk = V(:,startInd:endInd);
+
+                % Setup constraint that V(:,k)'*Q*V(:,k) >=0, for k in chunk in vectorized way
+                pr = pr.withPos(sum(V_chunk.*(Q*V_chunk)));
+            end
+
+
+%            % Constrain Q to be in dual of DD (non-vectorized version)
+%            n = length(Q);
+%            pr = pr.withPos(diag(Q));
+%            for i = 1:n
+%                for j = [1:i-1, i+1:n] % j \neq i
+%                    [pr,tau(i,j)] = pr.newFree(1);
+%                    pr = pr.withPos(tau(i,j) - Q(i,j));
+%                    pr = pr.withPos(tau(i,j) + Q(i,j));
+%                    pr = pr.withPos(Q(i,i) + Q(j,j) - 2*tau(i,j));
+%                end
+%            end
 
 
         end

util/isSOS.m

@@ -0,0 +1,35 @@
+function [issos, Q, phi] = isSOS(p)
+
+x = decomp(p);
+
+prog = spotsosprog;
+prog = prog.withIndeterminate(x);
+
+[prog,psos] = prog.newFreePoly(monomials(x,0:deg(p,x)));
+[prog,ind1] = prog.withSOS(p);
+
+prog = prog.withPolyEqs(p - psos);
+
+
+% Options
+options = spot_sdp_default_options();
+options.verbose = 2;
+options.solveroptions.MSK_IPAR_BI_CLEAN_OPTIMIZER = 'MSK_OPTIMIZER_INTPNT'; % Use just the interior point algorithm to clean up
+options.solveroptions.MSK_IPAR_INTPNT_BASIS = 'MSK_BI_NEVER'; % Don't use basis identification (it's slow)
+options.scale_monomials = false;
+
+
+
+sol = prog.minimize(0, @spot_mosek, options);
+
+if (sol.isPrimalFeasible) && (sol.isDualFeasible)
+    issos = true;
+    Q = sol.eval(sol.gramMatrices{ind1});
+    phi = sol.gramMonomials{ind1};
+else
+    issos = false;
+    Q = [];
+    phi = [];
+end
+
+end