Depending on the number of projections collected, the system of equations to be solved may be either over- or underdetermined. In our case, due to limitations on the amount of memory space, the system was underdetermined (W had fewer rows than columns). The problem was therefore reformulated as an optimization problem,
Here, 
denotes an estimate of the intensity vector
, 
is a LaGrange multiplier, and the matrix Q is chosen
in order to control the behavior of the solution. Setting the
derivative of the above to zero, we obtain the equation
For the purposes of this project, it was found that choosing Q as the identity matrix produced satisfactory results, although other choices were also investigated (see the section on Wiener filtering). By choosing Q as the identity matrix, we obtain the solution to (16) that has the smallest norm, i.e., a reasonably smooth solution.
Given the size of the system matrix 
(e.g., 
in the case of a 
image and as used in this study), it
is not feasible to use methods such as Gauss elimination or LU
decomposition to solve this system. This is due both to the large
number of operations required to invert or factor such a matrix, and to the fact that the inverted or factored matrix would be more
dense, i.e., utilize more memory than the original matrix. In fact, it
is desirable to use a method which does not even require that the
matrix product be computed explicitly. For this project, we
chose to use the conjugate gradient method [5,6], which
is an iterative method frequently used in solving large systems of
equations arising, e.g., in the context of numerical solution of
partial differential equations. The conjugate gradient method only
requires that one is able to compute the result of multiplying
by a vector. Thus, 
may be computed as 
, which avoids computing (and storing) the denser . Computation times on the Sparc Ultras on Owlnet were generally
less than 10 minutes (CPU time) for a reconstruction
from 120 projections.
