Using tikhonov regularization and LSQR to solve a linear set of equations

Hi all

I am trying to use Tikhonov regularization to minimize a linear set of equations. I take the generalized weighted minimization to be:

min( ||Ax-b||^2-lambda^2||Lx||^2 )      ,   [M,N]=size(A);

which can be formulated and solved in Matlab using LSQR (I typically increase the number of iterations):

x_estimate=lsqr([A;lambda*L],[b;0],[],1000);

Here the range of representative lambda, and subsequently the lambda representing the best weighting (reg_corner) between the residual norm |Ax-b|^2 and the other norm |Lx|^2, is determined using the function l_curve from regularization tools (RT: http://www2.imm.dtu.dk/~pcha/Regutools/)

[U,s,V] = csvd(A);
[reg_corner,~,~,lambda]=l_curve(U,s,b,'Tikh',L,V);
ix=find(lambda==reg_corner);
lambda=lambda(ix);

So far so good. This all works well with L=identity matrix.

I would then like to compare with other approaches such as including a smoothing requirement in place of the identity matrix (i.e. L-->smoothing) or the Total variation approach. I am having problems translating these concepts into the framework just explained though. Can anyone give me some insights on how to make these two approaches work in Matlab?

I have toyed around with a smoothing matrix, L, that was suggested to me, but it causes LSQR to have difficulties reaching to within the standard tolerance (1e-6). If d is the diagonal the matrix contains diagonal values of:

d-M=-1 , d-1=-1 , d=4 , d+1=-1 , d+M=-1

Can anyone tell me if this is indeed a good smoothing matrix or if there is an obvious reason why LSQR fails to converge properly for this matrix? Even if not converging, the results obtained using this L seem to mimic the true result much better than if L=identity though.

Also: are there better alternatives than using LSQR?

Thanks in advance!

Jakob