added lecture5_6 & ex05

exercises/ex05/cs_ex.m

+% In this script one solves the MDR (missing data recovery) problem
+% A data vector f represented by a sparse (few nonzeros xs) expansion
+% in a frame (here the discrete cosine transform DCT) Phi. f experiences
+% a loss of information, its last n-m components are lost in b. The
+% system Phi(1:m,:)*xr = b that would need to be solved to recover the
+% coefficients xr is underdetermined and thus has infinitely many
+% solutions. Since it is known that xs is sparse, one tries to
+% find a sparse xr vector by solving the optimization problem
+%
+%     min ||x||_0  subject to   Ax = b, ||.||_0 denotes no. of nonzeros
+%
+% This is a combinatorial problem and NP hard. One replaces the ||.||_0
+% norm by the ||.||_1 norm which encourages sparsity when minimized.
+%
+%     (MDR)     min ||x||_1  subject to   Ax = b  (x free)
+%
+% This problem can be rewritten as standard form LP and solved by
+% Mehrotra's interior point method. Sparsity should be exploited
+% Use m = 400, 350, 300
+%
+n = 500;			% data size
+m = 400;			% available data size; may be changed
+Phi = dctmtx(n);		% get the frame
+xs = zeros(n,1);		% initialize x*
+k = 10;				% number of nonzeros in x*
+rand('state',1011);		% for reproducibility
+p = randperm(60);		% indices of nonzeros in x*
+randn('state',11);		% for reproducibility
+xs(p(1:k)) = randn(k,1);	% random values for x*
+f = Phi*xs;			% complete data vector v*
+b = f(1:m);			% available data b
+A = Phi(1:m,:);			% matrix A
+subplot(4,1,1);			% put all graphs in one figure
+plot(1:n,f);			% plot complete data
+subplot(4,1,2);
+plot(1:n,[b;zeros(n-m,1)]);	% plot available data
+%
+% fill in here; set up the LP and solve with Matlab's linprog
+%
+subplot(4,1,3);
+plot(1:n,Phi*xr);		% plot recovered data
+subplot(4,1,4);
+plot(1:n,Phi*xr-f);		% plot error

exercises/ex05/meh.m

+% IPM for LP, solving
+% min  c'*x
+% s.t. A*x=b
+% with A of size m times n
+function [x,ff]=meh(A,c,b,x,m,n);
+format long
+inf=10^20;
+eta=.99;
+tol=10^(-8);
+mu=1;
+maxit=100;
+s=1*ones(n,1);
+la=zeros(m,1);
+e=ones(n,1);
+dxaf=zeros(n,1);
+dlaf=zeros(m,1);
+dsaf=dxaf;
+X=diag(x);
+S=diag(s);
+L=[zeros(n) A' eye(n);A zeros(m) zeros(m,n);S zeros(n,m) X];
+for k=1:maxit
+   if (mu>tol)
+   iter=k;
+mu
+     f(k)=c'*x;
+     X=diag(x);
+     S=diag(s);
+     rb=A*x-b;
+     rc=A'*la+s-c;
+     XSe=x.*s;
+
+if 0 % fully sparse
+     L=[zeros(n) A' eye(n);A zeros(m) zeros(m,n);S zeros(n,m) X];
+     %[dxaf;dlaf;dsaf]=L\[-rc;-rb;-XSe];
+     ll=L\[-rc;-rb;-XSe];
+     dxaf=ll(1:n);
+     dlaf=ll(n+1:n+m);
+     dsaf=ll(n+m+1:2*n+m);
+else % augmented system
+	 KKT = [-diag(s./x) A'; A zeros(m)];
+%	 lll = KKT\[-rc+XSe./x;-rb];
+	 [LL, UU, pp] = lu(KKT, 'vector'); % KKT matrix is symmetric, but LDL is not implemented in Octave (yet)
+	 lll = [-rc+XSe./x;-rb];
+	 lll = UU\(LL\(lll(pp)));
+     dxaf=lll(1:n);
+     dlaf=lll(n+1:n+m);
+	 dsaf = (-XSe-s.*dxaf)./x;
+end
+
+if 1 % predictor-corrector IPM
+     mx=-x./dxaf;
+     ms=-s./dsaf;
+     alpaf=inf;
+     aldaf=inf;
+     for i=1:n
+        if (dxaf(i)<0)
+           alpaf=min(alpaf,mx(i));
+        end 
+        if(dsaf(i)<0)
+           aldaf=min(aldaf,ms(i));
+        end
+     end
+     alpaf=min(1,alpaf);
+     aldaf=min(1,aldaf);
+     muaf(k)=(x+alpaf*dxaf)'*(s+aldaf*dsaf)/n;
+     sigema(k)=(muaf(k)/mu)^3;
+     Xaf=diag(dxaf);
+     Saf=diag(dsaf);
+     dd=-XSe-Xaf*Saf*e+sigema(k)*mu*e;
+
+if 0 % fully sparse
+     % [dx;dla;ds]=L\[-rc;-rb;dd];
+     ll=L\[-rc;-rb;dd];
+     dx=ll(1:n);
+     dla=ll(n+1:n+m);
+     ds=ll(n+m+1:2*n+m);
+else % augmented system
+	 lll = [-rc-dd./x;-rb];
+	 lll = UU\(LL\(lll(pp)));
+     dx=lll(1:n);
+     dla=lll(n+1:n+m);
+	 ds = (dd-s.*dx)./x;
+end
+
+else % predictor IPM
+	 dx = dxaf;
+	 dla = dlaf;
+	 ds = dsaf;
+end
+
+     mx=-x./dx;
+     ms=-s./ds;
+     alpm=inf;
+     aldm=inf;
+     for i=1:n
+        if (dx(i)<0)
+           alpm=min(alpm,mx(i));
+        end 
+        if(ds(i)<0)
+           aldm=min(aldm,ms(i));
+        end
+     end
+     alpm=min(1,alpm);
+     aldm=min(1,aldm);   
+     alpk=min(1,eta*alpm);
+     aldk=min(1,eta*aldm);
+     x=x+alpk*dx;
+     ls=[la;s]+aldk*[dla;ds];
+     la=ls(1:m);
+     s=ls(m+1:m+n);
+     mu=x'*s/n;
+   end
+end
+iter
+ff=c'*x;
+%figure;
+%plot(f),title('f'),xlabel('# of iteration');