Matlab for Convex Optimization & Euclidean Distance Geometry
From Wikimization
(Difference between revisions)
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