Matlab for Convex Optimization & Euclidean Distance Geometry

From Wikimization

Revision as of 18:03, 30 July 2008

Made by The MathWorks http://www.mathworks.com, Matlab is a high level programming language and graphical user interface for linear algebra.

isedm()

% Is real D a Euclidean Distance Matrix. -Jon Dattorro
%
% [Dclosest,X,isisnot,r] = isedm(D,tolerance,verbose,dimension,V)
%
% Returns: closest EDM in Schoenberg sense (default output),
%          a generating list X,
%          string 'is' or 'isnot' EDM,
%          actual affine dimension r of EDM output.
% Input: matrix D,
%        optional absolute numerical tolerance for EDM determination,
%        optional verbosity 'on' or 'off',
%        optional desired affine dim of generating list X output,
%        optional choice of 'Vn' auxiliary matrix (default) or 'V'.

function [Dclosest,X,isisnot,r] = isedm(D,tolerance_in,verbose,dim,V);

isisnot = 'is';
N = length(D);

if nargin < 2 | isempty(tolerance_in)
   tolerance_in = eps;
end
tolerance = max(tolerance_in, eps*N*norm(D));
if nargin < 3 | isempty(verbose)
   verbose = 'on';
end
if nargin < 5 | isempty(V)
   use = 'Vn';
else
   use = 'V';
end

% is empty
if N < 1
   if strcmp(verbose,'on'), disp('Input D is empty.'), end
   X = [ ];
   Dclosest = [ ];
   isisnot = 'isnot';
   r = [ ];
   return
end
% is square
if size(D,1) ~= size(D,2)
   if strcmp(verbose,'on'), disp('An EDM must be square.'), end
   X = [ ];
   Dclosest = [ ];
   isisnot = 'isnot';
   r = [ ];
   return
end
% is real
if ~isreal(D)
   if strcmp(verbose,'on'), disp('Because an EDM is real,'), end
   isisnot = 'isnot';
   D = real(D);
end

% is nonnegative
if sum(sum(chop(D,tolerance) < 0))
   isisnot = 'isnot';
   if strcmp(verbose,'on'), disp('Because an EDM is nonnegative,'),end
end
% is symmetric
if sum(sum(abs(chop((D - D')/2,tolerance)) > 0))
   isisnot = 'isnot';
   if strcmp(verbose,'on'), disp('Because an EDM is symmetric,'), end
   D = (D + D')/2;  % only required condition
end
% has zero diagonal
if sum(abs(diag(chop(D,tolerance))) > 0)
   isisnot = 'isnot';
   if strcmp(verbose,'on')
      disp('Because an EDM has zero main diagonal,')
   end
end
% is EDM
if strcmp(use,'Vn')
   VDV = -Vn(N)'*D*Vn(N);
else
   VDV = -Vm(N)'*D*Vm(N);
end
[Evecs Evals] = signeig(VDV);
if ~isempty(find(chop(diag(Evals),...
            max(tolerance_in,eps*N*normest(VDV))) < 0))
   isisnot = 'isnot';
   if strcmp(verbose,'on'), disp('Because -VDV < 0,'), end
end
if strcmp(verbose,'on')
   if strcmp(isisnot,'isnot')
      disp('matrix input is not EDM.')
   elseif tolerance_in == eps
      disp('Matrix input is EDM to machine precision.')
   else
      disp('Matrix input is EDM to specified tolerance.')
   end
end

% find generating list
r = max(find(chop(diag(Evals),...
        max(tolerance_in,eps*N*normest(VDV))) > 0));
if isempty(r)
   r = 0;
end
if nargin < 4 | isempty(dim)
   dim = r;
else
   dim = round(dim);
end
t = r;
r = min(r,dim);
if r == 0
   X = zeros(1,N);
else
   if strcmp(use,'Vn')
      X = [zeros(r,1) diag(sqrt(diag(Evals(1:r,1:r))))*Evecs(:,1:r)'];
   else
      X = [diag(sqrt(diag(Evals(1:r,1:r))))*Evecs(:,1:r)']/sqrt(2);
   end
end
if strcmp(isisnot,'isnot') | dim < t
   Dclosest = Dx(X);
else
   Dclosest = D;
end

Subroutines for isedm()

chop()

% zeroing entries below specified absolute tolerance threshold
% -Jon Dattorro
function Y = chop(A,tolerance)

R = real(A);
I = imag(A);

if nargin == 1
   tolerance = max(size(A))*norm(A)*eps;
end
idR = find(abs(R) < tolerance);
idI = find(abs(I) < tolerance);

R(idR) = 0;
I(idI) = 0;

Y = R + i*I;

Vn()

function y = Vn(N)

y = [-ones(1,N-1);
      eye(N-1)]/sqrt(2);

Vm()

% returns EDM V matrix
function V = Vm(n)

V = [eye(n)-ones(n,n)/n];

signeig()

% Sorts signed real part of eigenvalues
% and applies sort to values and vectors.
% [Q, lam] = signeig(A)
% -Jon Dattorro

function [Q, lam] = signeig(A);

[q l] = eig(A);

lam = diag(l);
[junk id] = sort(real(lam));
id = id(length(id):-1:1);
lam = diag(lam(id));
Q = q(:,id);

if nargout < 2
   Q = diag(lam);
end

Dx()

% Make EDM from point list
function D = Dx(X)

[n,N] = size(X);
one = ones(N,1);

del = diag(X'*X);
D = del*one' + one*del' - 2*X'*X;

conic independence

The recommended subroutine lp() is a linear program solver (simplex method) from Matlab's Optimization Toolbox v2.0 (R11). Later releases of Matlab replace lp() with linprog() (interior-point method) that we find quite inferior to lp() on an assortment of problems; indeed, inherent limitation of numerical precision of solution to 1E-8 in linprog() causes failure in programs previously working with lp().

Given an arbitrary set of directions, this c.i. subroutine removes the conically dependent members. Yet a conically independent set returned is not necessarily unique. In that case, if desired, the set returned may be altered by reordering the set input.

% Test for c.i. of arbitrary directions in rows or columns of X.
% -Jon Dattorro

function [Xci, indep_str, how_many_depend] = conici(X,rowORcol,tol);

if nargin < 3
   tol=max(size(X))*eps*norm(X);
end
if nargin < 2 | strcmp(rowORcol,'col')
   rowORcol = 'col';
   Xin = X;
elseif strcmp(rowORcol,'row')
   Xin = X';
else
   disp('Invalid rowORcol input.')
   return
end
[n, N] = size(Xin);

indep_str = 'conically independent';
how_many_depend = 0;
if rank(Xin) == N
   Xci = X;
   return
end

count = 1;
new_N = N;
% remove zero rows or columns
for i=1:N
   if chop(Xin(:,count),tol)==0
      how_many_depend = how_many_depend + 1;
      indep_str = 'conically Dependent';
      Xin(:,count) = [ ];
      new_N = new_N - 1;
   else
      count = count + 1;
   end
end
% remove conic dependencies
count = 1;
newer_N = new_N;
for i=1:new_N
   if newer_N > 1
      A = [Xin(:,1:count-1) Xin(:,count+1:newer_N); -eye(newer_N-1)];
      b = [Xin(:,count); zeros(newer_N-1,1)];
      [a, lambda, how] = lp(zeros(newer_N-1,1),A,b,[ ],[ ],[ ],n,-1);
      if ~strcmp(how,'infeasible')
         how_many_depend = how_many_depend + 1;
         indep_str = 'conically Dependent';
         Xin(:,count) = [ ];
         newer_N = newer_N - 1;
      else
         count = count + 1;
      end
   end
end
if strcmp(rowORcol,'col')
   Xci = Xin;
else
   Xci = Xin';
end

lp()

LP     Linear programming.
X=LP(f,A,b) solves the linear programming problem:

             min f'x    subject to:   Ax <= b
              x

X=LP(f,A,b,VLB,VUB) defines a set of lower and upper
bounds on the design variables, X, so that the solution is
always in the range VLB <= X <= VUB.

X=LP(f,A,b,VLB,VUB,X0) sets the initial starting point to X0.

X=LP(f,A,b,VLB,VUB,X0,N) indicates that the first N constraints
defined by A and b are equality constraints.

X=LP(f,A,b,VLB,VUB,X0,N,DISPLAY) controls the level of warning
messages displayed.  Warning messages can be turned off with
DISPLAY = -1.

[X,LAMBDA]=LP(f,A,b) returns the set of Lagrangian multipliers,
LAMBDA, at the solution.

[X,LAMBDA,HOW] = LP(f,A,b) also returns a string how that
indicates error conditions at the final iteration.

LP produces warning messages when the solution is either
unbounded or infeasible.

Map of the USA

EDM, mapusa()

(\S\ref{USA}) \begin{verbatim} % Find map of USA using only distance information. % -Jon Dattorro % Reconstruction from EDM. clear all; close all;

load usalo; % From Matlab Mapping Toolbox \end{verbatim}\vspace{-3.7mm} {\tt\%\href{http://convexoptimization.com/TOOLS/USALO}{http://convexoptimization.com/TOOLS/USALO}} {\tt\%\href{http://www.mathworks.com/support/solutions/data/1-12MDM2.html?solution=1-12MDM2}{mathworks.com/support/solutions/data/1-12MDM2.html?solution=1-12MDM2}} \begin{verbatim} % To speed-up execution (decimate map data), make % 'factor' bigger positive integer. factor = 1; Mg = 2*factor; % Relative decimation factors Ms = factor; Mu = 2*factor;

gtlakelat = decimate(gtlakelat,Mg); gtlakelon = decimate(gtlakelon,Mg); statelat = decimate(statelat,Ms); statelon = decimate(statelon,Ms); uslat = decimate(uslat,Mu); uslon = decimate(uslon,Mu);

lat = [gtlakelat; statelat; uslat]*pi/180; lon = [gtlakelon; statelon; uslon]*pi/180; phi = pi/2 - lat; theta = lon; x = sin(phi).*cos(theta); y = sin(phi).*sin(theta); z = cos(phi);\end{verbatim}\pagebreak \begin{verbatim} % plot original data plot3(x,y,z), axis equal, axis off

lengthNaN = length(lat); id = find(isfinite(x)); X = [x(id)'; y(id)'; z(id)']; N = length(X(1,:))

% Make the distance matrix clear gtlakelat gtlakelon statelat statelon clear factor x y z phi theta conus clear uslat uslon Mg Ms Mu lat lon D = diag(X'*X)*ones(1,N) + ones(N,1)*diag(X'*X)' - 2*X'*X;

% destroy input data clear X

Vn = [-ones(1,N-1); speye(N-1)]; VDV = (-Vn'*D*Vn)/2;

clear D Vn pack

[evec evals flag] = eigs(VDV, speye(size(VDV)), 10, 'LR'); if flag, disp('convergence problem'), return, end; evals = real(diag(evals));

index = find(abs(evals) > eps*normest(VDV)*N); n = sum(evals(index) > 0); Xs = [zeros(n,1) diag(sqrt(evals(index)))*evec(:,index)'];

warning off; Xsplot=zeros(3,lengthNaN)*(0/0); warning on; Xsplot(:,id) = Xs; figure(2)

% plot map found via EDM. plot3(Xsplot(1,:), Xsplot(2,:), Xsplot(3,:)) axis equal, axis off \end{verbatim}

\subsubsection{USA map input-data decimation, {\tt decimate()}\index{decimate()}} \begin{verbatim} function xd = decimate(x,m) roll = 0; rock = 1; for i=1:length(x)

  if isnan(x(i))
     roll = 0;
     xd(rock) = x(i);
     rock=rock+1;
  else
     if ~mod(roll,m)
        xd(rock) = x(i);
        rock=rock+1;
     end
     roll=roll+1;
  end

end xd = xd'; \end{verbatim}

\subsection{EDM using ordinal data, {\tt omapusa()}\index{omapusa()}}\label{mouseo} (\S\ref{isomapu}) \begin{verbatim} % Find map of USA using ordinal distance information. % -Jon Dattorro clear all; close all; load usalo; % From Matlab Mapping Toolbox \end{verbatim}\vspace{-3.5mm} {\tt\%\href{http://convexoptimization.com/TOOLS/USALO}{http://convexoptimization.com/TOOLS/USALO}} {\tt\%\href{http://www.mathworks.com/support/solutions/data/1-12MDM2.html?solution=1-12MDM2}{mathworks.com/support/solutions/data/1-12MDM2.html?solution=1-12MDM2}} \begin{verbatim} factor = 1; Mg = 2*factor; % Relative decimation factors Ms = factor; Mu = 2*factor;

gtlakelat = decimate(gtlakelat,Mg); gtlakelon = decimate(gtlakelon,Mg); statelat = decimate(statelat,Ms); statelon = decimate(statelon,Ms); uslat = decimate(uslat,Mu); uslon = decimate(uslon,Mu);

lat = [gtlakelat; statelat; uslat]*pi/180; lon = [gtlakelon; statelon; uslon]*pi/180; phi = pi/2 - lat; theta = lon; x = sin(phi).*cos(theta); y = sin(phi).*sin(theta); z = cos(phi);

% plot original data plot3(x,y,z), axis equal, axis off

lengthNaN = length(lat); id = find(isfinite(x)); X = [x(id)'; y(id)'; z(id)']; N = length(X(1,:))

% Make the distance matrix clear gtlakelat gtlakelon statelat clear statelon state stateborder greatlakes clear factor x y z phi theta conus clear uslat uslon Mg Ms Mu lat lon D = diag(X'*X)*ones(1,N) + ones(N,1)*diag(X'*X)' - 2*X'*X;

% ORDINAL MDS - vectorize D count = 1; M = (N*(N-1))/2; f = zeros(M,1); for j=2:N

  for i=1:j-1
        f(count) = D(i,j);
        count = count + 1;
  end

end % sorted is f(idx) [sorted idx] = sort(f); clear D sorted X f(idx)=((1:M).^2)/M^2;

% Create ordinal data matrix O = zeros(N,N); count = 1; for j=2:N

  for i=1:j-1
        O(i,j) = f(count);
        O(j,i) = f(count);
        count = count+1;
  end

end

clear f idx

Vn = sparse([-ones(1,N-1); eye(N-1)]); VOV = (-Vn'*O*Vn)/2;

clear O Vn pack

[evec evals flag] = eigs(VOV, speye(size(VOV)), 10, 'LR'); if flag, disp('convergence problem'), return, end; evals = real(diag(evals));

Xs = [zeros(3,1) diag(sqrt(evals(1:3)))*evec(:,1:3)'];

warning off; Xsplot=zeros(3,lengthNaN)*(0/0); warning on; Xsplot(:,id) = Xs; figure(2)

% plot map found via Ordinal MDS. plot3(Xsplot(1,:), Xsplot(2,:), Xsplot(3,:)) axis equal, axis off \end{verbatim}

\section{Rank reduction subroutine, {\tt RRf()}\index{RRf()}}\label{rrf} (\S\ref{rrpro}) \begin{verbatim} % Rank Reduction function -Jon Dattorro % Inputs are: % Xstar matrix, % affine equality constraint matrix A whose rows are in svec format. % % Tolerance scheme needs revision...

function X = RRf(Xstar,A); rand('seed',23); m = size(A,1); n = size(Xstar,1); if size(Xstar,1)~=size(Xstar,2)

  disp('Rank Reduction subroutine: Xstar not square'), pause

end toler = norm(eig(Xstar))*size(Xstar,1)*1e-9; if sum(chop(eig(Xstar),toler)<0) ~= 0

  disp('Rank Reduction subroutine: Xstar not PSD'), pause

end X = Xstar; for i=1:n

  [v,d]=signeig(X);
  d(find(d<0))=0;
  rho = rank(d);
  for l=1:rho
     R(:,l,i)=sqrt(d(l,l))*v(:,l);
  end
  % find Zi
  svectRAR=zeros(m,rho*(rho+1)/2);
  cumu=0;
  for j=1:m
     temp = R(:,1:rho,i)'*svectinv(A(j,:))*R(:,1:rho,i);
     svectRAR(j,:) = svect(temp)';
     cumu = cumu + abs(temp);
  end

  % try to find sparsity pattern for Z_i
  tolerance = norm(X,'fro')*size(X,1)*1e-9;
  Ztem = zeros(rho,rho);
  pattern = find(chop(cumu,tolerance)==0);
  if isempty(pattern) % if no sparsity, do random projection
     ranp = svect(2*(rand(rho,rho)-0.5));
     Z(1:rho,1:rho,i)...
      =svectinv((eye(rho*(rho+1)/2)-pinv(svectRAR)*svectRAR)*ranp);
  else
     disp('sparsity pattern found')
     Ztem(pattern)=1;
     Z(1:rho,1:rho,i) = Ztem;
  end
  phiZ = 1;
  toler = norm(eig(Z(1:rho,1:rho,i)))*rho*1e-9;
  if sum(chop(eig(Z(1:rho,1:rho,i)),toler)<0) ~= 0
     phiZ = -1;
  end
  B(:,:,i) = -phiZ*R(:,1:rho,i)*Z(1:rho,1:rho,i)*R(:,1:rho,i)';
  % calculate t_i^*
  t(i) = max(phiZ*eig(Z(1:rho,1:rho,i)))^-1;
  tolerance = norm(X,'fro')*size(X,1)*1e-6;
  if chop(Z(1:rho,1:rho,i),tolerance)==zeros(rho,rho)
     break
  else
     X = X + t(i)*B(:,:,i);
  end

end \end{verbatim}

\pagebreak \subsection{{\tt svect()}\index{svect()}} \begin{verbatim} % Map from symmetric matrix to vector % -Jon Dattorro

function y = svect(Y,N)

if nargin == 1

  N=size(Y,1);

end

y = zeros(N*(N+1)/2,1); count = 1; for j=1:N

  for i=1:j
     if i~=j
        y(count) = sqrt(2)*Y(i,j);
     else
        y(count) = Y(i,j);
     end
     count = count + 1;
  end

end \end{verbatim}

\newpage \subsection{{\tt svectinv()}\index{svectinv()}} \begin{verbatim} % convert vector into symmetric matrix. m is dim of matrix. % -Jon Dattorro function A = svectinv(y)

m = round((sqrt(8*length(y)+1)-1)/2); if length(y) ~= m*(m+1)/2

  disp('dimension error in svectinv()');
  pause

end

A = zeros(m,m); count = 1; for j=1:m

  for i=1:m
     if i<=j
        if i==j
           A(i,i) = y(count);
        else
           A(i,j) = y(count)/sqrt(2);
           A(j,i) = A(i,j);
        end
        count = count+1;
     end
  end

end \end{verbatim}

\newpage \section{Sturm's procedure\index{procedure!Sturm}\index{procedure!dyad-decomposition}}\label{spu} This is a demonstration program that can easily be transformed to a subroutine for decomposing positive semidefinite matrix $X\,$.\, This procedure provides a nonorthogonal alternative (\S\ref{spu2}) to eigen decomposition. That particular decomposition obtained is dependent on choice of matrix $A\,$.

\begin{verbatim} % Sturm procedure to find dyad-decomposition of X -Jon Dattorro clear all

N = 4; r = 2; X = 2*(rand(r,N)-0.5); X = X'*X;

t = null(svect(X)'); A = svectinv(t(:,1));

% Suppose given matrix A is positive semidefinite %[v,d] = signeig(X); %d(1,1)=0; d(2,2)=0; d(3,3)=pi; %A = v*d*v';

tol = 1e-8; Y = X; y = zeros(size(X,1),r); rho = r; for k=1:r

   [v,d] = signeig(Y);
   v = v*sqrt(chop(d,1e-14));
   viol = 0;
   j = [ ];
   for i=2:rho
       if chop((v(:,1)'*A*v(:,1))*(v(:,i)'*A*v(:,i)),tol) ~= 0
           viol = 1;
       end
       if (v(:,1)'*A*v(:,1))*(v(:,i)'*A*v(:,i)) < 0
           j = i;
           break
       end
   end
   if ~viol
       y(:,k) = v(:,1);
   else
    if isempty(j)
     disp('Sturm procedure taking default j'), j = 2; return
    end  % debug
    alpha = (-2*(v(:,1)'*A*v(:,j)) + sqrt((2*v(:,1)'*A*v(:,j)).^2 ...
     -4*(v(:,j)'*A*v(:,j))*(v(:,1)'*A*v(:,1))))/(2*(v(:,j)'*A*v(:,j)));
    y(:,k) = (v(:,1) + alpha*v(:,j))/sqrt(1+alpha^2);
    if chop(y(:,k)'*A*y(:,k),tol) ~= 0
     alpha = (-2*(v(:,1)'*A*v(:,j)) - sqrt((2*v(:,1)'*A*v(:,j)).^2 ...
      -4*(v(:,j)'*A*v(:,j))*(v(:,1)'*A*v(:,1))))/(2*(v(:,j)'*A*v(:,j)));
     y(:,k) = (v(:,1) + alpha*v(:,j))/sqrt(1+alpha^2);
     if chop(y(:,k)'*A*y(:,k),tol) ~= 0
      disp('Zero problem in Sturm!'), return
     end  % debug
    end
   end
   Y = Y - y(:,k)*y(:,k)';
   rho = rho - 1;

end z = zeros(size(y)); e = zeros(N,N); for i=1:r

   z(:,i) = y(:,i)/norm(y(:,i));
   e(i,i) = norm(y(:,i))^2;

end lam = diag(e); [junk id] = sort(real(lam)); id = id(length(id):-1:1); z = [z(:,id(1:r)) null(z')]  % Sturm e = diag(lam(id)) [v,d] = signeig(X)  % eigenvalue decomposition X-z*e*z' traceAX = trace(A*X) \end{verbatim}

\newpage \section{Convex Iteration demonstration}\label{berg} We demonstrate implementation of a rank constraint in a semidefinite Boolean feasibility problem from \S\ref{Booleanlt}. It requires {\tt CVX}, \cite{cvx} an intuitive {\sc Matlab} interface for interior-point method solvers.

There are a finite number \,\mbox{$2^{N=_{}50}\!\approx\!{\tt 1E15}$}\, of binary \hbox{vectors \,$x\,$}.\, The feasible set of semidefinite program (\ref{ll041}) is the intersection of an elliptope with \,\mbox{$M\!=\!10$}\, halfspaces in vectorized composite $_{}G\,$.\, Size of the optimal rank$_{}$-$1$ solution set is proportional to the positive factor scaling vector \,{\tt b\,}.\, The smaller that optimal Boolean solution set, the harder this problem is to solve; indeed, it can be made as small as one point. That scale factor and initial state of random number generators, making matrix \,{\tt A}\, and \hbox{vector \,{\tt b\,}},\, are selected to demonstrate Boolean solution to one instance in a few iterations \textbf{(}a few seconds\textbf{)}, whereas sequential binary search takes one hour to test \hbox{$25.7$ million} vectors before finding one Boolean solution feasible to nonconvex \hbox{problem (\ref{obp})}. (Other parameters can be selected to realize a reversal of these timings.)

\begin{verbatim} % Discrete optimization problem demo. % -Jon Dattorro, June 4, 2007 % Find x\in{-1,1}^N such that Ax <= b clear all; format short g; M = 10; N = 50; randn('state',0); rand('state',0); A = randn(M,N); b = rand(M,1)*5;

disp('Find binary solution by convex iteration:') tic Y = zeros(N+1); count = 1; traceGY = 1e15; cvx_precision([1e-12, 1e-4]); cvx_quiet(true);\end{verbatim}\pagebreak \begin{verbatim} while 1

   cvx_begin  % requires CVX Boyd
     variable X(N,N) symmetric;
     variable x(N,1);
     G = [X,  x;
          x', 1];
     minimize(trace(G*Y));
     diag(X) == 1;
     G == semidefinite(N+1);
     A*x <= b;
   cvx_end

   [v,d,q] = svd(G);
   Y = v(:,2:N+1)*v(:,2:N+1)';
   rankG = sum(diag(d) > max(diag(d))*1e-8)
   oldtrace = traceGY;
   traceGY = trace(G*Y)
   if rankG == 1
       break
   end

   if round((oldtrace - traceGY)*1e3) == 0
       disp('STALLED');disp(' ');
       Y = -v(:,2:N+1)*(v(:,2:N+1)' + randn(N,1)*v(:,1)');
   end
   count = count + 1;

end x count toc disp('Ax <= b , x\in{-1,1}^N')\end{verbatim}\pagebreak \begin{verbatim} disp(' ');disp('Combinatorial search for a feasible binary solution:') tic for i=1:2^N

   binary = str2num(dec2bin(i-1)');
   binary(find(~binary)) = -1;
   y = [-ones(N-length(binary),1); binary];
   if sum(A*y <= b) == M
       disp('Feasible binary solution found.')
       y
       break
   end

end toc \end{verbatim}

\newpage \sectionTemplate:\sc fast max cut\label{fmck} We use the graph generator (C program) {\sc rudy} written by Giovanni Rinaldi \cite{Rinaldi} which can be found at \,\href{http://convexoptimization.com/TOOLS/RUDY}{\tt http://convexoptimization.com/TOOLS/RUDY}\, together with graph data. (\S\ref{fastmc}) \begin{verbatim} % fast max cut, Jon Dattorro, July 2007, http://convexoptimization.com clear all; format short g; tic fid = fopen('graphs12','r'); average = 0; NN = 0; s = fgets(fid); cvx_precision([1e-12, 1e-4]); cvx_quiet(true); w = 1000; while s ~= -1

   s = str2num(s);
   N = s(1);
   A = zeros(N);
   for i=1:s(2)
       s = str2num(fgets(fid));
       A(s(1),s(2)) = s(3);
       A(s(2),s(1)) = s(3);
   end
   Q = (diag(A*ones(N,1)) - A)/4;
   W = zeros(N);
   traceXW = 1e15;
   while 1
       cvx_begin                     % CVX Boyd
          variable X(N,N) symmetric;
          maximize(trace(Q*X) - w*trace(W*X));
          X == semidefinite(N);
          diag(X) == 1;
       cvx_end
       [v,d,q] = svd(X);
       W = v(:,2:N)*v(:,2:N)';
       rankX = sum(diag(d) > max(diag(d))*1e-8)
       oldtrace = traceXW;
       traceXW = trace(X*W)
       if (rankX == 1)
           break
       end
       if round((oldtrace - traceXW)*1e3) <= 0
           disp('STALLED');disp(' ')
           W = -v(:,2:N)*(v(:,2:N)' + randn(N-1,1)*v(:,1)');
       end
   end
   x = sqrt(d(1,1))*v(:,1)
   disp(' ');
   disp('Combinatorial search for optimal binary solution...')
   maxim = -1e15;
   ymax = zeros(N,1);
   for i=1:2^N
       binary = str2num(dec2bin(i-1)');
       binary(find(~binary)) = -1;
       y = [-ones(N-length(binary),1); binary];
       if y'*Q*y > maxim
           maxim = y'*Q*y;
           ymax = y;
       end
   end
   if (maxim == 0) && (abs(trace(Q*X)) <= 1e-8)
       optimality_ratio = 1
   elseif maxim <= 0
       optimality_ratio = maxim/trace(Q*X)
   else
       optimality_ratio = trace(Q*X)/maxim
   end
   ymax
   average = average + optimality_ratio;
   NN = NN + 1
   running_average = average/NN
   toc, disp(' ')
   s = fgets(fid);

end \end{verbatim}

\section{Signal dropout problem}\label{sdpr} (\S\ref{ssdp}) Requires {\tt CVX}. \cite{cvx} \begin{verbatim} %signal dropout problem -Jon Dattorro clear all close all N = 500; k = 4;  % cardinality Phi = dctmtx(N)';  % DCT basis

successes = 0; total = 0; logNormOfDifference_rate = 0; SNR_rate = 0; for i=1:1

   close all
   SNR = 0;
   na = .02;
   SNR10 = 10;
   while SNR <= SNR10
       xs = zeros(N,1);        % initialize x*
       p = randperm(N);        % calls rand
       xs(p(1:k)) = 10^(-SNR10/20) + (1-10^(-SNR10/20))*rand(k,1);
       s = Phi*xs;
       noise = na*randn(size(s));
       ll = 130;
       ul = N-ll+1;
       SNR = 20*log10(norm([s(1:ll); s(ul:N)])/norm(noise));
   end

   normof = 1e15;
   count = 1;
   y = zeros(N,1);
   cvx_quiet(true);\end{verbatim}\pagebreak

\begin{verbatim}

   while 1
       cvx_begin
           variable x(N);
           f = Phi*x;
           minimize(x'*y + norm([f(1:ll)-(s(1:ll)+noise(1:ll));
                                 f(ul:N)-(s(ul:N)+noise(ul:N))]));
           x >= 0;
       cvx_end

       if ~(strcmp(cvx_status,'Solved') ||...
            strcmp(cvx_status,'Inaccurate/Solved'))
           disp('cit failed')
           return
       end

       cvx_begin
           variable y(N);
           minimize(x'*y);
           0 <= y; y <= 1;
           y'*ones(N,1) == N-k;
       cvx_end

       if ~(strcmp(cvx_status,'Solved') ||...
            strcmp(cvx_status,'Inaccurate/Solved'))
           disp('Fan failed')
           return
       end

       cardx = sum(x > max(x)*1e-8)
       traceXW = x'*y
       oldnorm = normof;
       normof = norm([f(1:ll)-(s(1:ll)+noise(1:ll));
                      f(ul:N)-(s(ul:N)+noise(ul:N))]);
       if (cardx <= k) && (abs(oldnorm - normof) <= 1e-8)
           break
       end
       count = count + 1;
   end
   count
   figure(1);plot([s(1:ll)+noise(1:ll); noise(ll+1:ul-1);
                   s(ul:N)+noise(ul:N)]);
   hold on;  plot([s(1:ll)+noise(1:ll); zeros(ul-ll-1,1);
                   s(ul:N)+noise(ul:N)],'k')
   V = axis;
   figure(2);plot(f,'r');
    hold on; plot([s(1:ll)+noise(1:ll); noise(ll+1:ul-1);
                   s(ul:N)+noise(ul:N)]);
   axis(V);
   logNormOfDifference = 20*log10(norm(f-s)/norm(s))
   SNR
   figure(3);plot(f,'r'); hold on; plot(s);axis(V);
   figure(4);plot(f-s);axis(V);
   temp = sort([find(chop(x,max(x)*1e-8)); zeros(k-cardx,1)])...
          - sort(p(1:k))'
   if ~sum(temp)
       successes = successes + 1;
   end
   total = total + 1
   success_avg = successes/total
   logNormOfDifference_rate = logNormOfDifference_rate...
                              + logNormOfDifference;
   SNR_rate = SNR_rate + SNR;
   logNormOfDifferenceAvg = logNormOfDifference_rate/total
   SNR_avg = SNR_rate/total, disp(' ')

end successes </pre>