On Diffusion Tensor Estimation

In this paper we propose a formal formulation for the estimation of Diffusion Tensors in the space of symmetric positive semidefinite (PSD) tensors. Traditionally, diffusion tensor model estimation has been carried out imposing tensor symmetry without constraints for negative eigenvalues. When diffusion weighted data does not follow the diffusion model, due to noise or signal drop, negative eigenvalues may arise. An estimation method that accounts for the positive definiteness is desirable to respect the underlying principle of diffusion. This paper proposes such an estimation method and provides a theoretical interpretation of the result. A closed-form solution is derived that is the optimal data-fit in the matrix 2-norm sense, removing the need for optimization-based tensor estimation.


I. INTRODUCTION
Diffusion tensors allow for the modeling and measurement of water diffusion inside the brain. They are used as indirect measures of brain connectivity, aligning with white matter fiber directions. A diffusion tensor is usually fitted to the diffusion weighted images (DWI) by solving the Stejskal-Tanner equation: (1) where S k are the measured diffusion weighted images in the gradient directions, g k , S 0 is the baseline image and D is the diffusion tensor. D is a symmetric positive definite matrix.
The classic method for fitting the tensor to the diffusion data is by means of linear regression, such that the following sum of squared differences error is minimized: (2) where S is the space of symmetric tensors. The minimization in the space of symmetric tensors is imposed by solving only for the 6 degrees of freedom of the symmetric D. However, there is no constraint that guarantees the resulting diffusion tensor to be PSD, i.e, the regression fitting is not computed in the manifold of PSD tensors.
Although indefinite diffusion tensors are non-physical, they commonly arise from diffusion tensor estimation in locations where the DWI signal is corrupted by noise and the degree of anisotropy is high [1], [2]. Indefiniteness is equivalent to σ(D) (the spectrum of D) possessing eigen-values of different signs (For D to be positive semidefinite σ(D) ≥ 0). Several authors have proposed methods to restrict the estimated tensor to the PSD manifold. To estimate a positive semidefinite diffusion tensor one can • project the tensor back onto the face of the cone of positive semidefinite tensors (at every iteration step in an iterative scheme), • assure positive semidefiniteness during estimation by enforcing the tensor flow to stay within the manifold of positive semidefinite tensors, or • assure by construction that the estimation can never leave the manifold of positive semidefinite tensors.
Projecting back onto the face of the cone can for example be accomplished by setting negative eigenvalues to zero (optimal in the Frobenius norm sense) or to their absolute values, which may not be the best choices [3]. The second approach has been investigated by Tschumperle and Deriche [4], who proposed a method for tensor estimation and regularization that ensures the positive-definiteness of the estimated diffusion tensor by forcing it to stay on the manifold of PSD tensors. (This is accomplished by a Riccati type equation and appropriate numerics that guarantee to preserve the positive-definiteness property.) Wang et al. [5] and Koay et al. [2] use Cholesky factorizations to ensure the positive-definiteness of the estimated tensor by construction. The methodology proposed in this paper relies on a symmetric square-root factorization of the diffusion tensor. In particular, • a gradient descent scheme is proposed that by construction estimates the diffusion tensor in the manifold of PSD tensors, • a formal interpretation of the least squares PSD tensor fitting in the 2-norm sense is presented, • and a closed-form solution for the 2-norm optimal data-fit is derived, allowing for a simple and computationally efficient estimation of positive definite tensors best matching the measured data.
Section II reviews and interprets projection based approaches for suppression of negative eigenvalues that are of relevance to the method proposed in this paper. Section III describes the gradient descent scheme based on the matrix square root decomposition, allowing for the estimation of guaranteed symmetric positive semidefinite tensors best fitting the measured data. Section IV presents analytic solutions to the symmetric PSD data-fitting problem in the case of six rotationally invariant gradients and theoretically analyzes the method. Results on real and synthetic datasets to illustrate our findings are given in Section V.

II. PSD TENSOR ESTIMATION BY PROJECTION
Finding the symmetric PSD tensor best fitting the measured data is a constrained optimization problem. The space of symmetric PSD tensors is non-linear; its faces (i.e., the boundary between the space of PSD and indefinite tensors), are given by the degenerate tensors, whose spectra contain at least one zero eigenvalue [6]. Those faces are known as the PSD cone faces. Unconstrained tensor estimation may yield tensors outside the PSD cone for noisy or very anisotropic data. To make these indefinite tensors positive semidefinite a (in some sense/norm optimal) tensor inside the PSD cone needs to be found. The optimal PSD tensor to an indefinite tensor will lie on the cone face, thus negative eigenvalues are set to zero.

A. Frobenius norm solution
Setting the negative eigenvalues of D to zero while keeping the remaining eigenvalues fixed, yields the Frobenius norm optimal solution [7], i.e., where .

B. Two-norm solution
The closest PSD tensor in the matrix 2-norm sense does not need to be unique [7]. Interestingly, the Frobenius norm solution is also one of the minimizers of the 2-norm, since D is symmetric by definition [7], i.e., where . To obtain a unique 2-norm optimal solution the optimization problem needs to be further constrained. One possibility is to search for a solution that best fits the measured data. This is the focus of Section III.

III. LEAST-SQUARES PSD CONSTRAINED TENSOR ESTIMATION
Solving Eq. (2) without restricting the space of admissible tensors S does not guarantee D to be positive semidefinite. A proper solution implies a problem formulation that assures a tensor remaining in the PSD space. For a given PSD matrix, A, there always exists a unique symmetric positive semidefinite matrix X such that A = XX, where X is called the principal square root of A [8]. In particular, this implies that A can always be decomposed into the square of a symmetric matrix (not necessarily positive semidefinite) X. Thus, by estimating the symmetric matrix X, A is assured to be symmetric positive semidefinite.
The symmetric PSD tensor D psd that optimally fits the measured data (in the sum of squared error sense) is given by (3) where D psd = XX. Estimating X instead of D psd directly renders the estimation problem quadratic (from linear) and does not generally allow for an analytic solution, but requires numerical optimization (for example by gradient descent or one of its flavors).
The gradient of the error E with respect to X is (4) where where X i is the i-th column vector of X and superscripts denote coordinates, i.e., . The gradient descent is then (5) Note, that the gradient (4) is symmetric. Thus D psd = XX will always be symmetric PSD. The induced gradient descent on D psd "lives" in PSD space.
The following Section shows that there is an analytic solution to the least-square problem given by Eq. (3) in the case where the number of gradients is equal to the number of unknowns (6 in the 3D case). The analytic solution consists of a linear fitting in the space of symmetric tensors followed by an eigenvalue correction.

IV. ANALYTIC SOLUTION
This Section presents analytic formulas to compute the 2-norm optimal symmetric PSD tensor best matching the measurement data when the number of unknowns is equal to the number of gradients or the symmetric tensor estimation matches the measurement data perfectly.
Minimizing Eq. (3), yields the 2-norm optimal solution that best matches the data 1 :
where the gradients g k are rotationally invariant. Then

A. 2D Example
To assess the behavior of the different estimation strategies, it is instructive to look at the 2 × 2 tensor case, since it allows for visualization in three dimensions. Given a tensor the faces of the cone are given by ac -b 2 = 0. Positive definite tensors are contained to the halfspace ac -b 2 > 0. Fig. 1 shows the surface traced out by the PSD cone faces. To evaluate our method, an indefinite tensor with eigenvalues σ(D) = {2,-1.5} and random eigenvectors was generated (diamond point). Diffusion weighted values have been obtained according to the model given by Eq. (1). Figs. 1a and 1b show the projections into the PSD cone achieved by the proposed method and the Frobenius norm approach as well as the the iterative solution obtained through the gradient descent given by Eq. 5. The gradient descent optimization is initialized by using a random tensor inside the cone. The 2-norm of the difference between the indefinite tensor D and the tensors of the PSD cone faces is shown as a colormap on the cone surface. One can see how the proposed solution lies on the minima line corresponding to the 2-norm solutions. While the Frobenius norm lies on the crest of minima values, the proposed solution is the only one that also fits the data best. Fig. 1b show the projection of the previous result in the a -b plane. The continuous and dashed contours represent the isolines on the cone of the 2-norm and the fitting error E respectively. It is clear to see how the proposed approach minimizes both metrics at the same time.

B. Real Case
This Section presents results of the proposed estimation method as applied to a real human brain dataset. The data corresponds to an LSDI (Line Scan Diffusion Imaging) acquisition acquired at Brigham's and Women's Hospital on a 1.5T GE scanner. A 6 gradient direction protocol with 2 baseline images was employed. The gradients are uniformly spaced fulfilling the conditions given by Eq. (8). Tensors have been estimated using the standard and the constrained Least Squares methods. Standard tensor estimation resulted in indefinite tensors in several isolated areas inside highly anisotropic structures like the corpus callosum as shown in Fig. 2. The proposed method leads to tensor estimations more consistent with the underlying data, as shown in Fig. 2 (the planar measure [9], c p = (λ 2 -λ 3 )/λ 1 , is more consistent along the corpus callosum for the proposed method).

VI. CONCLUSIONS
This paper investigated the estimation of symmetric positive semidefinite tensors. A gradient descent scheme for tensor estimation working within the manifold of symmetric positive semidefinite tensors was proposed. Further, an analytic solution was given for the special case of six rotationally invariant gradient directions in three dimensions and its optimality in the 2-norm sense was demonstrated. The analytic solution removes the need for iterative solution strategies and is computationally very efficient. and it follows that the global minimum is attained at for the optimal λ 1 , λ 2 , . Thus the solution is a minimizer in the 2-norm. However, the global minimum is only valid, if it results in a positive semidefinite tensor D psd , i.e., if and . Define Then, if and , the global minimum can be attained and , .
Assume, , then and is known. Then , (9) From Eq. (9)    Coronal slice of a real human brain. Indefinite tensors appear in the corpus callosum and yield higher than expected planar measures (see circles). The proposed method leads to a visibly smoother estimation of the planar measure within the highly anisotropic corpus callosum.