## The System Matrix

We can now summarize the assumptions in the two previous sections. Putting the image model (Eq. (26)) into the data model (Eq. (22)) reduces the problem to a set of linear equations:

where the elements of the system matrix are aji = T j drbi(r)¥j(r) j = 1,...,NLm;i = 1,...,P

A line integral model including only attenuation correction and normalization generates a sparse system matrix a with elements simple enough to be calculated on the fly. More accurate models that include scatter lead to densely populated matrices, which are complex to calculate. A practical algorithm then requires a compromise between accuracy, required storage, and speed. A useful approach is to factor a as a product of matrices, each of which models a specific aspect of the data acquisition .

A direct inversion of the linear system (Eq. (28)) with the < pj > replaced by the measured data pj is impractical for two reasons:

• The discrete system is ill-conditioned: the condition number of a is large(10). Consequently, the solution(11) of Eq. (28) is unstable for small perturbations pj- < pj > of the data. Ill-conditioning is the discrete equivalent of the ill-posedness of the inverse x-ray transform discussed in the section on the ill-posedness of the inverse X-ray transform (above).

• Numerically, the inversion of matrix a is hindered by its very large size (typically P = 106 unknowns and Nlor ^ 106 up to Nlor ^ 109 in 3D PET).

The first problem is solved by incorporating prior knowledge in a cost function. The second, numerical problem, is solved by optimizing the cost function by successive approximations.