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