Regularization with and total variation and regularizer using the CVX toolkit and automatic figure saving.

if ~exist('cvx_begin')
error('You first have to run cvx_setup!')
end

clear all % clears all variables/memory
close all % closes all figures

%%%% HERE YOU CAN SET PARAMETER VALUES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Reg=1: \|\nabla f\|_2^2 , Reg=0: \|\nabla f\|_1, Reg=2: \|\Delta f\|_2^2
REG=2; SIN=1;lambda=0.01; % Two Functions: SIN=0, step function, SIN=1, x^2 sin(x)

%to create all figures automatically
for REG=[0 1 2]
for SIN=[0 1]
for lambda=[10^-5 10^-4 10^-3 0.01 0.1 1]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% do not modify this part for the submission
% Generates training data
randn('state',0);%sum(100*clock));
m=200; n=4;
points = 2*rand(m,1)-1;

if(SIN)
Y = (points+0.5).^2.*sin(points) + 0.2*randn(m,1);
else
Y=(points>0) + 0.2*randn(m,1);
end

% generates the interpolation matrix
h=0.01; x=-1.1: h :1.1; % the grid
InterP=sparse(length(Y),length(x));
for i=1:length(Y)
index=find(points(i)>x);
pos1=max(index);
weight1=(x(pos1+1)-points(i))/h;
weight2=(points(i)-x(pos1))/h;
InterP(i,pos1)=weight1;
InterP(i,pos1+1)=weight2;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%This is the part which you should modify %%%%%%%%%%%%%%%%%%%%%%%%%

% generate the discretization of the first derivative
n=length(x); % number of variables in the problem
D=spdiags([-ones(n,1),ones(n,1)],[-1 0],n,n)/h;
D(1,:)=zeros(1,n); D(end,:)=zeros(1,n); %eliminates the underspecified boundary

Ds=spdiags([ones(n,1),-2.*ones(n,1),ones(n,1)],[-1 0 1],n,n)/h.^2;
%Ds(1,:)=zeros(1,n);
Ds(:,1)=zeros(1,n);
Ds(end,:)=zeros(1,n); %eliminates the underspecified boundary
%Ds(:,end)=zeros(1,n);

% D is the matrix for computing the first derivative

cvx_begin
variable f(n)
if(REG==1)
minimize( (InterP*f-Y)'*(InterP*f-Y) + lambda * (D*f)'*(D*f));
elseif(REG==0)
minimize( (InterP*f-Y)'*(InterP*f-Y) + lambda * norm(D*f,1));
else %REG==2
minimize( (InterP*f-Y)'*(InterP*f-Y) + lambda * (Ds*f)'*(Ds*f));
end
cvx_end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot solutions
figure
hold on
plot(x(2:end-1),f(2:end-1),'r-.','LineWidth',1.5);

pointsTest=-1.1:0.01:1.1;
plot(points,Y,'k.','MarkerSize',5); % Training data
if(SIN)
plot(pointsTest,(pointsTest+0.5).^2.*sin(pointsTest),'k-','LineWidth',1);
else
plot(pointsTest,(pointsTest>0),'k-','LineWidth',1);
end
h=legend('fitted function','Training data','True function');
set(h,'Location','NorthWest');
set(h,'interpreter','latex');
set(h,'FontSize',10);

grid on
if(REG==1)
if(SIN)
h=title(['Fitting: $x^2\sin(x)$ - Linear Spline: $\|\nabla f\|_2^2, \lambda=$',num2str(lambda)],'FontSize',10)
else
h=title(['Fitting step function - Linear Spline: $\|\nabla f\|_2^2, \lambda=$',num2str(lambda)],'FontSize',10)
end
elseif(REG==0)
if(SIN)
h=title(['Fitting: $x^2\sin(x)$ - Total Variation: $\|\nabla f\|_1, \lambda=$',num2str(lambda)],'FontSize',10)
else
h=title(['Fitting step function - Total Variation: $\|\nabla f\|_1, \lambda=$',num2str(lambda)],'FontSize',10)
end
else %REG==2
h=title(['Fitting: f(x) - Total Variation: $\|\Delta f\|_2^2, \lambda=$',num2str(lambda)],'FontSize',10)
end
set(h,'Interpreter','latex')
hold off

fOut = sprintf('totalVarReg-REG%d-SIN%d-lambda%0.5g.png',REG,SIN,lambda);
print('-dpng',fOut);

end;
end;
end;