Path: news.mathworks.com!not-for-mail
From: "sheng fang" <maelstromer@gmail.com>
Newsgroups: comp.soft-sys.matlab
Subject: unstable and un-convergence total variation regularization (with the code)
Date: Fri, 3 Oct 2008 03:49:01 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 56
Message-ID: <gc44nd$ejf$1@fred.mathworks.com>
Reply-To: "sheng fang" <maelstromer@gmail.com>
NNTP-Posting-Host: webapp-02-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1223005741 14959 172.30.248.37 (3 Oct 2008 03:49:01 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Fri, 3 Oct 2008 03:49:01 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1412410
Xref: news.mathworks.com comp.soft-sys.matlab:493421


HI,everyone:

I am working to solve a linear ill-posed problem Au=f using total variation regularization as follows:

min || A*u - f ||^2_L^2  +  lambda*TV(u)

However, the TV regularized solution is not stable. When the lambda is large, the solution contains large oscillation. When the lambda is small, the solution remains the same as the direct matrix inversion solution,i.e. u=f\A. Besides, the solution doesn't converge. Whatever lambda is, it is either one of the aforemention cases.

This really puzzles me, because the same code works well for denoising. will anyone give me some hints?
Thanks a lot!

Code:

% solve the WING problem using TV regularization:
% min || A*u - f ||^2_L^2  +  lambda*TV(u)
%% setup problem
n=64;% length of tested data.
t1=1/3;
t2=2/3;

% Set up matrix A.
A = zeros(n,n); h = 1/n;
sti = ((1:n)-0.5)*h;
for i=1:n
  A(i,:) = h*sti.*exp(-sti(i)*sti.^2);
end
% set up the idea solution;
I = find(t1 < sti & sti < t2);
u0 = zeros(n,1); u0(I) = sqrt(h)*ones(length(I),1);
% setup the right-hand 
f = sqrt(h)*0.5*(exp(-sti*t1^2)' - exp(-sti*t2^2)')./sti';

%% solve by TV regularized problem
% parmeters setup
ep2 = 1e-3;
dt = 0.02; % time step
lambda = 1;
nx = size(f,1);
NumSteps=400; % iteration number

u=zeros(size(f));% set initial value

for i=1:NumSteps,  
   % estimate derivatives
	u_x = (u([2:nx nx],:)-u([1 1:nx-1],:))/2;
	u_xx = u([2:nx nx],:)+u([1 1:nx-1],:)-2*u;
   % compute flow
    Num = ep2.*u_xx;
    Den = (ep2+u_x.^2).^(3/2);
    u_tv = Num./Den;
    u_fidelity = 2*A'*(f-A*u);
   % evolve image by dt
    u=u+dt*(lambda*u_tv+u_fidelity);  
end 
%% display
figure(2);plot(u,'r');hold on;plot(u0,'-.k');hold off;legend('tv','original');