SVD for imaging systems with discrete rotational symmetry

The singular value decomposition (SVD) of an imaging system is a computationally intensive calculation for tomographic imaging systems due to the large dimensionality of the system matrix. The computation often involves memory and storage requirements beyond those available to most end users. We have developed a method that reduces the dimension of the SVD problem towards the goal of making the calculation tractable for a standard desktop computer. In the presence of discrete rotational symmetry we show that the dimension of the SVD computation can be reduced by a factor equal to the number of collection angles for the tomographic system. In this paper we present the mathematical theory for our method, validate that our method produces the same results as standard SVD analysis, and finally apply our technique to the sensitivity matrix for a clinical CT system. The ability to compute the full singular value spectra and singular vectors could augment future work in system characterization, image-quality assessment and reconstruction techniques for tomographic imaging systems.


Introduction
Medical imaging systems can be characterized by a linear system operator that maps a continuous object distribution through the system onto an image plane. The singular value decomposition (SVD) of provides an orthonormal set of basis functions that have numerous applications toward system characterization and image-quality assessment. For example, SVD data can be used to compute the null functions for an imaging system [1,2]. There have also been numerous publications that discuss the use of SVD in reconstruction methods [3][4][5][6]. Additionally there is ongoing work investigating the use of singular vectors as channels for mathematical observer models [7,8].
The forward model for an imaging system is often approximated by a sensitivity matrix H that maps voxels in object space to detector pixels in image space. To minimize the error in this approximation, it is desirable to have a large number of voxels to sample object space causing the column dimension N of the sensitivity matrix to be large. The row dimension of the sensitivity matrix is especially large in tomographic systems where M is the product of the number detector elements P and the number of collection angles J. The consequent large dimensions of H for tomographic medical imaging systems are problematic when using conventional SVD algorithms because they require adequate computer memory to compute and store HH † or H † H.
In this paper, we present a method to reduce the dimension of the SVD computation. We show that for tomographic systems with discrete rotational symmetry the dimension of the problem can be reduced by a factor equal to the number of collection angles. The theory and mathematics of our method are presented in Sections 2 -5. In Section 6, we present results from a proof-of-concept experiment. These results verify that our reduced dimension algorithm produced the same data set as conventional SVD analysis. Finally, we show the singular value spectra and singular vector data that we obtained by applying the reduced dimension SVD to the sensitivity matrix for a simulated 3rd generation cone beam x-ray CT system.

The system operator and the symmetry operators
The symmetry operators we will be considering arise from rotations about an axis of the imaging system. If an object is described by a function f (r) of the three dimensional location vector r, and R j is the matrix that describes a rotation by the angle about the symmetry axis, then an operator acting in object space is given by . The symmetry group in this case is Z J , the cyclic group of order J. This group is represented in by the set of matrices {I,R 1 ,…,R J−1 } which satisfy R j R k = R j+k and . These matrices can all be expressed in terms of R 1 by . Similarly this cyclic group is represented in object space by the set of operators { , ,…, }, which satisfy the same multiplication rules as the rotation matrices.
Components of the data vector will be indexed by a vector m given by m = (p, j) = (p x , p y , j). The vector p identifies the location of a detector element in a given detector array, and the index j identifies the detector array. These J identical detector arrays are assumed to be spaced around the object at angles and to have identical apertures in front of them. If h p (r) is the detector sensitivity function for the detector located at position p on the detector with j = 0, then the detector sensitivity function for detector m is given by . The mean data vector has components given by (1) which we will also write as . If the number of detector elements in one detector array is P, then is a concatenation of J P-dimensional vectors : (2) Thus is an operator from some function space U to data space . We represent this symbolically by . The decomposition of given above corresponds to the decomposition of V into an orthogonal direct sum of subspaces V j , all of which are Pdimensional: (3) A different decomposition of V into subspaces arises from symmetry considerations and will be described below.
The components of the subvectors can be thought of as the result of an operator acting on the object function: (4) We will consider to be an operator that takes an object function and gives a vector in the subspace V j of the full data space In other words, the result of applying this operator is the vector (5) Physically, the operator would be obtained by covering all of the apertures except the j th one. We will assume that the range of is all of V j . We may then write the action of the system operator on an object function as (6) with and .
Corresponding to the rotations of the object there are rotations of the data given by matrices S j . The matrix S 1 is described by the equation (7) and the matrices S j satisfy . We may assume that the detectors are numbered in the clockwise direction so that . This equation says that we can get the projection on the j th detector array by rotating the object by an angle of , getting the projection on the 0 th detector array, and then moving this P-dimensional data vector back to the correct location in the JP-dimensional vector .
The operator satisfies the symmetry condition . In mathematical terms we say the operator intertwines the group representation { , ,…, } with the representation {I, S 1 ,…,S J−1 }. Intuitively the symmetry condition tells us that rotating the object around the system axis gives the same result as rotating the ring of detectors through the same angle in the opposite direction. This condition places restrictions on the form of the system operator. These restrictions can be used to reduce the dimensionality of the SVD calculation for by at least a factor of J. Before we show this, we will introduce the inner products that we will use for object space and data space.

Inner products and adjoint operators
The SVD of any operator involves the notion of an adjoint operator. To define the adjoint of an operator, we need an inner product in the object space and in the data space. The inner product in object space will be defined by (8) where S is a support region for all object functions and is invariant under the rotations in the symmetry group. For example, S could be a circular cylinder whose axis is the axis of the rotations in the group. In data space, we have the standard inner product (9) We will have occasion to make use of complex object functions and data vectors, which is why the complex conjugate operation appears in these inner products.
We define adjoint operators and by the equations and . The symmetry operators are unitary: and . The first equality in each case is the unitary property, while the second follows from the definitions of the symmetry operators.
We can also define the adjoint of the system operator, , by the equation . This operator is a map from data space to object space: . In terms of the detector sensitivity functions, this adjoint operator is given by (10) This operator is also called the backprojection operator for the imaging system described by . We can use the the system operator given above to write (11) Notice that we used the unitary property of the symmetry operators here. The adjoint operator also maps data space to object space and its null space is the orthogonal complement to V j . Symbolically, we write and nullspace . Next, we will consider a different decomposition of data space that arises from the symmetry group of the system.

Projection operators
In order to reduce the dimensionality of the SVD problem for , we need to decompose object functions into components that are eigenfunctions of the operators and data vectors into components that are eigenvectors of the S j matrices. This process is essentially a Fourier decomposition with respect to the group of rotations. To this end define operators and P k , for k = 0,1,…,J−1 via the equations [9] (12) (13) These operators satisfy the idempotent property: and P k P k = P k . They are also Hermitian: and . These two properties imply that these operators are orthogonal projections onto their ranges. If we apply a symmetry operator after a projection operator the result is (14) ( 15) This shows that the range of each projection operator consists of functions (or vectors) that are simultaneous eigenfunctions (or eigenvectors) of all of the symmetry operators. Further more, we have the decompositions The components in each of these decompositions are orthogonal to each other. If is the range of , and is the range of P j , then we have the corresponding orthogonal decompositions of object space and data space This suggests that we define the operators by , where the second equality here follows from the symmerty property of . These operators are maps from object space to data space with the following ranges and null spaces: and . The system operator can be expressed as a sum of these operators: (20) In this case, the decomposition of is in block diagonal form in the sense that maps to and vanishes on for j ≠ k.
The adjoint of maps data space to object space and has the following properties: and . The adjoint system operator also has block diagonal from: Then we find that the JP×JP matrix has block diagonal form as well: Each block in this decomposition maps to and vanishes on for k ≠ j. This fact, together with the orthogonality of the decompositions of object and data space, implies that the SVD problem for the system operator reduces to finding the eigenvalues and eigenvectors for the J matrices . These matrices map a P-dimensional space to itself and therefore the corresponding eigenvalue problems reduce to finding the eigenvalues for P×P matrices. The original JP×JP eigenvalue problem for the SVD of has now been reduced to a set of J P×P eigenvalue problems, which is more computationally tractable, especially for large values of J.

Reduction of dimension for SVD
We now want to find the eigenvalues and eigenvectors of the linear operator . This operator is given by , where we have used the symmetry property of and the orthogonal projection properties of P k . We may write this last form out as a triple sum: Now we let q = p−l and re-index the sum over l to get (24) Next we let r = j+p−q and re-index the sum over j to arrive at (25) By separating the exponential factor, we can separate the sums. The end result is . We know that any eigenvector of this operator must lie in , so we will write such a vector as P k g in the eigenvalue equation: . Since any vector in may be reproduced from its values on the first detector, we may assume that g has the following form Now we look at the eigenvalue equation and notice that, since vanishes on we must have . In this equation, J 2 is replaced with J because P k has a 1/J factor. These last two vectors are equal iff . This last equation can be regarded as a P×P linear system since the operator maps V 0 to V 0 .
The result can also be formulated in terms of a lower dimensional SVD problem. We start with and write g 0 = Zg so that Z is a P×M matrix. With g given as above, we also have g = Z † g 0 . Our eigenvalue equation is now . Multiplying both sides by Z, we get . Using the orthogonal projection properties of , we can write this as . Therefore, the eigenvalue problem for each k is equivalent to the SVD of the operator . Finally, note that so that we could also do the SVD for these last two operators.
To make this calculation explicit, we define the operator by . Then we have the SVD equations for this operator in the data space for a single detector: . This is an eigenvector equation for a P×P matrix. The eigenvalues are real and non-negative. To get the singular vectors in the full data space, we use: . The factor is there to preserve normalization. The singular functions f kl in object space are then given by: . This expression can be expanded in terms of the single-detector system operator by using . This expression show that we can generate the singular object function by backprojecting a single detector singular vector with and then applying the operator . This operator rotates the function through each of the J angles, multiplies each rotated function by a phase factor determined by k, and then sums. The symmetry properties of the resulting singular functions are determined by k.

Examples
A proof-of-concept experiment was conducted to compare the singular data from the reduced dimension SVD to the singular data from a standard dimension SVD computation. Each technique was used to calculate the SVD of a simulated 2-D parallel beam source x-ray CT system. In the absence of scattering, the system operator for an x-ray CT system can be linearized by making the approximations discussed in Barrett and Myers [2]. The discrete representation of this linear operator is the sensitivity matrix for the modeled system. The sensitivity matrix, H, for the simulated x-ray CT system was computed using a discreteto-discrete forward model to map from object space U to image space V. The object space was discretized into N = 256×256 voxels with an imposed circular field of view. The image space was measured with a P = 32 detector element camera collecting at J = 90 equally spaced angles over 360 degrees of rotation about the center of the field of view. As discussed in Section 5, the reduced dimension SVD therefore led to J SVD computations of the P×N matrices . In comparison, the standard dimension SVD consisted of a single SVD computation of the JP×N matrix H. In both cases, the algorithm used to compute singular data for either the reduced dimension or the standard dimension matrix was the SVD function in the MATLAB software package [10]. Computations for the proofof-concept experiment were performed on an iMac computer equipped with a 3.2 GHz Intel Core i3 processor and 8GB of DDR3 memory. Calculation of all singular values and associated singular vectors up to the rank of H on this machine required 15.86 minutes using the reduced dimension SVD algorithm and 473.55 minutes using the standard dimension SVD algorithm. Our reduced dimension algorithm was therefore approximately 30 times faster than the standard dimension algorithm at computing the full SVD of the simulated sensitivity matrix. The speed of our algorithm could be further improved by trivially distributing the SVD computation of the A k matrices onto the cores of a multi-core machine such as a Mac Pro computer. Figure 1 shows the largest 16 singular values computed using the reduced dimension SVD and the standard dimension SVD. A slice through the central plane of the object space vectors for the 10 largest singular values is presented in Fig. 2.
In Fig. 1 we can see that the some singular vectors occur in pairs with the same singular value (doublets), while others occur by themselves (singlets). The existence of this type of pattern would have been predicted from an analysis that made use of the full symmetry group of the system, which would include reflections as well as rotations [11][12][13]. By replacing the cyclic group of rotations used in this work with the dihedral group of rotations and reflections, we can show that the spaces spanned by the singular vectors that share a singular value are predicted to be either one-dimensional or two-dimensional. This prediction follows from the fact that the irreducible representations of the dihedral group are either one-dimensional or two-dimensional [9].
If we look at the standard dimension SVD data in column (a) of Fig. 2, we can see that singular vectors in the doublets are related to each other by a reflection through an axis of symmetry of the system. On the other hand, singular vectors in the singlets are invariant under reflections. This structure is also apparent in the data from the reduced dimension SVD data in columns (b) and (c) of Fig. 2. Since the operator HH † is a real matrix, the real and imaginary parts of the singular vectors that the symmetry analysis produced are separately real eigenvectors of this matrix with the same eigenvalue. Since doublet eigenvectors share an eigenvalue, each of the four real eigenvectors corresponding to a given doublet eigenvalue in columns (b) and (c) is only required to be a linear combination of the two eigenvectors for that doublet in column (a). This extra degree of freedom for the doublets explains why the vectors in the column (a) do not exactly match up with any of the vectors in columns (b) and (c) for the doublet eigenvalues.
In order to check the outputs from the reduced dimension algorithm we tested the orthonormality of the computed singular vectors. Orthonormality of singular vectors u n can be expressed mathematically in terms of the inner product and the Kronecker delta as [2] (26) Shown in Table 1 and Table 2 are the magnitudes of the inner products for selected singular vectors that were computed using the reduced dimension SVD algorithm. The magnitudes of the inner products of singular vectors with the same index are unity while the magnitudes of the inner products of singular vectors with differing indices are close to zero. The error in the inner products is due to the estimations we made when simulating H for the proof-ofconcept system. In our forward model, we used subvoxels and monte carlo techniques when quantifying the values of H at each collection angle. Consequently, our forward model does not strictly exhibit the property of discrete rotational symmetry discussed in Section 2. To further assess how this modeling error manifested itself in the SVD data, we compared singular vectors produced by the reduced dimension and standard dimension techniques. For this analysis we chose singlets u 1 , u 6 and u 15 . We subtracted the singlet output by the reduced dimension SVD from the singlet output by the standard SVD resulting in an image of the error for the singlet. The results of this subtraction are shown in Fig. 3. The structure of these images shows that that there are contributions from additional singular vectors within each singlet that are a consequence of not preserving discrete rotational symmetry in our forward model. Note that the magnitude of the inner product between two singular vectors of differing indices that came from a reduced dimension SVD computation of the same A k matrix (e.g. u 1 and u 6 ) are in fact zero to within the numerical precision of the computer.
The reduced dimension SVD was next used to calculate singular data for the sensitivity matrix of a simulated clinical x-ray CT system. The specifications for the simulated system were chosen to roughly approximate a 3rd generation Siemens Sensation 64 scanner (Siemens Healthcare) [14]. In order to shorten the computation time the number of slices for our modeled system was changed from 64 slices to 12 slices. The model consisted of a cone beam source with source to center of rotation distance of R F = 600 mm and a center of rotation to center of detector array distance of R D = 450 mm. The slice thickness for each of the 12 detector column was 0.5 mm and each detector row had 672 channels of 0.9 mm width. This detector geometry was modeled as a rectangular array with P = 8064 detector elements. The detector array collected data at J = 180 equally spaced angles over 360 degrees. Object space was centered at the center of rotation and discretized using N = 512 × 512 × 5 square voxels of dimension 0.6741 mm. This choice of voxel parameters ensured that a ray traveling from the source through the edge of object space was observed by the detector array. A circular field of view was also imposed within each 512 × 512 slice through object space to preserve rotational symmetry. An illustration of the geometry for the modeled system is shown in Fig. 4.
Using a geometric forward projector, the sensitivity matrix for the simulated Siemens scanner was calculated. Applying the reduced dimension algorithm resulted in J SVD computations of P × N matrices. The singular data for a given P × N matrix was computed using a power methods algorithm [3,15]. Figure 5 shows the singular value spectra that was computed for the simulated Siemens x-ray CT system. Figure 6 shows the central 512 × 512 slice through the object space singular vector associated with the singular value spectra for singular value indices 1 -10. Similarly, Fig. 7 and Fig. 8 show the central slice through the singular vectors associated with indices 88 -92 and 175 -179 respectively.

Conclusion
In this paper, we have shown that in the presence of discrete rotational symmetry, computation of the SVD for the imaging operator can be reduced by a factor equal to the number of collection angles for the system geometry. The results from our proof-of-concept experiment confirm that the reduced dimension SVD algorithm results in a decomposition that is analogous to that produced by a standard dimension SVD algorithm. However, data from the reduced SVD algorithm will be corrupted by approximations and estimates in the forward model that break the discrete rotational symmetry of the system operator. The proof-of-concept experiment also showed that the reduced dimension algorithm can significantly shorten the time required to compute an SVD. Using the same computing resources, the reduced SVD algorithm was able to calculate the singular value spectra and associated singular vectors up to the rank of our tested system matrix in 15.86 minutes compared to 463.55 minutes for a standard dimension SVD computation which is a factor of 30 times faster.
Our algorithm is especially useful for tomographic medical imaging systems in which the system matrix has large dimensions resulting in memory issues for a standard desktop computer running a SVD computation. As we have demonstrated, we were able to apply the reduced dimension SVD technique to the simulated system matrix of a 3rd generation clinical x-ray CT system resulting in the singular values and singular vectors associated with 180 values from the singular value spectra. If we were to run our algorithm for a longer time, we could compute more terms from the spectra yielding a data set that could be used to quantify image-quality assessment metrics such as the null space and measurement space for the imaging system [2]. Our reduced dimension SVD algorithm can therefore be used for tomographic medical imaging systems to output the higher order terms of the singular value spectra that would require extensive computing resources using standard SVD. Plots showing the largest 16 singular values computed using the standard dimension SVD and the reduced dimension SVD. The spectra output by each SVD technique is equivalent.  The difference between singular vectors output by the Standard SVD and Reduced SVD algorithms. An illustration of the geometry of a simulated clinical x-ray CT system for which we computed the SVD using the reduced dimension SVD algorithm. This system specifications were chosen to roughly approximate the Siemens Sensation 64 scanner (Media 1, Media 2, Media 3, Media 4, Media 5, Media 6, Media 7, Media 8, Media 9, and Media 10).

Fig. 5.
A plot showing the singular value spectra for the simulated Siemens x-ray CT system. These data were computed using the reduced dimension SVD analysis. Images of slices through the central plane of the object space vectors for the modeled Siemens x-ray CT systems associated with singular value indices 1-10. Images of slices through the central plane of the object space vectors for the modeled Siemens x-ray CT systems associated with singular value indices 88-92.   Table 2 The magnitudes of the inner products of singular vectors that were computed using the reduced dimension SVD. Shown in this table are for m = {9,10,…,16} and n = {1,2,…,16}