Matlab for Convex Optimization & Euclidean Distance Geometry

From Wikimization

(Difference between revisions)
Jump to: navigation, search
(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>...)
Line 165: Line 165:
V = [eye(n)-ones(n,n)/n];
V = [eye(n)-ones(n,n)/n];
-
\end{verbatim}
+
</pre>
-
\pagebreak
+
== signeig() ==
-
\subsubsection{{\tt signeig()}\index{signeig()}}
+
<pre>
-
\begin{verbatim}
+
% Sorts signed real part of eigenvalues
% Sorts signed real part of eigenvalues
% and applies sort to values and vectors.
% and applies sort to values and vectors.
Line 188: Line 187:
Q = diag(lam);
Q = diag(lam);
end
end
 +
</pre>
-
\end{verbatim}
+
== Dx() ==
-
 
+
<pre>
-
\subsubsection{{\tt Dx()}\index{Dx()}}
+
-
\begin{verbatim}
+
% Make EDM from point list
% Make EDM from point list
function D = Dx(X)
function D = Dx(X)
Line 201: Line 199:
del = diag(X'*X);
del = diag(X'*X);
D = del*one' + one*del' - 2*X'*X;
D = del*one' + one*del' - 2*X'*X;
-
\end{verbatim}
+
</pre>
-
\newpage
+
== conic independence ==
-
\section{conic independence, {\tt conici()}\index{conici()}}\label{conici}
+
The recommended subroutine lp() is a linear program solver (''simplex method'')
-
(\S\ref{conicindy}) The recommended subroutine {\tt $_{_{}}$lp()} (\S\ref{linp}) is a linear program solver \textbf{(}\emph{simplex method}\textbf{)}
+
from Matlab's ''Optimization Toolbox'' v2.0 (R11).
-
from {\sc Matlab}'s \emph{Optimization Toolbox} v2.0 (R11).\index{Matlab!lp() \emph{versus} linprog()}
+
Later releases of Matlab replace lp() with linprog() (interior-point method)
-
Later releases of {\sc Matlab} replace {\tt $_{_{}}$lp()} with {\tt $_{_{}}$linprog()} \textbf{(}interior-point method\textbf{)}
+
that we find quite inferior to lp() on an assortment of problems;
-
that we find quite inferior to {\tt $_{_{}}$lp()} on an assortment of problems;
+
indeed, inherent limitation of numerical precision of solution to 1E-8 in
-
indeed, inherent limitation of numerical precision of solution to {\tt 1E-8} in\index{precision!numerical}\index{accuracy}
+
linprog() causes failure in programs previously working with lp().
-
{\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.
+
Given an arbitrary set of directions, this c.i. subroutine removes the conically dependent members.
Yet a conically independent set returned
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.
is not necessarily unique. In that case, if desired, the set returned may be altered by reordering the set input.
-
\begin{verbatim}
+
<pre>
% Test for c.i. of arbitrary directions in rows or columns of X.
% Test for c.i. of arbitrary directions in rows or columns of X.
% -Jon Dattorro
% -Jon Dattorro
Line 243: Line 240:
return
return
end
end
-
\end{verbatim}
 
-
\begin{verbatim}
 
count = 1;
count = 1;
new_N = N;
new_N = N;
Line 282: Line 277:
Xci = Xin';
Xci = Xin';
end
end
 +
</pre>
-
\end{verbatim}
 
-
 
+
== lp() ==
-
\newpage
+
<pre>
-
\subsection{{\tt lp()}\index{lp()}\index{Lagrangian}}\label{linp}
+
-
\begin{verbatim}
+
LP Linear programming.
LP Linear programming.
X=LP(f,A,b) solves the linear programming problem:
X=LP(f,A,b) solves the linear programming problem:
Line 316: Line 309:
LP produces warning messages when the solution is either
LP produces warning messages when the solution is either
unbounded or infeasible.
unbounded or infeasible.
-
\end{verbatim}
+
</pre>
-
 
+
-
 
+
-
 
+
-
\newpage
+
== Map of the USA ==
-
\section{Map of the USA\index{map!USA}}
+
=== EDM, mapusa() ===
-
\subsection{EDM, {\tt mapusa()}\index{mapusa()}}\label{mouse}
+
(\S\ref{USA})
(\S\ref{USA})
\begin{verbatim}
\begin{verbatim}

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

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];

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>

Personal tools