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