Matlab for Convex Optimization & Euclidean Distance Geometry
From Wikimization
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.
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 .
This procedure provides a nonorthogonal alternative to eigen decomposition.
That particular decomposition obtained is dependent on choice of matrix .
% 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 of binary vectors
.
The feasible set of the semidefinite program from the book is the intersection of an elliptope
with halfspaces in vectorized composite
.
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 
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)