5.4 Inverse of L
5.4.3 Case 3: L is rectangular with fewer rows than columns
(AR−1)ˆξ−b2
2+λ2ξˆ2
2
, (43)
where we have the relationξˆ2
2 =kξk22. And ˆξ can be rewritten as
ξˆ=QTξ =QTLx=QTQRx=Rx ⇔ x=R−1ξ.ˆ
We can now solve expression (43) with respect to ˆξ and find the solution in terms of x: x=R−1ξˆ.
5.4.3 Case 3: L is rectangular with fewer rows than columns This section is based on theory derived by Per Christian Hansen [6].
In this case I considerL∈Rp×nbeing a rectangular matrix wherep < nand rank(L) =p. From the rank-nullity Theorem we know that these L-matrices have a non-empty null space. Therefore any Tikhonov solutionxto (40) can be splitted into two components: A component x0 ∈ N(L) that belongs to the null space of L and a component z that does not. Nowx equals z+x0 and the Tikhonov problem (40) becomes:
nkAz+Ax0−bk22+λ2kLzk22o.
Note that Lx0 = 0 since x0 ∈ N(L). By requiring that Az is orthogonal to Ax0 we can obtain two independent problems that can be solved with respect toz and x0:
minz
nkAz−bk22+λ2kLzk22o min
x0
kAx0−bk22.
Since Az ⊥ Ax0 we can split Rn into the two components N(L) and N(L)⊥A so that z∈ N(L)⊥A and x0 ∈ N(L).
The derivation of the pseudo inverse of L in this case is very technical and requires defi-nitions about oblique projectors and oblique pseudo inverses [6]. Therefore I will skip this part and only mention information most relevant to my project.
By introducing oblique projectors associated to N(L)⊥A and N(L) we can rewrite our minimization problems. Using the matrixW that spans the null space ofLwe can obtain an expression ofx0:
x0 =W(AW)†b.
From this definition of x0 we can derive an expression for the oblique projector that is associated to the null space ofL. Using this information we end up with an expression for L#such that our minimization problem for Tikhonov is on the form:
minξ
(AL#)ξ−b2
2+λ2kξk22
. (44)
The expression forL# will be given by:
L#= (I−W(AW)†A)L†,
where L† depends on the problem and the structure of L. In this case we consider L having full rank and therefore L can be splitted into L = (Llef t, Lright), where Llef t is square and has full rank. Thus,L† will in this case be on the form:
L†= L−1lef t 0
! .
And an expression for (AW)† can be obtained by QR-factorization. Decomposing AW intoQAW and RAW we get: AW =QAWRAW. In case 2 we derived the Moore-Penrose pseudo inverse for the matrix L which was given by L† = R−1QT. The same principle will work in this case and thus, (AW)† = R−1AWQTAW. Now we arrive at the following expression forL#:
L#= (I−W(R−1AWQTAW)A) L−1lef t 0
! .
We can now solve the Tikhonov problem (44) with respect toξ and find the solution in terms ofx: x=L#ξ.
For the second derivative matrix in two dimensions, L22, we encounter problems when separating the matrix intoLlef t andLright and inverting the quadratic matrixLlef t, since L22 has a structure as shown in the left part of figure 6. We see that most elements in the diagonal and upper diagonal are nonzero, while some elements are equal to zero.
ThusLlef t will be rank-deficient and have no finite inverse. ConsideringN discretization intervals for x, the zero elements appear in the diagonal as number N2−N −1 and in the upper diagonal as numberN2−N. To get rid of the zero elements in the diagonals we introduce a permutation matrix Π that repositions column numberN2 −N −1 and N2−N to the last columns in the matrix, while the lastN2−N + 1 columns are moved two columns to the left. This proces is illustrated in figure 6. I define L= LΠ and can now separate this matrix intoLlef t and Lright. L† will here be given by:
L†= L−1lef t 0
! .
5 Priorconditioner matrix L Iterative tomographic reconstruction To obtain the pseudo inverse with respect to L we permute L†: L† = ΠL†. Using this expression forL†we can follow the same principle as forL-matrices that didn’t have zeros in the diagonal to obtain the oblique pseudo inverseL#:
L#= (I −W(RAW−1 QTAW)A)Π L−1lef t 0
! .
This approach opens the possibility to work with the matrixL22. Note that this only works for this specificL-matrix. Further approaches could aim toward finding a procedure that works for arbitraryL-matrices.
0 10 20 30 40 50 60 70 80 90 100
nz = 1791 0
10
20
30
40
50
60
70
80
90
0 10 20 30 40 50 60 70 80 90 100
nz = 1791 0
10
20
30
40
50
60
70
80
90
Figure 6: Nonzero elements in the L22 matrix to the left and in the L†-matrix obtained by using a permutation matrix to the right.
6 Matlab implementation
The methods PCimmino and PKaczmarz are implemented in Matlab namedPCimmino.m andPKaczmarz.m. These methods use the functionget_l.m from the AIR Tools package to define the L-matrices for one-dimensional problems and the function get_L2.m to define theL-matrices for two-dimensional problems. Multiplications with the inverse and the transpose of the inverse of L in one dimension are performed using the functions lsolve.mand ltsolve.mfrom the AIR Tools package. For two dimensions Per Christian Hansen added modifications for the second derivative matrix as derived in section 5.4.3, resulting in the functionslsolve2D.mand ltsolve2D.m.
Since I don’t use large test problems, where the resolution is high, I define the matrix L#L#TAT explicitly to save computation time. A more general approach would be to calculateL#L#TAT and L#L#TrTi for PCimmino and PKaczmarz respectively in each iteration to save memory space. In this case it would be possible to use the method on larger test problems than possible with my approach. Note that for the last approach the computation time will be K times longer than for the my approach. Here K is the number of iterations for PCimmino and the number of sweeps over the rows of AL# for PKaczmarz. Also note that longer computation time would especially be a disadvantage for PCimmino since PCimmino, particularly for low noise levels, needs many iterations to obtain a good reconstruction. A more efficient way could therefore be to write code that uses the first approach for small and the last approach for large test problems.
As Matlab is optimized for matrix and vector operations but is slow at performing loops, PKaczmarz performs really slow since we need aforloop to sweep over the rows ofAL# and a for loop to perform a given number of iterations. Therefore the code could be improved by writing in a programming language that supports loops more effectively.
7 Performance of PKaczmarz & PCimmino Iterative tomographic reconstruction
7 Performance of PKaczmarz & PCimmino
To test the performance of PKaczmarz and PCimmino we created and used different test problems in one and two dimensions. In the following I name Kaczmarz’s and Cimmino’s method the standard version of PKaczmarz and PCimmino. In general the solution for a reconstruction problem in tomography is a two-dimensional image that represents a slice of a three-dimensional object. Therefore the main purpose of the one-dimensional problems is highlighting the difference between the priorconditioned methods and their standard versions. Another advantage is the simplicity of these problems which makes it easier to analyze the effect of using a priorconditioner matrix.
7.1 Test problems from the "AIR Tools" package
The following test problems are created using the "AIR Tools" package in Matlab [5]. Note that in all cases the datavectorb is obtained by multiplying the coefficient matrix A by the known solution x and adding an amount of noise; this phenomena is called "inverse crime" since the same model is used to compute the reconstruction and to get test data [6]. Therefore the following test problems may not be representative for how PKaczmarz and PCimmino operate on real data. However the purpose of the following tests is to give an impression of their performance; if the methods don’t yield good results for synthetic data they won’t yield good results for real data either.