Matlab for Convex Optimization & Euclidean Distance Geometry

From Wikimization

(Difference between revisions)
Jump to: navigation, search

66.92.48.168 (Talk)
(New page: Made by The MathWorks [http://www.mathworks.com http://www.mathworks.com], Matlab is a high level programming language and graphical user interface for linear algebra. == isedm() == <pre>...)
Next diff →

Revision as of 17:53, 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.

Contents

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];
\end{verbatim}

\pagebreak
\subsubsection{{\tt signeig()}\index{signeig()}}
\begin{verbatim}
% 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

\end{verbatim}

\subsubsection{{\tt Dx()}\index{Dx()}}
\begin{verbatim}
% 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;
\end{verbatim}


\newpage
\section{conic independence, {\tt conici()}\index{conici()}}\label{conici}
(\S\ref{conicindy}) The recommended subroutine {\tt $_{_{}}$lp()} (\S\ref{linp}) is a linear program solver \textbf{(}\emph{simplex method}\textbf{)}
from {\sc Matlab}'s \emph{Optimization Toolbox} v2.0 (R11).\index{Matlab!lp() \emph{versus} linprog()}
Later releases of {\sc Matlab} replace {\tt $_{_{}}$lp()} with {\tt $_{_{}}$linprog()} \textbf{(}interior-point method\textbf{)}
that we find quite inferior to {\tt $_{_{}}$lp()} on an assortment of problems;
indeed, inherent limitation of numerical precision of solution to {\tt 1E-8} in\index{precision!numerical}\index{accuracy}
{\tt $_{_{}}$linprog()} causes failure in programs previously working with {\tt $_{_{}}$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.

\begin{verbatim}
% 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
\end{verbatim}

\begin{verbatim}
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

\end{verbatim}


\newpage
\subsection{{\tt lp()}\index{lp()}\index{Lagrangian}}\label{linp}
\begin{verbatim}
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.
\end{verbatim}





\newpage
\section{Map of the USA\index{map!USA}}
\subsection{EDM, {\tt mapusa()}\index{mapusa()}}\label{mouse}
(\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
\section{{\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
Personal tools