# Toeplitz block toeplitz matrices

[ Follow Ups ] [ Post Followup ] [ Netlib Discussion Forum ] [ FAQ ]

Posted by Graham French on June 04, 1998 at 09:34:01:

The equation I'm trying to solve is

H[i,j] = 1/4pi * M[i,j,k,l]g[k,l]

where i,j,k,l are indices running from 1..N where N is the number of
pixels in the image.

H is the measured field at a point i,j in the detector film. g is the
magnetisation at point [k,l] in the sample being measured. M is a
geometric function.

I need to determine g. One way is to invert this matrix equation i.e

g = M^-1 x H but this requires too much memory for a large image.
Therefore I've decided to solve iterativley.

For a given sample geometry M has only to be calculated once and stored
for subsequent use in the iterations.

As you can see, the amount of memory required to do this will be
proportional to N^4. Therefore at present I am calculating M at each
point in every iteration which slows the whole process quite
considerably.

According to the papers M is a Toeplitz block Toeplitz matrix with only
N^2 DIFFERENT elements. Hence it is possible for the memory requirements
to be reduced.

the problem is that I don't know how to map the N^2 different elements
to the full N^4 array for the calculation.

Here is M in full

-arctan((sqr(a)*(k-i+0.5)*(l-j+0.5))/((d+t)*sqrt(sqr(a*(k-i+0.5))+sqr(a*
(l-j+0.5))+sqr(d+t))))

+

arctan((sqr(a)*(k-i+0.5)*(l-j+0.5))/(d*sqrt(sqr(a*(k-i+0.5))+sqr(a*(l-j+
0.5))+sqr(d))))

+

arctan((sqr(a)*(k-i+0.5)*(l-j-0.5))/((d+t)*sqrt(sqr(a*(k-i+0.5))+sqr(a*(
l-j-0.5))+sqr(d+t))))

-

arctan((sqr(a)*(k-i+0.5)*(l-j-0.5))/(d*sqrt(sqr(a*(k-i+0.5))+sqr(a*(l-j-
0.5))+sqr(d))))

+

arctan((sqr(a)*(k-i-0.5)*(l-j+0.5))/((d+t)*sqrt(sqr(a*(k-i-0.5))+sqr(a*(
l-j+0.5))+sqr(d+t))))

-

arctan((sqr(a)*(k-i-0.5)*(l-j+0.5))/(d*sqrt(sqr(a*(k-i-0.5))+sqr(a*(l-j+
0.5))+sqr(d))))

-

arctan((sqr(a)*(k-i-0.5)*(l-j-0.5))/((d+t)*sqrt(sqr(a*(k-i-0.5))+sqr(a*(
l-j-0.5))+sqr(d+t))))

+

arctan((sqr(a)*(k-i-0.5)*(l-j-0.5))/(d*sqrt(sqr(a*(k-i-0.5))+sqr(a*(l-j-
0.5))+sqr(d))))

where i,j,k,l are the indices as before and a,d,t are constants.

For N=2, M has the following results:

k=1 k=2

i=1 i=2 i=1 i=2

j=1 7.63 -9.67 | -9.67 7.63
l=1 |
j=2 -9.67 -3.01 | -3.01 -9.67
--------------------------------|--------------------
j=1 -9.67 -3.01 | -3.01 -9.67
l=2 |
j=2 7.63 -9.67 | -9.67 7.63

As you can see there are only N^2 different elements, those in the top
left sub matrix.. all other sub matrices are reflections.

I was wondering if there is a way to only store a N^2 array such as the
top left one and then map those values for different [k,l].

If it is any help, apparantly M is a symmetric semipositive definite
matrix, and the matrix equation above can be solved iteratively using
the conjugate gradient method.

I would be extremely grateful for any help you can give on this one.

Many thanks
Graham French