Matlab for Convex Optimization & Euclidean Distance Geometry

From Wikimization

(Difference between revisions)
Jump to: navigation, search
(Convex Iteration demonstration)
Current revision (13:07, 22 March 2021) (edit) (undo)
(added 1 norm distance matrix)
 
(109 intermediate revisions not shown.)
Line 1: Line 1:
-
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.
+
These MATLAB programs come from the book [http://meboo.convexoptimization.com/Meboo.html Convex Optimization & Euclidean Distance Geometry], by Dattorro, which is [http://meboo.convexoptimization.com/access.html available for download] online.
 +
 
 +
Made by The MathWorks [http://www.mathworks.com/support/books/book49103.html http://www.mathworks.com], MATLAB is a high level programming language and graphical user interface for linear algebra.
 +
 
 +
Some programs require an addon called [http://www.stanford.edu/~boyd/cvx CVX], an intuitive Matlab interface for interior-point method solvers.
== isedm() ==
== isedm() ==
<pre>
<pre>
-
% Is real D a Euclidean Distance Matrix. -Jon Dattorro
+
%Is real D a Euclidean Distance Matrix. -Jon Dattorro http://convexoptimization.com
%
%
-
% [Dclosest,X,isisnot,r] = isedm(D,tolerance,verbose,dimension,V)
+
%[Dclosest,X,isisnot,r] = isedm(D,tolerance,verbose,dimension,V)
%
%
-
% Returns: closest EDM in Schoenberg sense (default output),
+
%Returns: closest EDM in Schoenberg sense (default output),
-
% a generating list X,
+
% a generating list X,
-
% string 'is' or 'isnot' EDM,
+
% string 'is' or 'isnot' EDM,
-
% actual affine dimension r of EDM output.
+
% actual affine dimension r corresponding to EDM output.
-
% Input: matrix D,
+
%Input: candidate matrix D,
-
% optional absolute numerical tolerance for EDM determination,
+
% optional absolute numerical tolerance for EDM determination,
-
% optional verbosity 'on' or 'off',
+
% optional verbosity 'on' or 'off',
-
% optional desired affine dim of generating list X output,
+
% optional desired affine dimension of generating list X output,
-
% optional choice of 'Vn' auxiliary matrix (default) or 'V'.
+
% optional choice of 'Vn' auxiliary matrix (default) or 'V'.
function [Dclosest,X,isisnot,r] = isedm(D,tolerance_in,verbose,dim,V);
function [Dclosest,X,isisnot,r] = isedm(D,tolerance_in,verbose,dim,V);
-
isisnot = 'is';
+
isisnot = 'is';
-
N = length(D);
+
N = length(D);
if nargin < 2 | isempty(tolerance_in)
if nargin < 2 | isempty(tolerance_in)
Line 30: Line 34:
end
end
if nargin < 5 | isempty(V)
if nargin < 5 | isempty(V)
-
use = 'Vn';
+
V = 'Vn';
-
else
+
-
use = 'V';
+
end
end
-
% is empty
+
%is empty
-
if N < 1
+
if N < 1
if strcmp(verbose,'on'), disp('Input D is empty.'), end
if strcmp(verbose,'on'), disp('Input D is empty.'), end
X = [ ];
X = [ ];
Line 42: Line 44:
isisnot = 'isnot';
isisnot = 'isnot';
r = [ ];
r = [ ];
-
return
+
return
-
end
+
end
-
% is square
+
%is square
if size(D,1) ~= size(D,2)
if size(D,1) ~= size(D,2)
if strcmp(verbose,'on'), disp('An EDM must be square.'), end
if strcmp(verbose,'on'), disp('An EDM must be square.'), end
Line 53: Line 55:
return
return
end
end
-
% is real
+
%is real
if ~isreal(D)
if ~isreal(D)
if strcmp(verbose,'on'), disp('Because an EDM is real,'), end
if strcmp(verbose,'on'), disp('Because an EDM is real,'), end
Line 59: Line 61:
D = real(D);
D = real(D);
end
end
-
 
+
%is nonnegative
-
% is nonnegative
+
if sum(sum(chop(D,tolerance) < 0))
if sum(sum(chop(D,tolerance) < 0))
isisnot = 'isnot';
isisnot = 'isnot';
if strcmp(verbose,'on'), disp('Because an EDM is nonnegative,'),end
if strcmp(verbose,'on'), disp('Because an EDM is nonnegative,'),end
end
end
-
% is symmetric
+
%is symmetric
if sum(sum(abs(chop((D - D')/2,tolerance)) > 0))
if sum(sum(abs(chop((D - D')/2,tolerance)) > 0))
isisnot = 'isnot';
isisnot = 'isnot';
if strcmp(verbose,'on'), disp('Because an EDM is symmetric,'), end
if strcmp(verbose,'on'), disp('Because an EDM is symmetric,'), end
-
D = (D + D')/2; % only required condition
+
D = (D + D')/2; %only required condition
end
end
-
% has zero diagonal
+
%has zero diagonal
if sum(abs(diag(chop(D,tolerance))) > 0)
if sum(abs(diag(chop(D,tolerance))) > 0)
isisnot = 'isnot';
isisnot = 'isnot';
-
if strcmp(verbose,'on')
+
if strcmp(verbose,'on')
disp('Because an EDM has zero main diagonal,')
disp('Because an EDM has zero main diagonal,')
end
end
end
end
-
% is EDM
+
%is EDM
-
if strcmp(use,'Vn')
+
if strcmp(V,'Vn')
VDV = -Vn(N)'*D*Vn(N);
VDV = -Vn(N)'*D*Vn(N);
else
else
Line 90: Line 91:
end
end
if strcmp(verbose,'on')
if strcmp(verbose,'on')
-
if strcmp(isisnot,'isnot')
+
if strcmp(isisnot,'isnot')
disp('matrix input is not EDM.')
disp('matrix input is not EDM.')
elseif tolerance_in == eps
elseif tolerance_in == eps
Line 99: Line 100:
end
end
-
% find generating list
+
%find generating list
r = max(find(chop(diag(Evals),max(tolerance_in,eps*N*normest(VDV))) > 0));
r = max(find(chop(diag(Evals),max(tolerance_in,eps*N*normest(VDV))) > 0));
if isempty(r)
if isempty(r)
Line 111: Line 112:
t = r;
t = r;
r = min(r,dim);
r = min(r,dim);
-
if r == 0
+
if r == 0
X = zeros(1,N);
X = zeros(1,N);
else
else
-
if strcmp(use,'Vn')
+
if strcmp(V,'Vn')
X = [zeros(r,1) diag(sqrt(diag(Evals(1:r,1:r))))*Evecs(:,1:r)'];
X = [zeros(r,1) diag(sqrt(diag(Evals(1:r,1:r))))*Evecs(:,1:r)'];
else
else
Line 120: Line 121:
end
end
end
end
-
if strcmp(isisnot,'isnot') | dim < t
+
if strcmp(isisnot,'isnot') | dim < t
Dclosest = Dx(X);
Dclosest = Dx(X);
-
else
+
else
-
Dclosest = D;
+
Dclosest = D;
end
end
</pre>
</pre>
 +
 +
[http://www.convexoptimization.com/TOOLS/isedm.pdf Alencar, Bonates, Lavor, & Liberti]
 +
propose a more robust isedm() by replacing eigenvalue decomposition with singular value decomposition.
=== Subroutines for isedm() ===
=== Subroutines for isedm() ===
Line 195: Line 199:
one = ones(N,1);
one = ones(N,1);
-
del = diag(X'*X);
+
XTX = X'*X;
-
D = del*one' + one*del' - 2*X'*X;
+
del = diag(XTX);
 +
D = del*one' + one*del' - 2*XTX;
</pre>
</pre>
-
== conic independence ==
+
== taxicab (1 norm) distance matrix ==
-
The recommended subroutine lp() is a linear program solver (''simplex method'')
+
==== D1() ====
-
from Matlab's ''Optimization Toolbox'' v2.0 (R11).
+
<pre>
 +
%Distance matrix from point list in 1 norm
 +
function D = D1(X);
-
Later releases of Matlab replace lp() with linprog() (interior-point method)
+
[n,N] = size(X);
-
that we find quite inferior to lp() on an assortment of problems;
+
-
indeed, inherent limitation of numerical precision of solution to 1E-8 in
+
D = kron(eye(N),ones(1,n))*abs(kron(ones(1,N),X(:)) - kron(ones(N,1),X));
-
linprog() causes failure in programs previously working with lp().
+
</pre>
 +
 +
== conic independence ==
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.
Line 218: Line 226:
=== conici() ===
=== conici() ===
<pre>
<pre>
-
% Test for c.i. of arbitrary directions in rows or columns of X.
+
% Remove c.i. directions in rows or columns of X. -Jon Dattorro
-
% -Jon Dattorro
+
% function [Xci, indep_str, how_many_depend, kept] = conici(X,'rowORcol',tolerance);
-
function [Xci, indep_str, how_many_depend] = conici(X,rowORcol,tol);
+
function [Xci, indep_str, how_many_depend, kept] = conici(X,rowORcol,tol);
if nargin < 3
if nargin < 3
-
tol=max(size(X))*eps*norm(X);
+
tol=max(size(X))*eps*normest(X);
end
end
-
if nargin < 2 | strcmp(rowORcol,'col')
+
if nargin < 2 || strcmp(rowORcol,'col') || isempty(rowORcol)
rowORcol = 'col';
rowORcol = 'col';
Xin = X;
Xin = X;
elseif strcmp(rowORcol,'row')
elseif strcmp(rowORcol,'row')
Xin = X';
Xin = X';
-
else
+
else
disp('Invalid rowORcol input.')
disp('Invalid rowORcol input.')
return
return
Line 237: Line 245:
[n, N] = size(Xin);
[n, N] = size(Xin);
-
indep_str = 'conically independent';
+
indep_str = 'conically independent';
how_many_depend = 0;
how_many_depend = 0;
-
if rank(Xin) == N
+
kept = [1:N]';
-
Xci = X;
+
if rank(Xin,tol) == N
 +
Xci = chop(X,tol);
return
return
end
end
Line 246: Line 255:
count = 1;
count = 1;
new_N = N;
new_N = N;
-
% remove zero rows or columns
+
%remove zero columns
-
for i=1:N
+
for i=1:new_N
if chop(Xin(:,count),tol)==0
if chop(Xin(:,count),tol)==0
how_many_depend = how_many_depend + 1;
how_many_depend = how_many_depend + 1;
indep_str = 'conically Dependent';
indep_str = 'conically Dependent';
Xin(:,count) = [ ];
Xin(:,count) = [ ];
 +
kept(count) = [ ];
new_N = new_N - 1;
new_N = new_N - 1;
else
else
Line 257: Line 267:
end
end
end
end
-
% remove conic dependencies
+
%remove identical columns
 +
D = sqrt(Dx(Xin));
 +
T = tril(ones(new_N,new_N))*normest(D);
 +
ir = find(D+T < tol);
 +
if ~isempty(ir)
 +
[iir jir] = ind2sub(size(D),ir);
 +
indep_str = 'conically Dependent';
 +
sizebefore = size(Xin);
 +
Xin(:,jir) = [ ];
 +
sizeafter = size(Xin);
 +
kept(jir) = [ ];
 +
new_N = size(Xin,2);
 +
how_many_depend = how_many_depend + sizebefore(2)-sizeafter(2);
 +
end
 +
%remove conic dependencies
count = 1;
count = 1;
newer_N = new_N;
newer_N = new_N;
-
for i=1:new_N
+
for i=1:newer_N
if newer_N > 1
if newer_N > 1
A = [Xin(:,1:count-1) Xin(:,count+1:newer_N); -eye(newer_N-1)];
A = [Xin(:,1:count-1) Xin(:,count+1:newer_N); -eye(newer_N-1)];
Line 269: Line 293:
indep_str = 'conically Dependent';
indep_str = 'conically Dependent';
Xin(:,count) = [ ];
Xin(:,count) = [ ];
 +
kept(count) = [ ];
newer_N = newer_N - 1;
newer_N = newer_N - 1;
else
else
Line 276: Line 301:
end
end
if strcmp(rowORcol,'col')
if strcmp(rowORcol,'col')
-
Xci = Xin;
+
Xci = chop(Xin, tol);
else
else
-
Xci = Xin';
+
Xci = chop(Xin',tol);
end
end
</pre>
</pre>
- 
=== lp() ===
=== lp() ===
 +
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().
 +
 +
The source code is available here: [http://convexoptimization.com/TOOLS/lp.m lp.m] which calls [http://convexoptimization.com/TOOLS/myqpsubold.m myqpsubold.m]
<pre>
<pre>
LP Linear programming.
LP Linear programming.
Line 313: Line 347:
unbounded or infeasible.
unbounded or infeasible.
</pre>
</pre>
 +
== Map of the USA ==
== Map of the USA ==
Line 390: Line 425:
</pre>
</pre>
-
=== USA map input-data decimation, decimate() ===
+
==== USA map input-data decimation, decimate() ====
<pre>
<pre>
function xd = decimate(x,m)
function xd = decimate(x,m)
Line 410: Line 445:
xd = xd';
xd = xd';
</pre>
</pre>
- 
- 
=== EDM using ordinal data, omapusa() ===
=== EDM using ordinal data, omapusa() ===
Line 492: Line 525:
clear O Vn
clear O Vn
-
pack
 
[evec evals flag] = eigs(VOV, speye(size(VOV)), 10, 'LR');
[evec evals flag] = eigs(VOV, speye(size(VOV)), 10, 'LR');
Line 508: Line 540:
axis equal, axis off
axis equal, axis off
</pre>
</pre>
 +
== Rank reduction subroutine, RRf() ==
== Rank reduction subroutine, RRf() ==
Line 630: Line 663:
-
== Sturm's procedure ==
+
== Sturm & Zhang's procedure for constructing dyad-decomposition ==
This is a demonstration program that can easily be transformed
This is a demonstration program that can easily be transformed
to a subroutine for decomposing positive semidefinite matrix <math>\,X\,</math>.
to a subroutine for decomposing positive semidefinite matrix <math>\,X\,</math>.
Line 709: Line 742:
-
== Convex Iteration demonstration ==
+
== ''Convex Iteration'' demonstration - Boolean feasibility ==
We demonstrate implementation of a rank constraint in a semidefinite Boolean feasibility problem.
We demonstrate implementation of a rank constraint in a semidefinite Boolean feasibility problem.
-
It requires CVX, an intuitive Matlab interface for interior-point method solvers.
+
It requires [http://www.stanford.edu/~boyd/cvx CVX], an intuitive Matlab interface for interior-point method solvers.
-
There are a finite number <math>\,2^{N=_{}50}\!\approx\!1E15\,</math> of binary vectors <math>\,x\,</math>.
+
There are a finite number <math>\,2^{N=_{}50}\!\approx 1E15\,</math> of binary vectors <math>\,x\,</math>.
The feasible set of the semidefinite program from the book is the intersection of an elliptope
The feasible set of the semidefinite program from the book is the intersection of an elliptope
-
with <math>\,M\!=\!10\,</math> halfspaces in vectorized composite <math>G</math>.
+
with <math>\,M\!=10\,</math> halfspaces in vectorized composite <math>\,G\,</math>.
Size of the optimal rank-1 solution set
Size of the optimal rank-1 solution set
-
is proportional to the positive factor scaling vector b.
+
is proportional to the positive factor scaling vector '''b'''.
The smaller that optimal Boolean solution set, the harder this problem is to solve;
The smaller that optimal Boolean solution set, the harder this problem is to solve;
indeed, it can be made as small as one point.
indeed, it can be made as small as one point.
-
That scale factor and initial state of random number generators, making matrix A and vector b,
+
That scale factor and initial state of random number generators, making matrix '''A''' and vector '''b''',
are selected to demonstrate Boolean solution to one instance in a few iterations (a few seconds);
are selected to demonstrate Boolean solution to one instance in a few iterations (a few seconds);
whereas sequential binary search takes one hour to test 25.7 million vectors before finding
whereas sequential binary search takes one hour to test 25.7 million vectors before finding
-
one Boolean solution feasible to nonconvex problem from the book.
+
one Boolean solution feasible to the nonconvex problem from the book.
(Other parameters can be selected to realize a reversal of these timings.)
(Other parameters can be selected to realize a reversal of these timings.)
Line 754: Line 787:
while 1
while 1
-
cvx_begin % requires CVX Boyd
+
cvx_begin
variable X(N,N) symmetric;
variable X(N,N) symmetric;
variable x(N,1);
variable x(N,1);
Line 798: Line 831:
end
end
toc
toc
 +
</pre>
 +
 +
== ''Convex Iteration'' rank-1 ==
 +
This program demonstrates how a semidefinite problem with a rank-''r'' constraint is equivalently transformed into
 +
a problem sequence having a rank-1 constraint; discussed at the end of [http://meboo.convexoptimization.com/Meboo.html Chapter 4]. Requires [http://www.stanford.edu/~boyd/cvx CVX].
<pre>
<pre>
 +
%Convex Iteration Rank-1, -Jon Dattorro October 2007
 +
% find X
 +
% subject to A vec X = b
 +
% X positive semidefinite
 +
% rank X <= r
 +
% X\in\symm^N
 +
clear all
 +
N=26;
 +
M=10;
 +
randn('seed',100);
 +
A = randn(M,N*N);
 +
b = randn(M,1);
 +
 +
r = 2;
 +
W = ones(r*N);
 +
Q11 = zeros(N,N);
 +
Q22 = zeros(N,N);
 +
count = 1;
 +
traceGW = 1e15;
 +
normof = 1e15;
 +
w1 = 5;
 +
cvx_quiet(true);
 +
cvx_precision([1e-10 1e-4]);
 +
while 1
 +
cvx_begin
 +
variable L(r,1);
 +
X = L(1)*Q11 + L(2)*Q22;
 +
minimize(norm(A*X(:) - b));
 +
L >= 0;
 +
cvx_end
 +
 +
% find Q
 +
cvx_begin
 +
variable Q11(N,N) symmetric;
 +
variable Q12(N,N);
 +
variable Q22(N,N) symmetric;
 +
 +
trace(Q11) == 1;
 +
trace(Q22) == 1;
 +
trace(Q12) == 0;
 +
 +
G = [Q11, Q12;
 +
Q12', Q22];
 +
G == semidefinite(r*N);
 +
 +
X = L(1)*Q11 + L(2)*Q22;
 +
 +
minimize(w1*trace(G*W) + norm(A*X(:) - b));
 +
cvx_end
 +
 +
[v,d,q] = svd(full(G));
 +
W = v(:,2:r*N)*v(:,2:r*N)';
 +
rankG = sum(diag(d) > max(diag(d))*1e-6)
 +
oldtraceGW = traceGW;
 +
traceGW = trace(G*W)
 +
oldnorm = normof;
 +
normof = norm(A*X(:) - b)
 +
if (rankG == 1) && (norm(A*X(:) - b) < 1e-6)
 +
break
 +
end
 +
 +
if round((oldtraceGW - traceGW)*1e3) <= 0 && round((oldnorm - normof)*1e3) <= 0
 +
disp('STALLED');disp(' ')
 +
W = -v(:,2:r*N)*(v(:,2:r*N)' + randn(r*N-1,1)*v(:,1)');
 +
end
 +
count = count + 1;
 +
end
 +
count
 +
</pre>
 +
 +
==''Convex Iteration'' rank-1, 2013==
 +
<pre>
 +
%prototypical rank constrained sdp by rank-1 transformation, -Jon Dattorro October 2013
 +
% find X
 +
% subject to A vec X = b
 +
% X positive semidefinite
 +
% rank X <= r
 +
 +
% X is n x n symmetric
 +
clc;
 +
tachometer = 0;
 +
accelcount = 0;
 +
backoutcount = 0;
 +
iavg = 0;
 +
for exemplar = 1:1
 +
save exemplar exemplar iavg tachometer accelcount backoutcount
 +
clear all;
 +
load exemplar
 +
rand('twister',sum(100*clock));
 +
randn('state',sum(100*clock));
 +
m = 25; %given matrix A is m x n(n+1)/2
 +
n = 7;
 +
r = 3; %assumed rank of matrix
 +
starting_window_length = 25;
 +
quant = 1e6;
 +
So = rand(r,1);
 +
Uo = orth(randn(n,r));
 +
Xo = Uo*diag(So)*Uo';
 +
A = randn(m,n*(n+1)/2);
 +
b = A*svect2(Xo);
 +
 +
Y = eye(n*r);
 +
cvx_quiet('true');
 +
tic
 +
accelerant = 1;
 +
iavgadj = 0;
 +
Z = [];
 +
backout = false;
 +
i = 0;
 +
it = 0;
 +
while 1
 +
i = i + 1;
 +
it = it + 1;
 +
if i > 1 && isempty(strfind(cvx_status,'Solved'))
 +
temp = cvx_solver;
 +
if ~isempty(strfind(temp,'SeDuMi'))
 +
cvx_solver('SDPT3');
 +
else
 +
cvx_solver('SeDuMi');
 +
end
 +
iavgadj = iavgadj + 1;
 +
temp = cvx_solver;
 +
disp(sprintf('switching solver to %s\n',temp));
 +
else
 +
cvx_solver('SeDuMi');
 +
end
 +
Zold = Z;
 +
U = sparse(n,r);
 +
T = [];
 +
k = 1;
 +
for j=1:r
 +
idx(k) = 1 + length(T); % index Z matrix by block
 +
T = [T; U(:,j)];
 +
k = k + 1;
 +
end
 +
idx(k) = 1 + length(T);
 +
cvx_begin
 +
variable Z(n*r,n*r) symmetric;
 +
Z == semidefinite(n*r);
 +
XX = Z(idx(1):idx(2)-1,idx(1):idx(2)-1); % sum U_ii = X
 +
for j=2:r
 +
XX = XX + Z(idx(j):idx(j+1)-1,idx(j):idx(j+1)-1);
 +
end
 +
A*svect2(XX) == b;
 +
for j=1:r-1
 +
for l=j+1:r
 +
trace(Z(idx(j):idx(j+1)-1,idx(l):idx(l+1)-1)) == 0; % u_i u_j orthogonality
 +
end
 +
end
 +
minimize(trace(Y(:,:,it)'*Z));
 +
cvx_end
 +
if it == 1, clc; end
 +
disp(sprintf('exemplar = %d',exemplar))
 +
disp(sprintf('cvx_status = %s',cvx_status))
 +
if exemplar > 1
 +
disp(sprintf('average iterations = %d',round(iavg/(exemplar-1))));
 +
end
 +
[u s v] = svd(Z);
 +
if it > 1 && (isempty(strfind(cvx_status,'Solved')) || (accelerant > 1 && sum(s(ambig+1:end)) > darkmatter(it-1))) && ~tack
 +
it = it-1;
 +
backout = true;
 +
Z = Zold;
 +
[u s v] = svd(Z);
 +
end
 +
tack = false;
 +
temp = diag(s);
 +
coordinates = chop(temp(1:min(r*2,size(s,1))),quant^-1)
 +
ambig = max(1,length(find(abs(coordinates-coordinates(1)) < quant^-1)));
 +
Y(:,:,it+1) = u(:,ambig+1:end)*diag(sign(sum(u(:,ambig+1:end).*v(:,ambig+1:end))))*v(:,ambig+1:end)';
 +
Y(:,:,it+1) = (Y(:,:,it+1) + Y(:,:,it+1)')/2;
 +
darkmatter(it) = trace(Y(:,:,it+1)'*Z);
 +
disp(sprintf('traceYZ = %g',darkmatter(it)))
 +
if it > 1
 +
if it == 2, close all; end
 +
figure(1)
 +
if iavg, magi = floor(iavg/(exemplar-1)); else magi = starting_window_length; end
 +
if length(darkmatter) > magi
 +
plot(it-magi:it,darkmatter(end-magi:end));
 +
else
 +
plot(1:length(darkmatter),darkmatter);
 +
end
 +
set(gcf,'position',[400 430 256 256])
 +
pause(0.07);
 +
end
 +
if it > 2 % make a line specified by two points
 +
X(:,1) = svect(Y(:,:,it-1) + Y(:,:,it))/2;
 +
X(:,2) = svect(Y(:,:,it) + Y(:,:,it+1))/2;
 +
XVn = X*Vn(2);
 +
Px1 = X(:,1) + XVn*(XVn'*(svect(Y(:,:,it-1))-X(:,1)))/norm(XVn)^2;
 +
Px2 = X(:,1) + XVn*(XVn'*(svect(Y(:,:,it)) - X(:,1)))/norm(XVn)^2;
 +
Px3 = X(:,1) + XVn*(XVn'*(svect(Y(:,:,it+1))-X(:,1)))/norm(XVn)^2;
 +
straight = (norm(svect(Y(:,:,it-1))-Px1,1) + norm(svect(Y(:,:,it))-Px2,1) + norm(svect(Y(:,:,it+1))-Px3,1))/accelerant^2;
 +
disp(sprintf('straight = %g',straight));
 +
separation = norm(X(:,2) - X(:,1),1);
 +
disp(sprintf('separation %g',separation));
 +
if ~backout
 +
if straight < 1 && separation < 1
 +
accelerant = 2^-log10(straight);
 +
else
 +
accelerant = accelerant/2;
 +
end
 +
else
 +
disp('backing out')
 +
accelerant = accelerant/2;
 +
backoutcount = backoutcount + 1;
 +
end
 +
if accelerant > 1
 +
disp(sprintf('accelerant %g',accelerant))
 +
accelcount = accelcount + 1;
 +
elseif ~backout
 +
disp(' ')
 +
end
 +
if accelerant < 1, accelerant = 1; end
 +
Y(:,:,it+1) = Y(:,:,it+1) + svectinv((accelerant-1)*(X(:,2) - X(:,1)));
 +
end
 +
disp(' ')
 +
if length(find(coordinates)) <= 1
 +
[QQ DD VV] = svd(XX);
 +
eigens = diag(DD).*sign(sum(QQ.*VV)');
 +
if (sum(abs(A*svect(XX)-b)) <= quant^-1) && ~length(find(chop(eigens,quant^-1) < 0)) && (length(find(chop(eigens,quant^-1))) <= r)
 +
disp(sprintf('iterations = %d',i))
 +
iavg = max(0, iavg + i - iavgadj);
 +
toc
 +
break
 +
end
 +
else
 +
if iavg && ~mod(i,max(round(log10(quant)*iavg/(exemplar-1)),starting_window_length)) && separation < 1
 +
temp = randn(n*r,n*r);
 +
Y(:,:,it+1) = temp*temp';
 +
iavgadj = iavgadj + max(round(log10(quant)*iavg/(exemplar-1)),starting_window_length);
 +
tachometer = tachometer + 1;
 +
tack = true;
 +
elseif ~iavg && ~mod(i,round(log10(quant)*starting_window_length)) && separation < 1
 +
temp = randn(n*r,n*r);
 +
Y(:,:,it+1) = temp*temp';
 +
iavgadj = iavgadj + round(log10(quant)*starting_window_length);
 +
tachometer = tachometer + 1;
 +
tack = true;
 +
elseif ambig > 1 && ~sum(abs(coordinates(ambig+1:end)))
 +
temp = randn(n*r,n*r);
 +
Y(:,:,it+1) = temp*temp';
 +
iavgadj = iavgadj + 1;
 +
tachometer = tachometer + 1;
 +
tack = true;
 +
end
 +
end
 +
backout = false;
 +
end
 +
disp(sprintf('average iterations = %d',round(iavg/exemplar)))
 +
disp(sprintf('residual = %g',sum(abs(A*svect(XX)-b))))
 +
disp(sprintf('rankX = %g',length(find(chop(eigens,quant^-1)))))
 +
disp(sprintf('tack = %d percent',round(100*tachometer/exemplar)))
 +
if accelcount, disp(sprintf('backout/accelerant ratio = %g',backoutcount/accelcount)), end
 +
end
 +
</pre>
 +
 +
 +
=== svect2() ===
 +
<pre>
 +
%Map from symmetric matrix to vector
 +
% -Jon Dattorro
 +
 +
function y = svect2(Y)
 +
N = size(Y,1);
 +
ind = zeros(N*(N+1)/2,1);
 +
beta = zeros(N*(N+1)/2,1);
 +
count = 1;
 +
for j=1:N
 +
for i=1:j
 +
ind(count) = sub2ind(size(Y),i,j);
 +
if i==j
 +
beta(count) = 1;
 +
else
 +
beta(count) = sqrt(2);
 +
end
 +
count = count + 1;
 +
end
 +
end
 +
y = Y(ind).*beta;
 +
</pre>
 +
 +
===dvect3()===
 +
<pre>
 +
%Map from EDM to vector
 +
function y = dvect3(Y)
 +
N = size(Y,1);
 +
ind = zeros(N*(N-1)/2,1);
 +
beta = zeros(N*(N-1)/2,1);
 +
count = 1;
 +
for j=2:N
 +
for i=1:j-1
 +
ind(count) = sub2ind(size(Y),i,j);
 +
if i~=j
 +
beta(count) = 1;
 +
end
 +
count = count + 1;
 +
end
 +
end
 +
y = Y(ind).*beta;
 +
</pre>
 +
 +
=== Singular Value Decomposition (SVD) by rank-1 transformation ===
 +
Introduced at the end of [http://meboo.convexoptimization.com/Meboo.html Chapter 4] in 2014 book version.
 +
This Matlab program requires [http://www.cvxr.com CVX].
 +
<pre>
 +
%svd decomposition by rank-1 transformation -Jon Dattorro October 2013
 +
% X = U diag(S) V'
 +
clc;
 +
iavg = 0;
 +
for exemplar = 1:1
 +
save exemplar exemplar iavg
 +
clear all;
 +
load exemplar
 +
rand('twister',sum(100*clock));
 +
randn('state',sum(100*clock));
 +
m = 6; %given matrix A is m x n, rank r
 +
n = 7;
 +
r = 3; %assumed rank of matrix
 +
 +
quant = 1e6;
 +
Uo = orth(randn(m,r));
 +
So = rand(r,1);
 +
Vo = orth(randn(n,r));
 +
A = Uo*diag(So)*Vo'; % A = U diag(S) V'. Assume H = U diag(S)
 +
 +
Y = eye(2*m*r + n*r + r + 1);
 +
cvx_quiet('true');
 +
cvx_solver('sedumi');
 +
tic
 +
accelerant = 1;
 +
Z = [];
 +
backout = false;
 +
i = 0;
 +
it = 0;
 +
while 1
 +
i = i + 1;
 +
it = it + 1;
 +
Zold = Z;
 +
cvx_begin
 +
variable H(m,r);
 +
variables U(m,r) S(r,1) V(n,r);
 +
T = [];
 +
k = 1;
 +
idx(k) = 1; % index Z matrix by block
 +
k = k + 1;
 +
for j=1:r
 +
idx(k) = 2 + length(T);
 +
T = [T; H(:,j)];
 +
k = k + 1;
 +
end
 +
for j=1:r
 +
idx(k) = 2 + length(T);
 +
T = [T; U(:,j)];
 +
k = k + 1;
 +
end
 +
idx(k) = 2 + length(T);
 +
T = [T; S];
 +
k = k + 1;
 +
for j=1:r
 +
idx(k) = 2 + length(T);
 +
T = [T; V(:,j)];
 +
k = k + 1;
 +
end
 +
idx(k) = 2 + length(T);
 +
variable ZZT(length(T),length(T)) symmetric;
 +
Z = [1 T';
 +
T ZZT];
 +
for j=1:r
 +
% variables J58(n,n) J69(n,n) J710(n,n) symmetric; % H U' symmetry
 +
dvect3(Z(idx(j+1):idx(j+2)-1,idx(r+j+1):idx(r+j+2)-1)) == dvect3(Z(idx(j+1):idx(j+2)-1,idx(r+j+1):idx(r+j+2)-1)');
 +
end
 +
Z == semidefinite(2*m*r + n*r + r + 1);
 +
S >= 0;
 +
% J512 + J613 + J714 == A; % H V' = A
 +
AA = Z(idx(2):idx(3)-1,idx(2*r+3):idx(2*r+4)-1);
 +
for j=2:r
 +
AA = AA + Z(idx(j+1):idx(j+2)-1,idx(2*r+j+2):idx(2*r+j+3)-1);
 +
end
 +
AA == A;
 +
% Z=[1 H(:,1)' H(:,2)' H(:,3)' U(:,1)' U(:,2)' U(:,3)' S' V(:,1)' V(:,2)' V(:,3)'
 +
% H(:,1) J55 J56 J57 J58 J59 J510 J511 J512 J513 J514
 +
% H(:,2) J56' J66 J67 J68 J69 J610 J611 J612 J613 J614
 +
% H(:,3) J57' J67' J77 J78 J79 J710 J711 J712 J713 J714
 +
% U(:,1) J58' J68' J78' J88 J89 J810 J811 J812 J813 J814
 +
% U(:,2) J59' J69' J79' J89' J99 J910 J911 J912 J913 J914
 +
% U(:,3) J510' J610' J710' J810' J910' J1010 J1011 J1012 J1013 J1014
 +
% S J511' J611' J711' J811' J911' J1011' J1111 J1112 J1113 J1114
 +
% V(:,1) J512' J612' J712' J812' J912' J1012' J1112' J1212 J1213 J1214
 +
% V(:,2) J513' J613' J713' J813' J913' J1013' J1113' J1213' J1313 J1314
 +
% V(:,3) J514' J614' J714' J814' J914' J1014' J1114' J1214' J1314' J1414];
 +
% H == [J811(:,1) J911(:,2) J1011(:,3)]; % H = U diag(S)
 +
HH = [];
 +
for j=1:r
 +
HH = [HH Z(idx(r+j+1):idx(r+j+2)-1,idx(2*r+2)+j-1)];
 +
end
 +
HH == H;
 +
% trace(J58) == S(1,1); trace(J69) == S(2,1); trace(J710) == S(3,1); % H U' singular values
 +
for j=1:r
 +
trace(Z(idx(j+1):idx(j+2)-1,idx(r+j+1):idx(r+j+2)-1)) == S(j,1);
 +
end
 +
% trace(J56) == 0; trace(J57) == 0; trace(J67) == 0; % H orthogonality
 +
for j=2:r
 +
for l=j:r
 +
trace(Z(idx(j):idx(j+1)-1,idx(l+1):idx(l+2)-1)) == 0;
 +
end
 +
end
 +
% trace(J59) == 0; trace(J510) == 0; % U' H perpendicularity
 +
% trace(J68) == 0; trace(J610) == 0;
 +
% trace(J78) == 0; trace(J79) == 0;
 +
for j=1:r
 +
for l=1:r
 +
if l~=j
 +
trace(Z(idx(j+1):idx(j+2)-1,idx(r+l+1):idx(r+l+2)-1)) == 0;
 +
end
 +
end
 +
end
 +
% trace(J88) == 1; trace(J99) == 1; trace(J1010) == 1; %column normalization
 +
% trace(J1212) == 1; trace(J1313) == 1; trace(J1414) == 1;
 +
for j=1:r
 +
trace(Z(idx(r+j+1):idx(r+j+2)-1,idx(r+j+1):idx(r+j+2)-1)) == 1;
 +
trace(Z(idx(2*r+j+2):idx(2*r+j+3)-1,idx(2*r+j+2):idx(2*r+j+3)-1)) == 1;
 +
end
 +
% trace(J89) == 0; trace(J810) == 0; trace(J910) == 0; %orthogonality
 +
% trace(J1213) == 0; trace(J1214) == 0; trace(J1314) == 0;
 +
for j=2:r
 +
for l=j:r
 +
trace(Z(idx(r+j):idx(r+j+1)-1,idx(r+l+1):idx(r+l+2)-1)) == 0;
 +
trace(Z(idx(2*r+j+1):idx(2*r+j+2)-1,idx(2*r+l+2):idx(2*r+l+3)-1)) == 0;
 +
end
 +
end
 +
% trace(J55) == J1111(1,1); trace(J66) == J1111(2,2); trace(J77) == J1111(3,3); % H inner product
 +
for j=1:r
 +
trace(Z(idx(j+1):idx(j+2)-1,idx(j+1):idx(j+2)-1)) == Z(idx(2*r+2)+j-1,idx(2*r+2)+j-1);
 +
end
 +
minimize(trace(Y(:,:,it)'*Z));
 +
cvx_end
 +
if it == 1, clc; end
 +
disp(sprintf('exemplar = %d',exemplar))
 +
disp(sprintf('cvx_status = %s',cvx_status))
 +
if exemplar > 1
 +
disp(sprintf('average iterations = %d',round(iavg/(exemplar-1))));
 +
end
 +
[u s] = signeig(Z);
 +
if it > 1 && (isempty(strfind(cvx_status,'Solved')) || (accelerant > 1 && sum(s(ambig+1:end)) > darkmatter(it-1)))
 +
it = it-1;
 +
backout = true;
 +
Z = Zold;
 +
[u s] = signeig(Z);
 +
end
 +
temp = diag(s);
 +
coordinates = chop(temp(1:min(r*2,size(s,1))),quant^-1)
 +
ambig = length(find(abs(coordinates-coordinates(1)) < quant^-1));
 +
Y(:,:,it+1) = u(:,ambig+1:end)*u(:,ambig+1:end)';
 +
darkmatter(it) = trace(Y(:,:,it+1)'*Z);
 +
disp(sprintf('traceYZ = %g',darkmatter(it)))
 +
if it > 1
 +
if it == 2, close all; end
 +
figure(1)
 +
if iavg, magi = floor(iavg/(exemplar-1)); else magi = 22; end %22 is average number iterations for m=n=7 r=3
 +
if length(darkmatter) > magi
 +
plot(it-magi:it,darkmatter(end-magi:end));
 +
else
 +
plot(1:length(darkmatter),darkmatter);
 +
end
 +
set(gcf,'position',[400 430 256 256])
 +
pause(0.07);
 +
end
 +
if it > 2 % make a line specified by two points
 +
X(:,1) = svect(Y(:,:,it-1) + Y(:,:,it))/2;
 +
X(:,2) = svect(Y(:,:,it) + Y(:,:,it+1))/2;
 +
XVn = X*Vn(2);
 +
Px1 = X(:,1) + XVn*(XVn'*(svect(Y(:,:,it-1))-X(:,1)))/norm(XVn)^2;
 +
Px2 = X(:,1) + XVn*(XVn'*(svect(Y(:,:,it)) - X(:,1)))/norm(XVn)^2;
 +
Px3 = X(:,1) + XVn*(XVn'*(svect(Y(:,:,it+1))-X(:,1)))/norm(XVn)^2;
 +
straight = (norm(svect(Y(:,:,it-1))-Px1,1) + norm(svect(Y(:,:,it))-Px2,1) + norm(svect(Y(:,:,it+1))-Px3,1))/accelerant^2;
 +
disp(sprintf('straight = %g',straight));
 +
separation = norm(X(:,2) - X(:,1),1);
 +
disp(sprintf('separation %g',separation));
 +
if ~backout
 +
if straight < 1 && separation < 1
 +
accelerant = 2^-log10(straight);
 +
else
 +
accelerant = accelerant/2;
 +
end
 +
else
 +
disp('backing out')
 +
accelerant = accelerant/2;
 +
end
 +
if accelerant > 1
 +
disp(sprintf('accelerant %g',accelerant))
 +
elseif ~backout
 +
disp(' ')
 +
end
 +
if accelerant < 1, accelerant = 1; end
 +
Y(:,:,it+1) = Y(:,:,it+1) + svectinv((accelerant-1)*(X(:,2) - X(:,1)));
 +
end
 +
disp(' ')
 +
if length(find(coordinates)) <= 1
 +
disp(sprintf('iterations = %d',i))
 +
iavg = iavg + i;
 +
toc
 +
break
 +
end
 +
backout = false;
 +
end
 +
disp(sprintf('average iterations = %d',round(iavg/exemplar)))
 +
disp(sprintf('residual = %g',sum(sum(abs(A - U*diag(S)*V')))))
 +
disp(sprintf('singular value match = %g',sum(abs(sort(abs(S)) - sort(So)))))
 +
end
 +
</pre>
== fast max cut ==
== fast max cut ==
Line 805: Line 1,353:
written by Giovanni Rinaldi which can be found at
written by Giovanni Rinaldi which can be found at
[http://convexoptimization.com/TOOLS/RUDY http://convexoptimization.com/TOOLS/RUDY]
[http://convexoptimization.com/TOOLS/RUDY http://convexoptimization.com/TOOLS/RUDY]
-
together with graph data.
+
together with graph data. Requires [http://www.stanford.edu/~boyd/cvx CVX].
<pre>
<pre>
% fast max cut, Jon Dattorro, July 2007, http://convexoptimization.com
% fast max cut, Jon Dattorro, July 2007, http://convexoptimization.com
Line 830: Line 1,378:
traceXW = 1e15;
traceXW = 1e15;
while 1
while 1
-
cvx_begin % CVX Boyd
+
cvx_begin
variable X(N,N) symmetric;
variable X(N,N) symmetric;
maximize(trace(Q*X) - w*trace(W*X));
maximize(trace(Q*X) - w*trace(W*X));
Line 880: Line 1,428:
== Signal dropout problem ==
== Signal dropout problem ==
-
Requires CVX (Boyd).
+
Requires [http://www.stanford.edu/~boyd/cvx CVX].
<pre>
<pre>
%signal dropout problem -Jon Dattorro
%signal dropout problem -Jon Dattorro
Line 976: Line 1,524:
end
end
successes
successes
 +
</pre>
 +
 +
 +
== Compressive Sampling of Images by Convex Iteration <math>-</math> Shepp-Logan phantom ==
 +
<pre>
 +
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +
% Cardinality minimization by Convex Iteration.
 +
% -Jon Dattorro & Christine S. W. Law, July 2008.
 +
% Example. "Compressive sampling of a phantom" from Convex Optimization & Euclidean Distance Geometry.
 +
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +
clear all
 +
%%% Generate Shepp-Logan phantom %%%
 +
N = 256;
 +
image = chop(phantom(N));
 +
if mod(N,2), disp('N must be even for symmMap()'), return, end
 +
 +
%%% Generate K-space radial sampling mask with P lines %%%
 +
P = 10; %noninteger bends lines at origin.
 +
PHI = fftshift(symmMap(N,P)); %left justification of DC-centric radial sampling pattern
 +
 +
%%% Apply sampling mask %%%
 +
f = vect(real(ifft2(PHI.*fft2(image))));
 +
 +
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +
lambda = 1e8; % fidelity parameter on equality
 +
tol_inner = 1e-2; % inner loop exit tolerance
 +
tol_outer = 1e-5; % inner loop exit tolerance
 +
maxCGiter = 35; % max Conjugate Gradient iterations
 +
 +
%%% Build differential matrices %%%
 +
Psi1 = kron(speye(N) - spdiags(ones(N,1),-1,N,N), speye(N))';
 +
Psi2 = kron(speye(N), speye(N) - spdiags(ones(N,1),-1,N,N));
 +
Psi = [Psi1; Psi1'; Psi2; Psi2'];
 +
 +
close all
 +
figure, colormap(gray), axis image, axis off, set(gca,'position',[0 0 1 1]); %set(gcf,'position',[600 1300 512 512])
 +
 +
count = 0;
 +
u = vectinv(f);
 +
old_u = 0;
 +
last_u = 0;
 +
y = ones(4*N^2,1);
 +
card = 5092; % in vicinity of minimum cardinality which is 5092 at N=256; e.g., 8000 works
 +
tic
 +
while 1
 +
while 1
 +
count = count + 1;
 +
%%%%% r0 %%%%%
 +
tmp = Psi'*spdiags(y./(abs(Psi*u(:))+1e-3),0,4*N^2,4*N^2)*Psi;
 +
r = tmp*u(:) + lambda*vect(real(ifft2(PHI.*fft2(u - vectinv(f)))));
 +
 +
%%%%% Inversion via Conjugate Gradient %%%%%
 +
p = 0;
 +
beta = 0;
 +
for i=1:maxCGiter
 +
p = beta*p - r;
 +
Gp = tmp*p + lambda*vect(real(ifft2(PHI.*fft2(vectinv(p)))));
 +
rold = r'*r;
 +
alpha = rold/(p'*Gp);
 +
if norm(alpha*p)/norm(u(:)) < 1e-13
 +
break
 +
end
 +
u(:) = u(:) + alpha*p;
 +
r = r + alpha*Gp;
 +
beta = (r'*r)/rold;
 +
end
 +
if norm(u(:)) > 10*norm(f)
 +
disp('CG exit')
 +
disp(strcat('norm(u(:))=',num2str(norm(u(:)))))
 +
u = old_u;
 +
end
 +
%%%%% display intermediate image %%%%%
 +
imshow(real(u),[])
 +
pause(0.04)
 +
disp(count)
 +
 +
%%% test for fixed point %%%
 +
if norm(u(:) - old_u(:))/norm(u(:)) < tol_inner
 +
break
 +
end
 +
old_u = u;
 +
end
 +
%%% measure cardinality %%%
 +
err = 20*log10(norm(image-u,'fro')/norm(image,'fro'));
 +
disp(strcat('error=',num2str(err),' dB'))
 +
%%% new direction vector %%%
 +
[x_sorted, indices] = sort(abs(Psi*u(:)),'descend');
 +
cardest = 4*N^2 - length(find(x_sorted < 1e-3*max(x_sorted)));
 +
disp(strcat('cardinality=',num2str(cardest)))
 +
y = ones(4*N^2,1);
 +
y(indices(1:card)) = 0;
 +
 +
if norm(u(:) - last_u(:))/norm(u(:)) < tol_outer
 +
break
 +
end
 +
last_u = u;
 +
end
 +
toc
 +
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +
%%%% results %%%
 +
figure
 +
imshow(image,[])
 +
colormap(gray)
 +
axis image
 +
axis off
 +
 +
figure
 +
imshow(fftshift(PHI),[])
 +
colormap(gray)
 +
axis image
 +
axis off
 +
 +
figure
 +
imshow(abs(vectinv(f)),[])
 +
colormap(gray)
 +
axis image
 +
axis off
 +
 +
figure
 +
imshow(real(u),[])
 +
colormap(gray)
 +
axis image
 +
axis off
 +
 +
disp(strcat('ratio_Fourier_samples=',num2str(sum(sum(PHI))/N^2)))
 +
</pre>
 +
 +
=== symmMap() ===
 +
<pre>
 +
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +
% -Jon Dattorro July 2008
 +
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +
function SM2 = symmMap(N0,P)
 +
N = ceil(N0*sqrt(2)/2)*2;
 +
sampleMap = zeros(N,N);
 +
slice = zeros(N/2+1,2);
 +
rot_slice = zeros(N/2+1,2);
 +
dTheta = pi/P;
 +
for l=1:N/2+1
 +
x = N/2+1;
 +
y = l;
 +
slice(l,1) = x;
 +
slice(l,2) = y;
 +
end
 +
for p = 0:2*P
 +
theta = p*dTheta;
 +
for n = 1:N/2
 +
x = slice(n,1)-N/2-1;
 +
y = slice(n,2)-N/2-1;
 +
x2 = round( x*cos(theta)+y*sin(theta)+N/2+1);
 +
y2 = round(-x*sin(theta)+y*cos(theta)+N/2+1);
 +
sampleMap(y2,x2) = 1;
 +
end
 +
end
 +
SM = sampleMap(1:N,1:N);
 +
SM(N/2+1,N/2+1) = 1;
 +
 +
Nc = N/2;
 +
N0c = N0/2;
 +
 +
SM2 = SM(Nc-N0c+1:Nc+N0c, Nc-N0c+1:Nc+N0c);
 +
 +
% make vertically and horizontally symmetric
 +
[N1,N1] = size(SM2);
 +
Xi = fliplr(eye(N1-1));
 +
SM2 = round((SM2 + [ SM2(1,1) SM2(1,2:N1)*Xi;
 +
Xi*SM2(2:N1,1) Xi*SM2(2:N1,2:N1)*Xi])/2);
 +
</pre>
 +
 +
=== vect() ===
 +
<pre>
 +
function y = vect(A);
 +
y = A(:);
 +
</pre>
 +
 +
=== vectinv() ===
 +
<pre>
 +
%convert vector into matrix. m is dim of matrix.
 +
function A = vectinv(y)
 +
 +
m = round(sqrt(length(y)));
 +
if length(y) ~= m^2
 +
disp('dimension error in vectinv()');
 +
return
 +
end
 +
 +
A = reshape(y,m,m);
 +
</pre>
 +
<br>
 +
== High-order polynomials ==
 +
<pre>clear all
 +
E = zeros(7,7); E(1,5)=1; E(2,6)=1; E(3,4)=1; E = (E+E')/2;
 +
W1 = rand(7); W2 = rand(4);
 +
cvx_quiet('true')
 +
cvx_precision('high')
 +
while 1
 +
cvx_begin
 +
variable A(3,3) symmetric;
 +
variable C(6,6) symmetric;
 +
variable b(3);
 +
G = [A b; b' 1];
 +
G == semidefinite(4);
 +
X = [C [diag(A);b]; [diag(A)' b'] 1];
 +
X == semidefinite(7);
 +
sum(b) == 1;
 +
b >= 0;
 +
tr(X*E) == 4/27;
 +
minimize(trace(X*W1) + trace(G*W2));
 +
cvx_end
 +
 +
[v,d,q] = svd(G);
 +
W2 = v(:,2:4)*v(:,2:4)';
 +
rankG = sum(diag(d) > max(diag(d))*1e-8)
 +
 +
[v,d,q] = svd(X);
 +
W1 = v(:,2:7)*v(:,2:7)';
 +
rankX = sum(diag(d) > max(diag(d))*1e-8)
 +
 +
if (rankX + rankG) == 2, break, end
 +
end
 +
x = chop(G(1,4),1e-8), y = chop(G(2,4),1e-8), z = chop(G(3,4),1e-8)
</pre>
</pre>

Current revision

These MATLAB programs come from the book Convex Optimization & Euclidean Distance Geometry, by Dattorro, which is available for download online.

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

Some programs require an addon called CVX, an intuitive Matlab interface for interior-point method solvers.

Contents

isedm()

%Is real D a Euclidean Distance Matrix. -Jon Dattorro http://convexoptimization.com
%
%[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 corresponding to EDM output.
%Input: candidate matrix D, 
%       optional absolute numerical tolerance for EDM determination, 
%       optional verbosity 'on' or 'off',
%       optional desired affine dimension 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)
   V = 'Vn';
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(V,'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(V,'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

Alencar, Bonates, Lavor, & Liberti propose a more robust isedm() by replacing eigenvalue decomposition with singular value decomposition.

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

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


taxicab (1 norm) distance matrix

D1()

%Distance matrix from point list in 1 norm
function D = D1(X);

[n,N] = size(X);

D = kron(eye(N),ones(1,n))*abs(kron(ones(1,N),X(:)) - kron(ones(N,1),X));


conic independence

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.

conici()

% Remove c.i. directions in rows or columns of X.  -Jon Dattorro 
% function [Xci, indep_str, how_many_depend, kept] = conici(X,'rowORcol',tolerance);

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

if nargin < 3
   tol=max(size(X))*eps*normest(X);
end
if nargin < 2 || strcmp(rowORcol,'col') || isempty(rowORcol)
   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;
kept = [1:N]';
if rank(Xin,tol) == N 
   Xci = chop(X,tol); 
   return
end

count = 1;
new_N = N;
%remove zero columns
for i=1:new_N
   if chop(Xin(:,count),tol)==0
      how_many_depend = how_many_depend + 1;
      indep_str = 'conically Dependent';
      Xin(:,count) = [ ];
      kept(count) = [ ];
      new_N = new_N - 1;
   else
      count = count + 1;
   end
end
%remove identical columns
D = sqrt(Dx(Xin));
T = tril(ones(new_N,new_N))*normest(D);
ir = find(D+T < tol);
if ~isempty(ir)
    [iir jir] = ind2sub(size(D),ir);
    indep_str = 'conically Dependent';
    sizebefore = size(Xin);
    Xin(:,jir) = [ ];
    sizeafter = size(Xin);
    kept(jir) = [ ];
    new_N = size(Xin,2);
    how_many_depend = how_many_depend + sizebefore(2)-sizeafter(2);
end
%remove conic dependencies
count = 1;
newer_N = new_N;
for i=1:newer_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) = [ ];
         kept(count) = [ ];
         newer_N = newer_N - 1;
      else
         count = count + 1;
      end
   end
end
if strcmp(rowORcol,'col')
   Xci = chop(Xin, tol); 
else
   Xci = chop(Xin',tol);
end

lp()

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().

The source code is available here: lp.m which calls myqpsubold.m

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()

http://convexoptimization.com/TOOLS/USALO

mathworks.com/support/solutions/data/1-12MDM2.html?solution=1-12MDM2

% Find map of USA using only distance information.
% -Jon Dattorro
% Reconstruction from EDM.
clear all;
close all;

load usalo; % From Matlab Mapping Toolbox

% 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);

% 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

USA map input-data decimation, decimate()

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

EDM using ordinal data, omapusa()

http://convexoptimization.com/TOOLS/USALO

mathworks.com/support/solutions/data/1-12MDM2.html?solution=1-12MDM2

% Find map of USA using ordinal distance information.
% -Jon Dattorro
clear all;
close all;

load usalo; % From Matlab Mapping Toolbox

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

[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


Rank reduction subroutine, RRf()

% 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

svect()

% 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

svectinv()

% 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


Sturm & Zhang's procedure for constructing dyad-decomposition

This is a demonstration program that can easily be transformed to a subroutine for decomposing positive semidefinite matrix LaTeX: \,X\,.

This procedure provides a nonorthogonal alternative to eigen decomposition.

That particular decomposition obtained is dependent on choice of matrix LaTeX: \,A\,.

% 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)


Convex Iteration demonstration - Boolean feasibility

We demonstrate implementation of a rank constraint in a semidefinite Boolean feasibility problem.

It requires CVX, an intuitive Matlab interface for interior-point method solvers.

There are a finite number LaTeX: \,2^{N=_{}50}\!\approx 1E15\, of binary vectors LaTeX: \,x\,.

The feasible set of the semidefinite program from the book is the intersection of an elliptope with LaTeX: \,M\!=10\, halfspaces in vectorized composite LaTeX: \,G\,.

Size of the optimal rank-1 solution set is proportional to the positive factor scaling vector 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 A and vector b, are selected to demonstrate Boolean solution to one instance in a few iterations (a few seconds);

whereas sequential binary search takes one hour to test 25.7 million vectors before finding one Boolean solution feasible to the nonconvex problem from the book.

(Other parameters can be selected to realize a reversal of these timings.)

% 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);

while 1
    cvx_begin  
      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')

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

Convex Iteration rank-1

This program demonstrates how a semidefinite problem with a rank-r constraint is equivalently transformed into a problem sequence having a rank-1 constraint; discussed at the end of Chapter 4. Requires CVX.

%Convex Iteration Rank-1,  -Jon Dattorro  October 2007
% find        X
% subject to  A vec X = b
%             X positive semidefinite
%             rank X <= r

% X\in\symm^N
clear all
N=26;
M=10;
randn('seed',100);
A = randn(M,N*N);
b = randn(M,1);

r = 2;
W = ones(r*N);
Q11 = zeros(N,N);
Q22 = zeros(N,N);
count = 1;
traceGW = 1e15;
normof = 1e15;
w1 = 5;
cvx_quiet(true);
cvx_precision([1e-10 1e-4]);
while 1
    cvx_begin           
       variable L(r,1);
       X = L(1)*Q11 + L(2)*Q22;
       minimize(norm(A*X(:) - b));
       L >= 0;
    cvx_end

    % find Q
    cvx_begin
       variable Q11(N,N) symmetric;
       variable Q12(N,N);
       variable Q22(N,N) symmetric;
      
       trace(Q11) == 1;
       trace(Q22) == 1;
       trace(Q12) == 0;

       G = [Q11,  Q12;
            Q12', Q22];
       G == semidefinite(r*N);

       X = L(1)*Q11 + L(2)*Q22;

       minimize(w1*trace(G*W) + norm(A*X(:) - b));
    cvx_end

    [v,d,q] = svd(full(G));
    W = v(:,2:r*N)*v(:,2:r*N)';
    rankG = sum(diag(d) > max(diag(d))*1e-6)
    oldtraceGW = traceGW;
    traceGW = trace(G*W)
    oldnorm = normof;
    normof = norm(A*X(:) - b)
    if (rankG == 1) && (norm(A*X(:) - b) < 1e-6)
       break
    end

    if round((oldtraceGW - traceGW)*1e3) <= 0 && round((oldnorm - normof)*1e3) <= 0
        disp('STALLED');disp(' ')
        W = -v(:,2:r*N)*(v(:,2:r*N)' + randn(r*N-1,1)*v(:,1)'); 
    end
    count = count + 1;
end
count

Convex Iteration rank-1, 2013

%prototypical rank constrained sdp by rank-1 transformation,  -Jon Dattorro  October 2013
% find        X
% subject to  A vec X = b
%             X positive semidefinite
%             rank X <= r

% X is n x n symmetric
clc;
tachometer = 0;
accelcount = 0;
backoutcount = 0;
iavg = 0;
for exemplar = 1:1
   save exemplar exemplar iavg tachometer accelcount backoutcount
   clear all; 
   load exemplar
   rand('twister',sum(100*clock));
   randn('state',sum(100*clock));
   m = 25;                 %given matrix A is m x n(n+1)/2
   n = 7;
   r = 3;                  %assumed rank of matrix
   starting_window_length = 25;
   quant = 1e6;
   So = rand(r,1);
   Uo = orth(randn(n,r));
   Xo = Uo*diag(So)*Uo';            
   A = randn(m,n*(n+1)/2);
   b = A*svect2(Xo);
      
   Y = eye(n*r);
   cvx_quiet('true');
   tic
   accelerant = 1;
   iavgadj = 0;
   Z = [];
   backout = false;
   i = 0;
   it = 0;
   while 1
      i = i + 1;
      it = it + 1;
      if i > 1  &&  isempty(strfind(cvx_status,'Solved'))
         temp = cvx_solver;
         if ~isempty(strfind(temp,'SeDuMi')) 
            cvx_solver('SDPT3');
         else
            cvx_solver('SeDuMi');
         end
         iavgadj = iavgadj + 1;  
         temp = cvx_solver;
         disp(sprintf('switching solver to %s\n',temp));
      else
         cvx_solver('SeDuMi');
      end
      Zold = Z;
      U = sparse(n,r);
      T = [];
      k = 1;
      for j=1:r
         idx(k) = 1 + length(T);                                           % index Z matrix by block
         T = [T; U(:,j)];
         k = k + 1;
      end
      idx(k) = 1 + length(T);
      cvx_begin
         variable Z(n*r,n*r) symmetric;
         Z == semidefinite(n*r);
         XX = Z(idx(1):idx(2)-1,idx(1):idx(2)-1);                          % sum U_ii = X
         for j=2:r
            XX = XX + Z(idx(j):idx(j+1)-1,idx(j):idx(j+1)-1);
         end
         A*svect2(XX) == b;
         for j=1:r-1
            for l=j+1:r
               trace(Z(idx(j):idx(j+1)-1,idx(l):idx(l+1)-1)) == 0;         % u_i u_j orthogonality
            end
         end
         minimize(trace(Y(:,:,it)'*Z));
      cvx_end
      if it == 1, clc; end
      disp(sprintf('exemplar = %d',exemplar))
      disp(sprintf('cvx_status = %s',cvx_status))
      if exemplar > 1  
         disp(sprintf('average iterations = %d',round(iavg/(exemplar-1))));
      end
      [u s v] = svd(Z);
      if it > 1  &&  (isempty(strfind(cvx_status,'Solved'))  ||  (accelerant > 1  &&  sum(s(ambig+1:end)) > darkmatter(it-1)))  &&  ~tack
         it = it-1;
         backout = true;
         Z = Zold;
         [u s v] = svd(Z);
      end
      tack = false;
      temp = diag(s);
      coordinates = chop(temp(1:min(r*2,size(s,1))),quant^-1)
      ambig = max(1,length(find(abs(coordinates-coordinates(1)) < quant^-1)));
      Y(:,:,it+1) = u(:,ambig+1:end)*diag(sign(sum(u(:,ambig+1:end).*v(:,ambig+1:end))))*v(:,ambig+1:end)';
      Y(:,:,it+1) = (Y(:,:,it+1) + Y(:,:,it+1)')/2;
      darkmatter(it) = trace(Y(:,:,it+1)'*Z);
      disp(sprintf('traceYZ = %g',darkmatter(it)))
      if it > 1
         if it == 2, close all; end
         figure(1)
         if iavg, magi = floor(iavg/(exemplar-1)); else magi = starting_window_length; end
         if length(darkmatter) > magi
            plot(it-magi:it,darkmatter(end-magi:end));
         else
            plot(1:length(darkmatter),darkmatter);
         end
         set(gcf,'position',[400 430 256 256])
         pause(0.07);
      end
      if it > 2  % make a line specified by two points
         X(:,1) = svect(Y(:,:,it-1) + Y(:,:,it))/2;
         X(:,2) = svect(Y(:,:,it) + Y(:,:,it+1))/2;
         XVn = X*Vn(2);
         Px1 = X(:,1) + XVn*(XVn'*(svect(Y(:,:,it-1))-X(:,1)))/norm(XVn)^2;
         Px2 = X(:,1) + XVn*(XVn'*(svect(Y(:,:,it)) - X(:,1)))/norm(XVn)^2;
         Px3 = X(:,1) + XVn*(XVn'*(svect(Y(:,:,it+1))-X(:,1)))/norm(XVn)^2;
         straight = (norm(svect(Y(:,:,it-1))-Px1,1) + norm(svect(Y(:,:,it))-Px2,1) + norm(svect(Y(:,:,it+1))-Px3,1))/accelerant^2;
         disp(sprintf('straight = %g',straight));
         separation = norm(X(:,2) - X(:,1),1);
         disp(sprintf('separation %g',separation));
         if ~backout
            if straight < 1  &&  separation < 1
               accelerant = 2^-log10(straight);
            else
               accelerant = accelerant/2;
            end
         else 
            disp('backing out')
            accelerant = accelerant/2;
            backoutcount = backoutcount + 1;
         end
         if accelerant > 1 
            disp(sprintf('accelerant %g',accelerant))
            accelcount = accelcount + 1;
         elseif ~backout
            disp(' ')
         end 
         if accelerant < 1, accelerant = 1; end
         Y(:,:,it+1) = Y(:,:,it+1) + svectinv((accelerant-1)*(X(:,2) - X(:,1)));
      end      
      disp(' ')
      if length(find(coordinates)) <= 1
         [QQ DD VV] = svd(XX);
         eigens = diag(DD).*sign(sum(QQ.*VV)');
         if (sum(abs(A*svect(XX)-b)) <= quant^-1) && ~length(find(chop(eigens,quant^-1) < 0)) && (length(find(chop(eigens,quant^-1))) <= r)
            disp(sprintf('iterations = %d',i))
            iavg = max(0, iavg + i - iavgadj);  
            toc
            break
         end
      else
         if iavg  &&  ~mod(i,max(round(log10(quant)*iavg/(exemplar-1)),starting_window_length))  &&  separation < 1 
            temp = randn(n*r,n*r);
            Y(:,:,it+1) = temp*temp';
            iavgadj = iavgadj + max(round(log10(quant)*iavg/(exemplar-1)),starting_window_length);  
            tachometer = tachometer + 1;
            tack = true;
         elseif ~iavg  &&  ~mod(i,round(log10(quant)*starting_window_length))  &&  separation < 1 
            temp = randn(n*r,n*r);
            Y(:,:,it+1) = temp*temp';
            iavgadj = iavgadj + round(log10(quant)*starting_window_length); 
            tachometer = tachometer + 1;
            tack = true;
         elseif ambig > 1  &&  ~sum(abs(coordinates(ambig+1:end)))
            temp = randn(n*r,n*r);
            Y(:,:,it+1) = temp*temp';
            iavgadj = iavgadj + 1;  
            tachometer = tachometer + 1;
            tack = true;
         end
      end
      backout = false;
   end
   disp(sprintf('average iterations = %d',round(iavg/exemplar)))
   disp(sprintf('residual = %g',sum(abs(A*svect(XX)-b))))
   disp(sprintf('rankX = %g',length(find(chop(eigens,quant^-1)))))
   disp(sprintf('tack = %d percent',round(100*tachometer/exemplar)))
   if accelcount, disp(sprintf('backout/accelerant ratio = %g',backoutcount/accelcount)), end
end


svect2()

%Map from symmetric matrix to vector
% -Jon Dattorro 

function y = svect2(Y)
N = size(Y,1);
ind  = zeros(N*(N+1)/2,1);
beta = zeros(N*(N+1)/2,1);
count = 1;
for j=1:N
   for i=1:j
         ind(count) = sub2ind(size(Y),i,j);
      if i==j
         beta(count) = 1;
      else
         beta(count) = sqrt(2);
      end
      count = count + 1;
   end
end 
y = Y(ind).*beta;

dvect3()

%Map from EDM to vector
function y = dvect3(Y)
N = size(Y,1);
ind  = zeros(N*(N-1)/2,1);
beta = zeros(N*(N-1)/2,1);
count = 1;
for j=2:N
   for i=1:j-1
         ind(count) = sub2ind(size(Y),i,j);
      if i~=j
         beta(count) = 1;
      end
      count = count + 1;
   end
end 
y = Y(ind).*beta;

Singular Value Decomposition (SVD) by rank-1 transformation

Introduced at the end of Chapter 4 in 2014 book version. This Matlab program requires CVX.

%svd decomposition by rank-1 transformation  -Jon Dattorro  October 2013
% X = U diag(S) V'
clc;
iavg = 0;
for exemplar = 1:1
   save exemplar exemplar iavg
   clear all; 
   load exemplar
   rand('twister',sum(100*clock));
   randn('state',sum(100*clock));
   m = 6;                  %given matrix A is m x n, rank r
   n = 7;
   r = 3;                  %assumed rank of matrix
   
   quant = 1e6;
   Uo = orth(randn(m,r));
   So = rand(r,1);       
   Vo = orth(randn(n,r));
   A = Uo*diag(So)*Vo';    % A = U diag(S) V'.   Assume  H = U diag(S)
      
   Y = eye(2*m*r + n*r + r + 1);
   cvx_quiet('true');
   cvx_solver('sedumi');
   tic
   accelerant = 1;
   Z = [];
   backout = false;
   i = 0;
   it = 0;
   while 1
      i = i + 1;
      it = it + 1;
      Zold = Z;
      cvx_begin
         variable      H(m,r);
         variables U(m,r) S(r,1) V(n,r);
         T = [];
         k = 1;
         idx(k) = 1;    % index Z matrix by block
         k = k + 1;
         for j=1:r
            idx(k) = 2 + length(T);
            T = [T; H(:,j)];
            k = k + 1;
         end
         for j=1:r
            idx(k) = 2 + length(T);
            T = [T; U(:,j)];
            k = k + 1;
         end
         idx(k) = 2 + length(T);
         T = [T; S];
         k = k + 1;
         for j=1:r
            idx(k) = 2 + length(T);
            T = [T; V(:,j)];
            k = k + 1;
         end
         idx(k) = 2 + length(T);
         variable ZZT(length(T),length(T)) symmetric;
         Z = [1 T';
              T ZZT];
         for j=1:r
%           variables J58(n,n) J69(n,n) J710(n,n) symmetric;                                    % H U' symmetry
            dvect3(Z(idx(j+1):idx(j+2)-1,idx(r+j+1):idx(r+j+2)-1)) == dvect3(Z(idx(j+1):idx(j+2)-1,idx(r+j+1):idx(r+j+2)-1)');
         end
         Z == semidefinite(2*m*r + n*r + r + 1);
         S >= 0;
%        J512 + J613 + J714 == A;                                                               % H V' = A
         AA = Z(idx(2):idx(3)-1,idx(2*r+3):idx(2*r+4)-1);
         for j=2:r
            AA = AA + Z(idx(j+1):idx(j+2)-1,idx(2*r+j+2):idx(2*r+j+3)-1);
         end
         AA == A;
%        Z=[1        H(:,1)' H(:,2)' H(:,3)' U(:,1)' U(:,2)' U(:,3)' S'     V(:,1)' V(:,2)' V(:,3)'
%           H(:,1)   J55     J56     J57     J58     J59     J510    J511   J512    J513    J514  
%           H(:,2)   J56'    J66     J67     J68     J69     J610    J611   J612    J613    J614  
%           H(:,3)   J57'    J67'    J77     J78     J79     J710    J711   J712    J713    J714  
%           U(:,1)   J58'    J68'    J78'    J88     J89     J810    J811   J812    J813    J814 
%           U(:,2)   J59'    J69'    J79'    J89'    J99     J910    J911   J912    J913    J914  
%           U(:,3)   J510'   J610'   J710'   J810'   J910'   J1010   J1011  J1012   J1013   J1014   
%           S        J511'   J611'   J711'   J811'   J911'   J1011'  J1111  J1112   J1113   J1114 
%           V(:,1)   J512'   J612'   J712'   J812'   J912'   J1012'  J1112' J1212   J1213   J1214  
%           V(:,2)   J513'   J613'   J713'   J813'   J913'   J1013'  J1113' J1213'  J1313   J1314  
%           V(:,3)   J514'   J614'   J714'   J814'   J914'   J1014'  J1114' J1214'  J1314'  J1414];
%        H == [J811(:,1) J911(:,2) J1011(:,3)];                                                 % H = U diag(S)
         HH = [];
         for j=1:r
            HH = [HH Z(idx(r+j+1):idx(r+j+2)-1,idx(2*r+2)+j-1)]; 
         end
         HH == H;
%        trace(J58) == S(1,1);  trace(J69) == S(2,1);  trace(J710) == S(3,1);                   % H U' singular values
         for j=1:r
            trace(Z(idx(j+1):idx(j+2)-1,idx(r+j+1):idx(r+j+2)-1)) == S(j,1);
         end         
%        trace(J56) == 0;  trace(J57) == 0;  trace(J67) == 0;                                   % H orthogonality
         for j=2:r
            for l=j:r
               trace(Z(idx(j):idx(j+1)-1,idx(l+1):idx(l+2)-1)) == 0;
            end
         end
%        trace(J59) == 0; trace(J510) == 0;                                                     % U' H perpendicularity
%        trace(J68) == 0; trace(J610) == 0;
%        trace(J78) == 0; trace(J79)  == 0;
         for j=1:r
            for l=1:r
               if l~=j
                  trace(Z(idx(j+1):idx(j+2)-1,idx(r+l+1):idx(r+l+2)-1)) == 0;
               end
            end
         end
%        trace(J88)   == 1;  trace(J99)   == 1;  trace(J1010) == 1;                             %column normalization
%        trace(J1212) == 1;  trace(J1313) == 1;  trace(J1414) == 1;                                   
         for j=1:r
            trace(Z(idx(r+j+1):idx(r+j+2)-1,idx(r+j+1):idx(r+j+2)-1)) == 1;
            trace(Z(idx(2*r+j+2):idx(2*r+j+3)-1,idx(2*r+j+2):idx(2*r+j+3)-1)) == 1;
         end
%        trace(J89)   == 0;  trace(J810)  == 0;  trace(J910)  == 0;                             %orthogonality
%        trace(J1213) == 0;  trace(J1214) == 0;  trace(J1314) == 0;
         for j=2:r
            for l=j:r
               trace(Z(idx(r+j):idx(r+j+1)-1,idx(r+l+1):idx(r+l+2)-1)) == 0;
               trace(Z(idx(2*r+j+1):idx(2*r+j+2)-1,idx(2*r+l+2):idx(2*r+l+3)-1)) == 0;
            end
         end
%        trace(J55) == J1111(1,1);  trace(J66) == J1111(2,2);  trace(J77) == J1111(3,3);        % H inner product  
         for j=1:r
            trace(Z(idx(j+1):idx(j+2)-1,idx(j+1):idx(j+2)-1)) == Z(idx(2*r+2)+j-1,idx(2*r+2)+j-1);
         end
         minimize(trace(Y(:,:,it)'*Z));
      cvx_end
      if it == 1, clc; end
      disp(sprintf('exemplar = %d',exemplar))
      disp(sprintf('cvx_status = %s',cvx_status))
      if exemplar > 1  
         disp(sprintf('average iterations = %d',round(iavg/(exemplar-1))));
      end
      [u s] = signeig(Z);
      if it > 1  &&  (isempty(strfind(cvx_status,'Solved'))  ||  (accelerant > 1  &&  sum(s(ambig+1:end)) > darkmatter(it-1)))
         it = it-1;
         backout = true;
         Z = Zold;
         [u s] = signeig(Z);
      end   
      temp = diag(s);
      coordinates = chop(temp(1:min(r*2,size(s,1))),quant^-1)
      ambig = length(find(abs(coordinates-coordinates(1)) < quant^-1));
      Y(:,:,it+1) = u(:,ambig+1:end)*u(:,ambig+1:end)';
      darkmatter(it) = trace(Y(:,:,it+1)'*Z);
      disp(sprintf('traceYZ = %g',darkmatter(it)))
      if it > 1
         if it == 2, close all; end
         figure(1)
         if iavg, magi = floor(iavg/(exemplar-1)); else magi = 22; end   %22 is average number iterations for m=n=7 r=3
         if length(darkmatter) > magi
            plot(it-magi:it,darkmatter(end-magi:end));
         else
            plot(1:length(darkmatter),darkmatter);
         end
         set(gcf,'position',[400 430 256 256])
         pause(0.07);
      end
      if it > 2  % make a line specified by two points
         X(:,1) = svect(Y(:,:,it-1) + Y(:,:,it))/2;
         X(:,2) = svect(Y(:,:,it) + Y(:,:,it+1))/2;
         XVn = X*Vn(2);
         Px1 = X(:,1) + XVn*(XVn'*(svect(Y(:,:,it-1))-X(:,1)))/norm(XVn)^2;
         Px2 = X(:,1) + XVn*(XVn'*(svect(Y(:,:,it)) - X(:,1)))/norm(XVn)^2;
         Px3 = X(:,1) + XVn*(XVn'*(svect(Y(:,:,it+1))-X(:,1)))/norm(XVn)^2;
         straight = (norm(svect(Y(:,:,it-1))-Px1,1) + norm(svect(Y(:,:,it))-Px2,1) + norm(svect(Y(:,:,it+1))-Px3,1))/accelerant^2;
         disp(sprintf('straight = %g',straight));
         separation = norm(X(:,2) - X(:,1),1);
         disp(sprintf('separation %g',separation));
         if ~backout
            if straight < 1  &&  separation < 1
               accelerant = 2^-log10(straight);
            else
               accelerant = accelerant/2;
            end
         else 
            disp('backing out')
            accelerant = accelerant/2;
         end
         if accelerant > 1 
            disp(sprintf('accelerant %g',accelerant))
         elseif ~backout
            disp(' ')
         end 
         if accelerant < 1, accelerant = 1; end
         Y(:,:,it+1) = Y(:,:,it+1) + svectinv((accelerant-1)*(X(:,2) - X(:,1)));
      end      
      disp(' ')
      if length(find(coordinates)) <= 1
         disp(sprintf('iterations = %d',i))
         iavg = iavg + i;
         toc
         break
      end
      backout = false;
   end
   disp(sprintf('average iterations = %d',round(iavg/exemplar)))
   disp(sprintf('residual = %g',sum(sum(abs(A - U*diag(S)*V')))))
   disp(sprintf('singular value match = %g',sum(abs(sort(abs(S)) - sort(So)))))
end

fast max cut

We use the graph generator (C program) rudy written by Giovanni Rinaldi which can be found at http://convexoptimization.com/TOOLS/RUDY together with graph data. Requires CVX.

% 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                     
           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

Signal dropout problem

Requires CVX.

%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);

    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


Compressive Sampling of Images by Convex Iteration LaTeX: - Shepp-Logan phantom

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Cardinality minimization by Convex Iteration.
%   -Jon Dattorro & Christine S. W. Law, July 2008.
%   Example. "Compressive sampling of a phantom" from Convex Optimization & Euclidean Distance Geometry. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
%%% Generate Shepp-Logan phantom %%%
N = 256;
image = chop(phantom(N));
if  mod(N,2), disp('N must be even for symmMap()'), return, end

%%% Generate K-space radial sampling mask with P lines %%%
P = 10;                         %noninteger bends lines at origin.
PHI = fftshift(symmMap(N,P));   %left justification of DC-centric radial sampling pattern

%%% Apply sampling mask %%%
f = vect(real(ifft2(PHI.*fft2(image))));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lambda    = 1e8; 	% fidelity parameter on equality
tol_inner = 1e-2; 	% inner loop exit tolerance
tol_outer = 1e-5; 	% inner loop exit tolerance
maxCGiter = 35;   	% max Conjugate Gradient iterations

%%% Build differential matrices %%%
Psi1 = kron(speye(N) - spdiags(ones(N,1),-1,N,N), speye(N))';
Psi2 = kron(speye(N), speye(N) - spdiags(ones(N,1),-1,N,N));
Psi = [Psi1; Psi1'; Psi2; Psi2']; 

close all
figure, colormap(gray), axis image, axis off, set(gca,'position',[0 0 1 1]); %set(gcf,'position',[600 1300 512 512])

count = 0;
u = vectinv(f);
old_u = 0;
last_u = 0;
y = ones(4*N^2,1);
card = 5092;        % in vicinity of minimum cardinality which is 5092 at N=256; e.g., 8000 works
tic
while 1
    while 1
        count = count + 1;
        %%%%% r0 %%%%%
        tmp = Psi'*spdiags(y./(abs(Psi*u(:))+1e-3),0,4*N^2,4*N^2)*Psi;
        r = tmp*u(:) + lambda*vect(real(ifft2(PHI.*fft2(u - vectinv(f)))));

        %%%%% Inversion via Conjugate Gradient %%%%%
        p = 0;
        beta = 0;
        for i=1:maxCGiter
            p = beta*p - r;
            Gp = tmp*p + lambda*vect(real(ifft2(PHI.*fft2(vectinv(p)))));
            rold = r'*r;
            alpha = rold/(p'*Gp);
            if norm(alpha*p)/norm(u(:)) < 1e-13
                break
            end
            u(:) = u(:) + alpha*p;
            r = r + alpha*Gp;
            beta = (r'*r)/rold;
        end
        if norm(u(:)) > 10*norm(f)
            disp('CG exit')
            disp(strcat('norm(u(:))=',num2str(norm(u(:)))))
            u = old_u;
        end
        %%%%% display intermediate image %%%%%
        imshow(real(u),[])
        pause(0.04)
        disp(count)

        %%% test for fixed point %%%
        if norm(u(:) - old_u(:))/norm(u(:)) < tol_inner
            break
        end
        old_u = u;
    end
    %%% measure cardinality %%%
    err = 20*log10(norm(image-u,'fro')/norm(image,'fro'));
    disp(strcat('error=',num2str(err),' dB'))
    %%% new direction vector %%%
    [x_sorted, indices] = sort(abs(Psi*u(:)),'descend');  
    cardest = 4*N^2 - length(find(x_sorted < 1e-3*max(x_sorted)));
    disp(strcat('cardinality=',num2str(cardest)))
    y = ones(4*N^2,1);
    y(indices(1:card)) = 0;
    
    if norm(u(:) - last_u(:))/norm(u(:)) < tol_outer
        break
    end
    last_u = u;
end
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% results %%%
figure
imshow(image,[])
colormap(gray)
axis image
axis off

figure
imshow(fftshift(PHI),[])
colormap(gray)
axis image
axis off

figure
imshow(abs(vectinv(f)),[])
colormap(gray)
axis image
axis off

figure
imshow(real(u),[])
colormap(gray)
axis image
axis off

disp(strcat('ratio_Fourier_samples=',num2str(sum(sum(PHI))/N^2)))

symmMap()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   -Jon Dattorro July 2008
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function SM2 = symmMap(N0,P)
N = ceil(N0*sqrt(2)/2)*2;
sampleMap = zeros(N,N);
slice     = zeros(N/2+1,2);
rot_slice = zeros(N/2+1,2);
dTheta = pi/P;
for l=1:N/2+1
    x = N/2+1;
    y = l;
    slice(l,1) = x;
    slice(l,2) = y;    
end
for p = 0:2*P
   theta = p*dTheta;
   for n = 1:N/2
       x = slice(n,1)-N/2-1;
       y = slice(n,2)-N/2-1;
       x2 = round( x*cos(theta)+y*sin(theta)+N/2+1);
       y2 = round(-x*sin(theta)+y*cos(theta)+N/2+1);
       sampleMap(y2,x2) = 1;
    end
end
SM = sampleMap(1:N,1:N);
SM(N/2+1,N/2+1) = 1;

Nc  = N/2;
N0c = N0/2;

SM2 = SM(Nc-N0c+1:Nc+N0c, Nc-N0c+1:Nc+N0c);

% make vertically and horizontally symmetric
[N1,N1] = size(SM2);
Xi = fliplr(eye(N1-1));
SM2 = round((SM2 + [   SM2(1,1)        SM2(1,2:N1)*Xi;
                    Xi*SM2(2:N1,1)  Xi*SM2(2:N1,2:N1)*Xi])/2);

vect()

function y = vect(A);
y = A(:);

vectinv()

%convert vector into matrix.  m is dim of matrix.
function A = vectinv(y)

m = round(sqrt(length(y)));
if length(y) ~= m^2
   disp('dimension error in vectinv()');
   return
end

A = reshape(y,m,m);


High-order polynomials

clear all
E = zeros(7,7);  E(1,5)=1; E(2,6)=1; E(3,4)=1;  E = (E+E')/2;
W1 = rand(7);  W2 = rand(4);
cvx_quiet('true')
cvx_precision('high')
while 1
    cvx_begin
        variable A(3,3) symmetric;
        variable C(6,6) symmetric;
        variable b(3);
        G = [A b; b' 1];
        G == semidefinite(4);
        X = [C [diag(A);b]; [diag(A)' b'] 1];
        X == semidefinite(7);
        sum(b) == 1;
        b >= 0;
        tr(X*E) == 4/27;
        minimize(trace(X*W1) + trace(G*W2));
    cvx_end

    [v,d,q] = svd(G);
    W2 = v(:,2:4)*v(:,2:4)';
    rankG = sum(diag(d) > max(diag(d))*1e-8)

    [v,d,q] = svd(X);
    W1 = v(:,2:7)*v(:,2:7)';
    rankX = sum(diag(d) > max(diag(d))*1e-8)    

    if (rankX + rankG) == 2, break, end
end
x = chop(G(1,4),1e-8),  y = chop(G(2,4),1e-8),  z = chop(G(3,4),1e-8)
Personal tools