Thread Subject: How to add a total variation term to my tomographic reconstruction problem?

Hi there,

I have got a very simple model to do a tomographic reconstruction:

arg min_x ||Ax - y||^2

in which, A is the system matrix and y is the true projections.

We can solve the problem using Matlab Quasi-Newton (fminunc) or lsqr function; however, may I ask how to add a total variation penalty to the objective function please?

arg min_x ||Ax - y||^2 + TV(x)

Is there any exsisting implementation in Matlab I can use?

For my understanding, there are two way to construct this TV term:
1.
[D1,D2,D3] = gradient(x);
TV = sum((D1.^2+D2.^2+D3.^2).^.5);

2.
D = del2(x);
TV = sum(abs(D));

For example, if my cost function is:

C = sum(sum(sum((Ax-y).^2)));

should I just add TV term to my cost function directly? Like:

C_TV = sum(sum(sum((Ax-y).^2))) + TV;

If I use 'fminunc' to perform the optimisation, then I want to provide the gradient as well. I have got the gradient of my cost function C quite straightforward, but how to get the gradient of TV let's say TV'?

Cheers,
Aaronne

"Aaronne" wrote in message <is82ro$r6m$1@newscl01ah.mathworks.com>...
> Hi there,
>
> I have got a very simple model to do a tomographic reconstruction:
>
> arg min_x ||Ax - y||^2
>
> in which, A is the system matrix and y is the true projections.
>
> We can solve the problem using Matlab Quasi-Newton (fminunc) or lsqr function; however, may I ask how to add a total variation penalty to the objective function please?
>
> arg min_x ||Ax - y||^2 + TV(x)
>
> Is there any exsisting implementation in Matlab I can use?
>
> For my understanding, there are two way to construct this TV term:
> 1.
> [D1,D2,D3] = gradient(x);
> TV = sum((D1.^2+D2.^2+D3.^2).^.5);
>
> 2.
> D = del2(x);
> TV = sum(abs(D));

That's not my understanding. Normally, I see TV defined as follows

[D1,D2,D3] = gradient(x);
Dx=[ D1(:) ; D2(:) ; D3(:) ];
TV = norm( Dx , 1);

As defined this way, TV is not differentiable whenever any D1(i)=0 and similarly for D2 and D3. So, you can't take its gradient there or rely on fminunc, etc...

You could make a differentiable approximation to TV

TV = sum((D1.^2+D2.^2+D3.^2+delta).^.5);

for some small tuning parameter delta>0

You could also reformulate your cost function as a constrained problem by introducing the unknown vector of variables z


C_TV = sum(sum(sum((Ax-y).^2))) + sum(z);

-z<=Dx
z<=Dx


and then use FMINCON.
 
> If I use 'fminunc' to perform the optimisation, then I want to provide the gradient as well. I have got the gradient of my cost function C quite straightforward, but how to get the gradient of TV let's say TV'?
=======================

As I said, the TV function is not differentiable. However, whichever of the above solultions you choose, you will need to recognize that the expression Dx from above can be implemented as a sparse matrix multiplication

Dx=C*x(:)

where C is a sparse differencing matrix. You will need to compute C and its transpose in whatever you do. For example if x is a 5x5 image, then


>> n=5; Diff= diag(ones(1,n-1),1)-eye(n); Diff(end,:)=[]

Diff =

    -1 1 0 0 0
     0 -1 1 0 0
     0 0 -1 1 0
     0 0 0 -1 1


>> C=[kron(Diff,speye(5)) ; kron(speye(5),Diff)];

Hi Matt,

Thanks a lot for your reply! And your point is really useful that TV is not differentiable; however, for my understanding you give the differentiable approximation to TV, which is called Huber function regression (Or maybe I am wrong please correct me.)

Anyway, we can make a differentiable approximation to TV:

% The x is an image
delta = 0.05;
[D1,D2,D3] = gradient(x);
TV = sum((D1.^2+D2.^2+D3.^2+delta).^.5);

% Cost function
C_TV = sum(sum(sum((Ax-y).^2))) + TV;

Maybe I am stupid but I still can't see how to use the TV approximation to get the differentiation of TV. And I have two questions now:

1. How to differentiate TV? Using the form:
TV = sum((D1.^2+D2.^2+D3.^2+delta).^.5);

2. I don't understand the Dx=C*x(:) you are talking about.
I think it is a practical way to find the gradient of the image x, but if we can use
[D1,D2,D3] = gradient(x);
Dx=[ D1(:) ; D2(:) ; D3(:) ];
Then why do we need to calculate Dx by Dx=C*x(:)?
And also how does it related to the differentiation of TV?

At this stage I want to perform my tomographic reconstruction using 'fminunc'. And I may also look into 'fmincon' later.

I really appreciate all your help! Thanks very much and have a nice day.


Best Regards,
Aaronne


> That's not my understanding. Normally, I see TV defined as follows
>
> [D1,D2,D3] = gradient(x);
> Dx=[ D1(:) ; D2(:) ; D3(:) ];
> TV = norm( Dx , 1);
>
> As defined this way, TV is not differentiable whenever any D1(i)=0 and similarly for D2 and D3. So, you can't take its gradient there or rely on fminunc, etc...
>
> You could make a differentiable approximation to TV
>
> TV = sum((D1.^2+D2.^2+D3.^2+delta).^.5);
>
> for some small tuning parameter delta>0
>
> You could also reformulate your cost function as a constrained problem by introducing the unknown vector of variables z
>
>
> C_TV = sum(sum(sum((Ax-y).^2))) + sum(z);
>
> -z<=Dx
> z<=Dx
>
>
> and then use FMINCON.
>
> > If I use 'fminunc' to perform the optimisation, then I want to provide the gradient as well. I have got the gradient of my cost function C quite straightforward, but how to get the gradient of TV let's say TV'?
> =======================
>
> As I said, the TV function is not differentiable. However, whichever of the above solultions you choose, you will need to recognize that the expression Dx from above can be implemented as a sparse matrix multiplication
>
> Dx=C*x(:)
>
> where C is a sparse differencing matrix. You will need to compute C and its transpose in whatever you do. For example if x is a 5x5 image, then
>
>
> >> n=5; Diff= diag(ones(1,n-1),1)-eye(n); Diff(end,:)=[]
>
> Diff =
>
> -1 1 0 0 0
> 0 -1 1 0 0
> 0 0 -1 1 0
> 0 0 0 -1 1
>
>
> >> C=[kron(Diff,speye(5)) ; kron(speye(5),Diff)];

"Aaronne" wrote in message <isabd0$gaa$1@newscl01ah.mathworks.com>...
>
> Maybe I am stupid but I still can't see how to use the TV approximation to get the differentiation of TV. And I have two questions now:
>
> 1. How to differentiate TV? Using the form:
> TV = sum((D1.^2+D2.^2+D3.^2+delta).^.5);
>
> 2. I don't understand the Dx=C*x(:) you are talking about.
> I think it is a practical way to find the gradient of the image x, but if we can use
> [D1,D2,D3] = gradient(x);
> Dx=[ D1(:) ; D2(:) ; D3(:) ];
> Then why do we need to calculate Dx by Dx=C*x(:)?
> And also how does it related to the differentiation of TV?

The gradient of a function f(C*x) where C is a matrix and f is any differntiable vector-valued function,

  f:R^m-->R^n

can be obtained from the chain rule as C.'*g(C*x) where g(z) is the gradient of f(z) with respect to z.

Applying this to your case, the expression z=C*x is equivalent to

[D1,D2,D3]=gradient(x);
z=[D1(:),D2(:),D(3:)];

and g(z)= z/norm(z+delta).

Unfortunately, however, MATLAB doesn't give you a function equivalent to the operation C.'*y, so you have no clear way of completing the evaluation of the desired expression C.'g(C*x). You have to somehow re-express the gradient() function in matrix form...

"Matt J" wrote in message <isahpt$25q$1@newscl01ah.mathworks.com>...
>
> and g(z)= z/norm(z+delta).
>

Sorry, that should be

g(z)=z/norm([z;sqrt(delta)]);

"Matt J" wrote in message <isahpt$25q$1@newscl01ah.mathworks.com>...
>
> Unfortunately, however, MATLAB doesn't give you a function equivalent to the operation C.'*y, so you have no clear way of completing the evaluation of the desired expression C.'g(C*x). You have to somehow re-express the gradient() function in matrix form...
====================

As a matter of personal interest, I tested below 2 different ways of computing TV and its gradient for a 128x256x64 image, x.

One method implements the differencing operator C in sparse numeric form. The other method decomposes C into components Cx, Cy, and Cz, and implements these components using my KronProd class

http://www.mathworks.com/matlabcentral/fileexchange/25969-efficient-object-oriented-kronecker-product-manipulation

It turns out (see below) that putting C in sparse numeric form is faster by a few factors, however, it requires 155MB of memory whereas using KronProd cuts it down to about 0.75MB. I don't know how big your actual images are, or how memory-conservative you were hoping to be...

%% data
x=rand(128,256,64);

delta=.01;

tv=@(z) norm([z(:);sqrt(delta)]);

[M,N,P]=size(x);

%1D gradient operators along in numeric form
[~,Ax]=gradient(speye(M));
[~,Ay]=gradient(speye(N));
[~,Az]=gradient(speye(P));


%3D gradient operators along x,y,z as KronProd objects
 Cx=KronProd({Ax,1,1},1:3,[M,N,P]);
 Cy=KronProd({1,Ay,1},1:3,[M,N,P]);
 Cz=KronProd({1,1,Az},1:3,[M,N,P]);

 
%convert the gradient operators to sparse numeric form and stack
 C=[sparse(Cx);sparse(Cy);sparse(Cz)];


%% Test 1 --- using sparse numeric C

tic;

  z=C*x(:);
  
  TV=tv(z);
  
  gradientTV = C.'*z /TV;
   
     gradientTV =reshape( gradientTV, [M,N,P]);
  
toc;%Elapsed time is 0.265185 seconds.



%% Test 2 --- using KronProd

tic;

  z1=Cx*x;
  z2=Cy*x;
  z3=Cz*x;
 
  TV =tv(cat(4,z1,z2,z3));
  
  gradientTV =(Cx.'*z1 + Cy.'*z2 +Cz.'*z3)/TV;

toc; %Elapsed time is 0.800754 seconds.



%%memory consumption

>> Whos C Cx Cy Cz
  Name Size Kilobytes Class Attributes
                                                                
  C 6291456x2097152 155649 double sparse
  Cx 2097152x2097152 129 KronProd
  Cy 2097152x2097152 513 KronProd
  Cz 2097152x2097152 33 KronProd