Non-Euclidean statistics for covariance matrices, with applications to diffusion tensor imaging
aa r X i v : . [ s t a t . A P ] O c t The Annals of Applied Statistics (cid:13)
Institute of Mathematical Statistics, 2009
NON-EUCLIDEAN STATISTICS FOR COVARIANCE MATRICES,WITH APPLICATIONS TO DIFFUSION TENSOR IMAGING By Ian L. Dryden, Alexey Koloydenko and Diwei Zhou
University of South Carolina, Royal Holloway, University of London andUniversity of Nottingham
The statistical analysis of covariance matrix data is consideredand, in particular, methodology is discussed which takes into accountthe non-Euclidean nature of the space of positive semi-definite sym-metric matrices. The main motivation for the work is the analysis ofdiffusion tensors in medical image analysis. The primary focus is onestimation of a mean covariance matrix and, in particular, on the useof Procrustes size-and-shape space. Comparisons are made with otherestimation techniques, including using the matrix logarithm, matrixsquare root and Cholesky decomposition. Applications to diffusiontensor imaging are considered and, in particular, a new measure offractional anisotropy called Procrustes Anisotropy is discussed.
1. Introduction.
The statistical analysis of covariance matrices occursin many important applications, for example, in diffusion tensor imaging[Alexander (2005); Schwartzman, Dougherty and Taylor (2008)] or longitu-dinal data analysis [Daniels and Pourahmadi (2002)]. We wish to considerthe situation where the data at hand are sample covariance matrices, andwe wish to estimate the population covariance matrix and carry out statisti-cal inference. An example application is in diffusion tensor imaging where adiffusion tensor is a covariance matrix related to the molecular displacementat a particular voxel in the brain, as described in Section 2.If a sample of covariance matrices is available, we wish to estimate anaverage covariance matrix, or we may wish to interpolate in space betweentwo or more estimated covariance matrices, or we may wish to carry outtests for equality of mean covariance matrices in different groups.
Received May 2008; revised March 2009. Supported by a Leverhulme Research Fellowship and a Marie Curie Research Trainingaward.
Key words and phrases.
Anisotropy, Cholesky, geodesic, matrix logarithm, principalcomponents, Procrustes, Riemannian, shape, size, Wishart.
This is an electronic reprint of the original article published by theInstitute of Mathematical Statistics in
The Annals of Applied Statistics ,2009, Vol. 3, No. 3, 1102–1123. This reprint differs from the original in paginationand typographic detail. 1
I. L. DRYDEN, A. KOLOYDENKO AND D. ZHOU
The usual approach to estimating mean covariance matrices in Statistics isto assume a scaled Wishart distribution for the data, and then the maximumlikelihood estimator (m.l.e.) of the population covariance matrix is the arith-metic mean of the sample covariance matrices. The estimator can be formu-lated as a least squares estimator using Euclidean distance. However, sincethe space of positive semi-definite symmetric matrices is a non-Euclideanspace, it is more natural to use alternative distances. In Section 3 we definewhat is meant by a mean covariance matrix in a non-Euclidean space, us-ing the Fr´echet mean. We then review some recently proposed techniquesbased on matrix logarithms and also consider estimators based on matrixdecompositions, such as the Cholesky decomposition and the matrix squareroot.In Section 4 we introduce an alternative approach to the statistical anal-ysis of covariance matrices using the Kendall’s (1989) size-and-shape space.Distances, minimal geodesics, sample Fr´echet means, tangent spaces andpractical estimators based on Procrustes analysis are all discussed. We in-vestigate properties of the estimators, including consistency.In Section 5 we compare the various choices of metrics and their prop-erties. We investigate measures of anisotropy and discuss the deficient rankcase in particular. We consider the motivating applications in Section 6where the analysis of diffusion tensor images and a simulation study areinvestigated. Finally, we conclude with a brief discussion.
2. Diffusion tensor imaging.
In medical image analysis a particular typeof covariance matrix arises in diffusion weighted imaging called a diffusiontensor. The diffusion tensor is a 3 × x ∈ R has probabilitydensity function f ( x ) = 1(2 π ) / | Σ | / exp (cid:18) − x T Σ − x (cid:19) . The convention is to call D = Σ / g , g , . . . , g m ∈ RP withscanner gradient parameter b , where RP is the real projective space ofaxial directions (with g j ≡ − g j , k g j k = 1). The data at a voxel consist ofsignals ( Z , Z , . . . , Z m ) which are related to the Fourier transform of the ON-EUCLIDEAN STATISTICS FOR COVARIANCE MATRICES Fig. 1.
Visualization of a diffusion tensor as an ellipsoid. The principal axis is alsodisplayed. displacements in axial direction g j ∈ RP , j = 1 , . . . , m , and the reading Z is obtained with no gradient ( b = 0). The Fourier transform in axial direction g ∈ RP of the multivariate Gaussian displacement density is given by F ( g ) = Z exp( i √ bg ) f ( x ) dx = exp( − bg T Dg ) , and the theoretical model for the signals is Z j = Z F ( g j ) = Z exp( − bg T j Dg j ) , j = 1 , . . . , m. There are a variety of methods available for estimating D from the data( Z , Z , . . . , Z m ) at each voxel [see Alexander (2005)], including least squaresregression and Bayesian estimation [e.g., Zhou et al. (2008)]. Noise modelsinclude log-Gaussian, Gaussian and, more recently, Rician noise [Wang et al.(2004); Fillard et al. (2007); Basu, Fletcher and Whitaker (2006)]. A commonmethod for visualizing a diffusion tensor is an ellipsoid with principal axesgiven by the eigenvectors of D , and lengths of axes proportional to √ λ i , i =1 , , . An example is given in Figure 1.If a sample of diffusion tensors is available, we may wish to estimate anaverage diffusion tensor matrix, investigate the structure of variability indiffusion tensors or interpolate at higher spatial resolution between two ormore estimated diffusion tensor matrices.In diffusion tensor imaging a strongly anisotropic diffusion tensor indicatesa strong direction of white matter fiber tracts, and plots of measures ofanisotropy are very useful to neurologists. A measure that is very commonlyused in diffusion tensor imaging is Fractional Anistropy,
F A = ( kk − k X i =1 ( λ i − ¯ λ ) (cid:30) k X i =1 λ i ) / , where 0 ≤ F A ≤ λ i are the eigenvalues of the diffusion tensor matrix.Note that F A ≈ λ ≫ λ i ≈ , i > F A = 0 for isotropy. In diffusion tensor imaging k = 3. I. L. DRYDEN, A. KOLOYDENKO AND D. ZHOU
Fig. 2.
An FA map from a slice in a human brain. Lighter values indicate higher FA.
In Figure 2 we see a plot of FA from an example healthy human brain. Wefocus on the small inset region in the box, and we would like to interpolatethe displayed image to a finer scale. We return to this application in Section6.3.
3. Covariance matrix estimation.
Euclidean distance.
Let us consider n sample covariance matrices(symmetric and positive semi-definite k × k matrices) S , . . . , S n which areour data (or sufficient statistics). We assume that the S i are independentand identically distributed (i.i.d.) from a distribution with mean covariancematrix Σ, although we shall elaborate more later in Section 3.2 about whatis meant by a “mean covariance matrix.” The main aim is to estimate Σ.More complicated modeling scenarios are also of interest, but for now wejust concentrate on estimating the mean covariance matrix Σ.The most common approach is to assume i.i.d. scaled Wishart distribu-tions for S i with E ( S i ) = Σ, and the m.l.e. for Σ is ˆΣ E = n P ni =1 S i . Thisestimator can also be obtained if using a least squares approach by minimiz-ing the sum of square Euclidean distances. The Euclidean distance betweentwo matrices is given by d E ( S , S ) = k S − S k = q trace { ( S − S ) T ( S − S ) } , (1)where k X k = q trace( X T X ) is the Euclidean norm (also known as the Frobe-nius norm). The least squares estimator is given byˆΣ E = arg inf Σ n X i =1 k S i − Σ k . However, the space of positive semi-definite symmetric matrices is a non-Euclidean space and other choices of distance are more natural. One par-ticular drawback with Euclidean distance is when extrapolating beyond the
ON-EUCLIDEAN STATISTICS FOR COVARIANCE MATRICES data, nonpositive semi-definite estimates can be obtained. There are otherdrawbacks when interpolating covariance matrices, as we shall see in ourapplications in Section 6.3.2. The Fr´echet mean.
When using a non-Euclidean distance d ( · ) wemust define what is meant by a “mean covariance matrix.” Consider a prob-ability distribution for a k × k covariance matrix S on a Riemannian metricspace with density f ( S ). The Fr´echet (1948) mean Σ is defined asΣ = arg inf Σ Z d ( S, Σ) f ( S ) dS, and is also known as the Karcher mean [Karcher (1977)]. The Fr´echet meanneed not be unique in general, although for many distributions it will be.Provided the distribution is supported only on the geodesic ball of radius r ,such that the geodesic ball of radius 2 r is regular [i.e., supremum of sectionalcurvatures is less than ( π/ (2 r )) ], then the Fr´echet mean Σ is unique [Le(1995)]. The support to ensure uniqueness can be very large. For example,for Euclidean spaces (with sectional curvature zero), or for non-Euclideanspaces with negative sectional curvature, the Fr´echet mean is always unique.If we have a sample S , . . . , S n of i.i.d. observations available, then thesample Fr´echet mean is calculated by findingˆΣ = arg inf Σ n X i =1 d ( S i , Σ) . Uniqueness of the sample Fr´echet mean can also be determined from theresult of Le (1995).3.3.
Non-Euclidean covariance estimators.
A recently derived approachto covariance matrix estimation is to use matrix logarithms. We write thelogarithm of a positive definite covariance matrix S as follows. Let S = U Λ U T be the usual spectral decomposition, with U ∈ O ( k ) an orthogonalmatrix and Λ diagonal with strictly positive entries. Let log Λ be a diagonalmatrix with logarithm of the diagonal elements of Λ on the diagonal. Thelogarithm of S is given by log S = U (log Λ) U T and likewise the exponentialof the matrix S is exp S = U (exp Λ) U T . Arsigny et al. (2007) propose theuse of the log-Euclidean distance, where Euclidean distance between thelogarithm of covariance matrices is used for statistical analysis, that is, d L ( S , S ) = k log( S ) − log( S ) k . (2)An estimator for the mean population covariance matrix using this ap-proach is given byˆΣ L = exp ( arg inf Σ n X i =1 k log S i − log Σ k ) = exp ( n n X i =1 log S i ) . I. L. DRYDEN, A. KOLOYDENKO AND D. ZHOU
Using this metric avoids extrapolation problems into matrices with nega-tive eigenvalues, but it cannot deal with positive semi-definite matrices ofdeficient rank.A further logarithm-based estimator uses a Riemannian metric in thespace of square symmetric positive definite matrices d R ( S , S ) = k log( S − / S S − / ) k . (3)The estimator (sample Fr´echet mean) is given byˆΣ R = arg inf Σ n X i =1 k log( S − / i Σ S − / i ) k , which has been explored by Pennec, Fillard and Ayache (2006), Moakher(2005), Schwartzman (2006), Lenglet, Rousson and Deriche (2006) andFletcher and Joshi (2007). The estimate can be obtained using a gradi-ent descent algorithm [e.g., see Pennec (1999); Pennec, Fillard and Ayache(2006)]. Note that this Riemannian metric space has negative sectional cur-vature and so the population and sample Fr´echet means are unique in thiscase.Alternatively, one can use a reparameterization of the covariance matrix,such as the Cholesky decomposition [Wang et al. (2004)], where S i = L i L T i and L i = chol( S i ) is lower triangular with positive diagonal entries. TheCholesky distance is given by d C ( S , S ) = k chol( S ) − chol( S ) k . (4)A least squares estimator can be obtained fromˆΣ C = ˆ∆ C ˆ∆ T C , where ˆ∆ C = arg inf ∆ ( n n X i =1 k L i − ∆ k ) = 1 n n X i =1 L i . An equivalent model-based approach would use an independent Gaussianperturbation model for the lower triangular part of L i , with mean given bythe lower triangular part of ∆ C , and so ˆ∆ C is the m.l.e. of ∆ C under thismodel. Hence, in this approach the averaging is carried out on a square roottype-scale, which would indeed be the case for k = 1 dimensional case wherethe estimate of variance would be the square of the mean of the samplestandard deviations.An alternative decomposition is the matrix square root where S / = U Λ / U T , which has not been used in this context before as far as we areaware. The distance is given by d H ( S , S ) = k S / − S / k . (5) ON-EUCLIDEAN STATISTICS FOR COVARIANCE MATRICES A least squares estimator can be obtained fromˆΣ H = ˆ∆ H ˆ∆ T H , where ˆ∆ H = arg inf ∆ ( n X i =1 k S / i − ∆ k ) = 1 n n X i =1 S / i . However, because L i RR T L T i = L i L T i for R ∈ O ( k ) , another new alterna-tive is to relax the lower triangular or square root parameterizations andmatch the initial decompositions closer in terms of Euclidean distance byoptimizing over rotations and reflections. This idea provides the rationalefor the main approaches in this paper.
4. Procrustes size-and-shape analysis.
Non-Euclidean size-and-shape metric.
The non-Euclidean size-and-shape metric between two k × k covariance matrices S and S is definedas d S ( S , S ) = inf R ∈ O ( k ) k L − L R k , (6)where L i is a decomposition of S i such that S i = L i L T i , i = 1 ,
2. For example,we could have the Cholesky decomposition L i = chol( S i ) , i = 1 , , which islower triangular with positive diagonal elements, or we could consider thematrix square root L = S / = U Λ / U T , where S = U Λ U T is the spectraldecomposition. Note that S = ( L R )( L R ) T for any R ∈ O ( k ), and so thedistance involves matching L optimally, in a least-squares sense, to L byrotation and reflection. Since S = LL T , then the decomposition is repre-sented by an equivalence class { LR : R ∈ O ( k ) } . For practical computationwe often need to choose a representative from this class, called an icon, andin our computations we shall choose the Cholesky decomposition.The Procrustes solution for matching L to L isˆ R = arg inf R ∈ O ( k ) k L − L R k (7) = U W T , where L T1 L = W Λ U T , U, W ∈ O ( k ) , and Λ is a diagonal matrix of positive singular values [e.g., see Mardia, Kentand Bibby (1979), page 416].This metric has been used previously in the analysis of point set configura-tions where invariance under translation, rotation and reflection is required.Size-and-shape spaces were introduced by Le (1988) and Kendall (1989) aspart of the pioneering work on the shape analysis of landmark data [cf.Kendall (1984)]. The detailed geometry of these spaces is given by Kendallet al. [(1999), pages 254–264], and, in particular, the size-and-shape space isa cone with a warped-product metric and has positive sectional curvature. I. L. DRYDEN, A. KOLOYDENKO AND D. ZHOU
Equation (6) is a Riemannian metric in the reflection size-and-shape spaceof ( k + 1)-points in k dimensions [Dryden and Mardia (1998), Chapter 8].In particular, d S ( · ) is the reflection size-and-shape distance between the( k + 1) × k configurations H T L and H T L , where H is the k × ( k + 1)Helmert sub-matrix [Dryden and Mardia (1998), page 34] which has j throw given by( h j , . . . , h j | {z } j times , − jh j , , . . . , | {z } k − j times ) , h j = −{ j ( j + 1) } − / , for j = 1 , . . . , k .Hence, the statistical analysis of covariance matrices can be consideredequivalent to the dual problem of analyzing reflection size-and-shapes.4.2. Minimal geodesic and tangent space.
Let us consider the minimalgeodesic path through the reflection size-and-shapes of L and L in thereflection size-and-shape space, where L i L T i = S i , i = 1 ,
2. Following an ar-gument similar to that for the minimal geodesics in shape spaces [Kendall etal. (1999)], this minimal geodesic can be isometrically expressed as L + tT, where T are the horizontal tangent co-ordinates of L with pole L . Kendallet al. [(1999), Section 11.2] discuss size-and-shape spaces without reflectioninvariance, however, the results with reflection invariance are similar, as re-flection does not change the local geometry.The horizontal tangent coordinates satisfy L T T = T L T1 [Kendall et al.(1999), page 258]. Explicitly, the horizontal tangent coordinates are givenby T = L ˆ R − L , ˆ R = inf R ∈ O ( k ) k L − L R k , where ˆ R is the Procrustes match of L onto L given in (7). So, the geodesicpath starting at L and ending at L is given by w L + w L ˆ R, where w + w = 1 , w i ≥ , i = 1 ,
2, and ˆ R is given in (7). Minimal geodesicsare useful in applications for interpolating between two covariance matrices,in regression modeling of a series of covariance matrices, and for extrapola-tion and prediction.Tangent spaces are very useful in practical applications, where one usesEuclidean distances in the tangent space as approximations to the non-Euclidean metrics in the size-and-shape space itself. Such constructions areuseful for approximate multivariate normal based inference, dimension re-duction using principal components analysis and large sample asymptoticdistributions. ON-EUCLIDEAN STATISTICS FOR COVARIANCE MATRICES Procrustes mean covariance matrix.
Let S , . . . , S n be a sample of n positive semi-definite covariance matrices each of size k × k from a dis-tribution with density f ( S ), and we work with the Procrustes metric (6) inorder to estimate the Fr´echet mean covariance matrix Σ. We assume that f ( S ) leads to a unique Fr´echet mean (see Section 3.2).The sample Fr´echet mean is calculated by findingˆΣ S = arg inf Σ n X i =1 d S ( S i , Σ) . In the dual size-and-shape formulation we can write this asˆΣ S = ˆ∆ S ˆ∆ T S , where ˆ∆ S = arg inf ∆ n X i =1 inf R i ∈ O ( k ) k H T L i R i − H T ∆ k . (8)The solution can be found using the Generalized Procrustes Algorithm[Gower (1975); Dryden and Mardia (1998), page 90], which is available inthe shapes library (written by the first author of this paper) in R [R De-velopment Core Team (2007)]. Note that if the data lie within a geodesicball of radius r such that the geodesic ball of radius 2 r is regular [Le (1995);Kendall (1990)], then the algorithm finds the global unique minimum solu-tion to (8). This condition can be checked for any dataset and, in practice,the algorithm works very well indeed.4.4. Tangent space inference.
If the variability in the data is not toolarge, then we can project the data into the tangent space and carry out theusual Euclidean based inference in that space.Consider a sample S , . . . , S n of covariance matrices with sample Fr´echetmean ˆΣ S and tangent space coordinates with pole ˆΣ S = ˆ∆ S ˆ∆ T S given by V i = ˆ∆ S − L i ˆ R i , where ˆ R i is the Procrustes rotation for matching L i to ˆ∆ S , i = 1 , . . . , n, and S i = L i L T i , i = 1 , S v = 1 n n X i =1 vec( V i ) vec( V i ) T , where vec is the vectorize operation. The principal component (PC) load-ings are given by ˆ γ j , j = 1 , . . . , p , the eigenvectors of S v corresponding tothe eigenvalues ˆ λ ≥ ˆ λ ≥ · · · ≥ ˆ λ p >
0, where p is the number of nonzeroeigenvalues. The PC score for the i th individual on PC j is given by s ij = ˆ γ T j vec( V i ) , i = 1 , . . . , n, j = 1 , . . . , p. I. L. DRYDEN, A. KOLOYDENKO AND D. ZHOU
In general, p = min( n − , k ( k + 1) / j th PC can beexamined by evaluatingΣ( c ) = ( ˆ∆ S + c vec − k (ˆ λ / j ˆ γ j ))( ˆ∆ S + c vec − k (ˆ λ / j ˆ γ j )) T for various c [often in the range c ∈ ( − , − k (vec( V )) = V for a k × k matrix V .Tangent space inference can proceed on the first p PC scores, or possiblyin lower dimensions if desired. For example, Hotelling’s T test can be carriedout to examine group differences, or regression models could be developedfor investigating the PC scores as responses versus various covariates. Weshall consider principal components analysis of covariance matrices in anapplication in Section 6.2.4.5. Consistency.
Le (1995, 2001) and Bhattacharya and Patrangenaru(2003, 2005) provide consistency results for Riemannian manifolds, whichcan be applied directly to our situation. Consider a distribution F on thespace of covariance matrices which has size-and-shape Fr´echet mean Σ S . Let S , . . . , S n be i.i.d. from F , such that they lie within a geodesic ball B r suchthat B r is regular. ThenˆΣ S P → Σ S , as n → ∞ , where Σ S is unique. In addition, we can derive a central limit theorem resultas in Bhattacharya and Patrangenaru (2005), where the tangent coordinateshave an approximate multivariate normal distribution for large n . Hence,confidence regions based on the bootstrap can be obtained, as in Amaral,Dryden and Wood (2007) and Bhattacharya and Patrangenaru (2003, 2005).4.6. Scale invariance.
In some applications it may be of interest to con-sider invariance over isotropic scaling of the covariance matrix. In this casewe could consider the representation of the covariance matrix using Kendall’sreflection shape space, with the shape metric given by the full Procrustesdistance d F ( S , S ) = inf R ∈ O ( k ) ,β> (cid:13)(cid:13)(cid:13)(cid:13) L k L k − βL R (cid:13)(cid:13)(cid:13)(cid:13) , (9)where S i = L i L T i , i = 1 ,
2, and β > S , . . . , S n , which is scaleinvariant and based on the full Procrustes mean shape (extrinsic mean), isˆΣ F = ˆ∆ F ˆ∆ T F , where ˆ∆ F = arg inf ∆ n X i =1 n inf R i ∈ O ( k ) , β i > k β i L i R i − ∆ k o , ON-EUCLIDEAN STATISTICS FOR COVARIANCE MATRICES Table 1
Notation and definitions of the distances and estimators
Name Notation Form Estimator Equation
Euclidean d E ( S , S ) k S − S k ˆΣ E (1)Log-Euclidean d L ( S , S ) k log( S ) − log( S ) k ˆΣ L (2)Riemannian d R ( S , S ) k log( S − / S S − / ) k ˆΣ R (3)Cholesky d C ( S , S ) k chol( S ) − chol( S ) k ˆΣ C (4)Root Euclidean d H ( S , S ) k S / − S / k ˆΣ H (5)Procrustes size-and-shape d S ( S , S ) inf R ∈ O ( k ) k L − L R k ˆΣ S (6)Full Procrustes shape d F ( S , S ) inf R ∈ O ( k ) ,β> k L k L k − βL R k ˆΣ F (9)Power Euclidean d A ( S , S ) α k S α − S α k ˆΣ A (10) and S i = L i L T i , i = 1 , . . . , n , and β i > shapes library in R. Tangent space inference can then proceed in an analogousmanner to that of Section 4.4.
5. Comparison of approaches.
Choice of metrics.
In applications there are several choices of dis-tances between covariance matrices that one could consider. For complete-ness we list the metrics and the estimators considered in this paper in Table1, and we discuss briefly some of their properties.Estimators ˆΣ E , ˆΣ C , ˆΣ H , ˆΣ L , ˆΣ A are straightforward to compute using arith-metic averages. The Procrustes based estimators ˆΣ S , ˆΣ F involve the use ofthe Generalized Procrustes Algorithm, which works very well in practice.The Riemannian metric estimator ˆΣ R uses a gradient descent algorithmwhich is guaranteed to converge [see Pennec (1999); Pennec, Fillard andAyache (2006)].All these distances except d C are invariant under simultaneous rotationand reflection of S and S , that is, the distances are unchanged by replacingboth S i by V S i V T , V ∈ O ( k ) , i = 1 ,
2. Metrics d L ( · ) , d R ( · ) , d F ( · ) are invari-ant under simultaneous scaling of S i , i = 1 ,
2, that is, replacing both S i by βS i , β >
0. Metric d R ( · ) is also affine invariant, that is, the distances areunchanged by replacing both S i by AS i A T , i = 1 , , where A is a general k × k full rank matrix. Metrics d L ( · ) , d R ( · ) have the property that d ( A, I k ) = d ( A − , I k ) , where I k is the k × k identity matrix.Metrics d L ( · ) , d R ( · ) , d F ( · ) are not valid for comparing rank deficient co-variance matrices. Finally, there are problems with extrapolation with met-ric d E ( · ): extrapolate too far and the matrices are no longer positive semi-definite. I. L. DRYDEN, A. KOLOYDENKO AND D. ZHOU
Anisotropy.
In some applications a measure of anisotropy of thecovariance matrix may be required, and in Section 2 we described the com-monly used FA measure. An alternative is to use the full Procrustes shapedistance to isotropy and we have PA = s kk − d F ( I k , S ) = s kk − R ∈ O ( k ) ,β ∈ R + (cid:13)(cid:13)(cid:13)(cid:13) I k √ k − β chol( S ) R (cid:13)(cid:13)(cid:13)(cid:13) , = ( kk − k X i =1 ( √ λ i − √ λ ) (cid:30) k X i =1 λ i ) / , where √ λ = k P √ λ i . Note that the maximal value of d F distance fromisotropy to the rank 1 covariance matrix is p ( k − /k, which follows from Le(1992). We include the scale factor when defining the Procrustes Anisotropy(PA), and so 0 ≤ PA ≤
1, with PA = 0 indicating isotropy, and PA ≈ d L or d R is the geodesic anisotropy GA = ( k X i =1 (log λ i − log λ ) ) / , where 0 ≤ GA < ∞ [Arsigny et al. (2007); Fillard et al. (2007); Fletcher andJoshi (2007)], which has been used in diffusion tensor analysis in medicalimaging with k = 3.5.3. Deficient rank case.
In some applications covariance matrices areclose to being deficient in rank. For example, when FA or PA are equalto 1, then the covariance matrix is of rank 1. The Procrustes metrics caneasily deal with deficient rank matrices, which is a strong advantage of theapproach. Indeed, Kendall’s (1984, 1989) original motivation for developinghis theory of shape was to investigate rank 1 configurations in the contextof detecting “flat” (collinear) triangles in archeology.The use of ˆΣ L and ˆΣ R has strong connections with the use of Bookstein’s(1986) hyperbolic shape space and Le and Small’s (1999) simplex shapespace, and such spaces cannot deal with deficient rank configurations.The use of the Cholesky decomposition has strong connections with Book-stein coordinates and Goodall–Mardia coordinates in shape analysis, whereone registers configurations on a common baseline [Bookstein (1986); Goodalland Mardia (1992)]. For small variability the baseline registration methodand Procrustes superimposition techniques are similar, and there is an ap-proximate linear relationship between the two [Kent (1994)]. In shape anal-ysis edge superimposition techniques can be very unreliable if the baselineis very small in length, which would correspond to very small variability in ON-EUCLIDEAN STATISTICS FOR COVARIANCE MATRICES Fig. 3.
Four different geodesic paths between the two tensors. The geodesic paths areobtained using d E ( · ) (1st row), d L ( · ) (2nd row), d C ( · ) (3rd row) and d S ( · ) (4th row). particular diagonal elements of the covariance matrix in the current context.Cholesky methods would be unreliable in such cases. Also, Bookstein coor-dinates induce correlations in the shape variables and, hence, estimation ofcovariance structure is biased [Kent (1994)]. Hence, in general, Procrustestechniques are preferred over edge superimposition techniques in shape anal-ysis. Hence, this would mean in the current context that the Procrustes ap-proaches of this paper should be preferred to inference using the Choleskydecomposition.
6. Applications.
Interpolation of covariance matrices.
Frequently in diffusion tensorimaging one wishes to carry out interpolation between tensors. When thetensors are quite different, interpolation using different metrics can lead tovery different results. For example, consider Figure 3, where four differentgeodesic paths are plotted between two tensors. Arsigny et al. (2007) notethat the Euclidean metric is prone to swelling, which is seen in this example.Also, the log-Euclidean metric gives strong weight to small volumes. In thisexample the Cholesky and Procrustes size-and-shape paths look rather dif-ferent, due to the extra rotation in the Procrustes method. From a variety ofexamples it does seem clear that the Euclidean metric is very problematic,especially due to the swelling of the volume. In general, the log-Euclideanand Procrustes size-and-shape methods seem preferable.In some applications, for example, fiber tracking, we may need to inter-polate between several covariance matrices on a grid, in which case we canuse weighted Fr´echet meansˆΣ = arg inf Σ n X i =1 w i d ( S i , Σ) , n X i =1 w i = 1 , I. L. DRYDEN, A. KOLOYDENKO AND D. ZHOU
Fig. 4.
Demonstration of PCA for covariance matrices. The true geodesic path is givenin the penultimate row (black). We then add noise in the three initial rows (red). Then weestimate the mean and find the first principal component (yellow), displayed in the bottomrow. where the weights w i are proportional to a function of the distance (e.g.,inverse distance or Kriging based weights).6.2. Principal components analysis of diffusion tensors.
We consider nowan example estimating the principal geodesics of the covariance matrices S , . . . , S n using the Procrustes size-and-shape metric. The data are dis-played in Figure 4 and here k = 3. We consider a true geodesic path (black)and evaluate 11 equally spaced covariance matrices along this path. Wethen add noise for three separate realizations of noisy paths (in red). Thenoise is independent and identically distributed Gaussian and is added inthe dual space of the tangent coordinates. First, the overall Procrustes size-and-shape mean ˆΣ S is computed based on all the data ( n = 33), and thenthe Procrustes size-and-shape tangent space co-ordinates are obtained. Thefirst principal component loadings are computed and projected back to givean estimated minimal geodesic in the covariance matrix space. We plot thispath in yellow by displaying 11 covariance matrices along the path. As wewould expect, the first principal component path bears a strong similarityto the true geodesic path. The percentages of variability explained by thefirst three PCs are as follows: PC1 (72.0%), PC2 (8.8%), PC3 (6.5%).The data can also be seen in the dual Procrustes space of 4 points in k = 3dimensions in Figure 5. We also see the data after applying the Procrustesfitting, we show the effects of the first three principal components, and alsothe plot of the first three PC scores. ON-EUCLIDEAN STATISTICS FOR COVARIANCE MATRICES Fig. 5. (top left) The noisy configurations in the dual space of k + 1 = 4 points in k = 3 dimensions. For each configuration point 1 is colored black, point 2 is red, point 3 is greenand point 4 is blue, and the points in a configuration are joined by lines. (top right) TheProcrustes registered data, after removing translation, rotation and reflection. (bottom left)The Procrustes mean size-and-shape, with vectors drawn along the directions of the firstthree PCs (PC1—black, PC2—red, PC3—green). (bottom right) The first three PC scores.The points are colored by the position along the true geodesic from left to right (black, red,green, blue, cyan, purple, yellow, grey, black, red, green). Interpolation.
We consider the interpolation of part of the brainimage in Figure 2. In Figure 6(a) we see the original FA image, and in Figure6(b) and (c) we see interpolated images using size-and-shape distance. Theinterpolation is carried out at two equally spaced points between voxels,and Figure 6(b) shows the FA image from the interpolation and Figure 6(c)shows the PA image. In the bottom right plot of Figure 6 we highlightthe selected regions in the box. It is clear that the interpolated images aresmoother, and it is clear from the anisotropy maps of the interpolated datathat the cingulum (cg) is distinct from the corpus callosum (cc).6.4.
Anisotropy.
As a final application we consider some diffusion ten-sors obtained from diffusion weighted images in the brain. In Figure 7 wesee a coronal slice from the brain with the 3 × I. L. DRYDEN, A. KOLOYDENKO AND D. ZHOU
Fig. 6.
FA maps from the original (a) and interpolated (b) data. In (c) the PA mapis displayed, and in (a1) , (b1) , (c1) we see the zoomed in regions marked in (a) , (b) , (c) respectively. internal capsule and on the lower right we see the superior fronto-occipitalfasciculus. Fig. 7.
In the upper plots we see the anisotropy measures (left) FA, (middle) PA, (right)GA. In the lower plot we see the diffusion tensors, which have been scaled to have volumeproportional to √ F A . ON-EUCLIDEAN STATISTICS FOR COVARIANCE MATRICES At first sight all three measures appear broadly similar. However, thePA image offers more contrast than the FA image in the highly anisotropicregion—the corpus callosum. Also, the GA image has rather fewer brighterareas than PA or FA. Due to the improved contrast, we believe PA is slightlypreferable in this example.6.5.
Simulation study.
Finally, we consider a simulation study to com-pare the different estimators. We consider the problem of estimating a pop-ulation covariance matrix Ω from a random sample of k × k covariance ma-trices S , . . . , S n .We consider a random sample generated as follows. Let ∆ = chol(Ω) and X i be a random matrix with i.i.d. entries with E [( X i ) jl ] = 0 , var(( X i ) jl ) = σ , i = 1 , . . . , n ; j = 1 , . . . , k ; l = 1 , . . . , k . We take S i = (∆ + X i )(∆ + X i ) T , i = 1 , . . . , n. We shall consider four error models:I. Gaussian square root: ( X i ) jl are i.i.d. N (0 , σ ) for j = 1 , . . . , k ; l = 1 , . . . , k .II. Gaussian Cholesky: ( X i ) jl are i.i.d. N (0 , σ ) for j ≤ k and zero other-wise.III. Log-Gaussian: i.i.d. Gaussian errors N (0 , σ ) are added to the matrixlogarithm of ∆ to give Y , and then the matrix exponential of Y Y T istaken.IV. Student’s t with 3 degrees of freedom: ( X i ) jl are i.i.d. ( σ/ √ t for j = 1 , . . . , k ; l = 1 , . . . , k .We consider the performance in a simulation study, with 1000 MonteCarlo simulations. The results are presented in Tables 2 and 3 for two choicesof population covariance matrix. We took k = 3 and n = 10 ,
30. In order toinvestigate the efficiency of the estimators, we use three measures: estimatedmean square error between the estimate and the matrix Ω with metrics d E ( · ) , d S ( · ) and the estimated risk from using Stein loss [James and Stein(1961)] which is given by L ( S , S ) = trace( S S − ) − log det( S S − ) − k, where det( · ) is the determinant. Clearly the efficiency of the methods de-pends strongly on the Ω and the error distribution.Consider the first case where the mean has λ = 1 , λ = 0 . , λ = 0 . n = 10, with ˆΣ H performing the best. For n = 30 either ˆΣ H or ˆΣ S arebetter, with ˆΣ E performing least well. For model II with Gaussian errorsadded in the Cholesky decomposition we see that ˆΣ C is the best, althoughthe other estimators are quite similar, with the exception of ˆΣ E which is I. L. DRYDEN, A. KOLOYDENKO AND D. ZHOU worse. For model III with Gaussian errors on the matrix logarithm scale allestimators are quite similar, as the variability is rather small. The estimateˆΣ R is slightly better here than the others. For model IV with Student’s t errors we see that ˆΣ H and ˆΣ S are slightly better on the whole, although ˆΣ E is again the worst performer. Table 2
Measures of efficiency, with k = 3 and σ = 0 . . RMSE is the root mean square errorusing either the Euclidean norm or the Procrustes size-and-shape norm, and “Stein”refers to the risk using the Stein loss function. The smallest value in each row ishighlighted in bold. The mean has parameters λ = 1 , λ = 0 . , λ = 0 . . The errordistributions for Models I–IV are Gaussian (square root), Gaussian (Cholesky),log-Gaussian and Student’s t , respectively ˆΣ E ˆΣ C ˆΣ S ˆΣ H ˆΣ L ˆΣ R ˆΣ F I n = 10 RMSE( d E ) 0.1136 0.1057 0.104 d S ) 0.0911 0.082 0.0802 n = 30 RMSE( d E ) 0.0788 0.0669 0.0626 d S ) 0.0691 0.0516 n = 10 RMSE( d E ) 0.0973 d S ) 0.0797 n = 30 RMSE( d E ) 0.0641 d S ) 0.0585 n = 10 RMSE( d E ) 0.0338 0.0333 0.0336 0.0335 0.0333 d S ) 0.0195 0.0193 0.0194 0.0194 0.0192 n = 30 RMSE( d E ) 0.0329 0.0324 0.0327 0.0327 0.0324 d S ) 0.0187 0.0184 0.0185 0.0185 0.0183 n = 10 RMSE( d E ) 0.119 0.1012 0.1006 d S ) 0.1202 0.082 0.0818 n = 30 RMSE( d E ) 0.081 0.0618 0.0598 d S ) 0.0828 0.0489 Table 3
Measures of efficiency, with k = 3 and σ = 0 . . RMSE is the root mean square errorusing either the Euclidean norm or the Procrustes size-and-shape norm, and “Stein”refers to the risk using the Stein loss function. The smallest value in each row ishighlighted in bold. The mean has parameters λ = 1 , λ = 0 . , λ = 0 . . The errordistributions for Models I–IV are Gaussian (square root), Guassian (Cholesky),log-Gaussian and Student’s t , respectively ˆΣ E ˆΣ C ˆΣ S ˆΣ H ˆΣ L ˆΣ R ˆΣ F I n = 10 RMSE( d E ) 0.0999 0.2696 0.0894 d S ) 0.2091 0.2172 0.1424 0.1491 n = 30 RMSE( d E ) 0.0708 0.2836 0.0552 d S ) 0.2064 0.2112 0.1301 0.1388 v 0.3484 0.1317Stein 53.3301 25.8512 22.2974 25.378 n = 10 RMSE( d E ) 0.0907 0.4879 0.0844 d S ) 0.1669 0.3571 0.1139 0.1176 n = 30 RMSE( d E ) 0.0606 0.5151 0.0509 d S ) 0.1632 0.3369 0.1022 0.1067 n = 10 RMSE( d E ) 0.0315 0.0312 0.0313 0.0313 d S ) 0.0162 0.016 0.0161 0.0161 0.016 n = 30 RMSE( d E ) 0.031 0.0307 0.0309 0.0309 0.0306 d S ) 0.0156 0.0154 0.0155 0.0155 0.0154 IV n = 10 RMSE( d E ) 0.1055 0.2519 0.0848 d S ) 0.2187 0.197 0.1253 0.1301 n = 30 RMSE( d E ) 0.0755 0.2628 0.0523 d S ) 0.2098 0.186 0.1089 0.1161 In Table 3 we now consider the case λ = 1 , λ = 0 . , λ = 0 . C and Σ R can behave quite poorly in this example, when using RMSE ( d E ) or RMSE ( d S ) for assessment. This is particularly noticeable in the simulationsfor models I, II and IV. The better estimators are generally ˆΣ H , ˆΣ S and ˆΣ L ,with ˆΣ E a little inferior. I. L. DRYDEN, A. KOLOYDENKO AND D. ZHOU
Overall, in these and other simulations ˆΣ H , ˆΣ S and ˆΣ L have performedconsistently well.
7. Discussion.
In this paper we have introduced new methods and re-viewed recent developments for estimating a mean covariance matrix wherethe data are covariance matrices. Such a situation appears to be increasinglycommon in applications.Another possible metric is the power Euclidean metric d A ( S , S ) = 1 α k S α − S α k , (10)where S α = U Λ α U T . We have considered α ∈ { / , } earlier. As α → , themetric approaches the log-Euclidean metric. We could consider any nonzero α ∈ R depending on the situation, and the estimate of the covariance matrixwould beˆΣ A = ( ˆ∆ A ) /α , where ˆ∆ A = arg inf ∆ ( n X i =1 k S αi − ∆ k ) = 1 n n X i =1 S αi . For positive α the estimators become more resistant to outliers for smaller α , and for larger α the estimators become less resistant to outliers. Fornegative α one is working with powers of the inverse covariance matrix.Also, one could include the Procrustes registration if required. The resultingfractional anisotropy measure using the power metric (10) is given by FA ( α ) = ( kk − k X i =1 ( λ αi − λ α ) (cid:30) k X i =1 λ αi ) / , and λ α = k P ki =1 λ αi . A practical visualization tool is to vary α in order fora neurologist to help interpret the white fiber tracts in the images.We have provided some new methods for estimation of covariance matri-ces which are themselves rooted in statistical shape analysis. Making thisconnection also means that methodology developed from covariance ma-trix analysis could also be useful for applications in shape analysis. Thereis much current interest in high-dimensional covariance matrices [cf. Bickeland Levine (2008)], where k ≫ n . Sparsity and banding structure often areexploited to improve estimation of the covariance matrix or its inverse. Mak-ing connections with the large amount of activity in this field should alsolead to new insights in high-dimensional shape analysis [e.g., see Dryden(2005)].Note that the methods of this paper also have potential applications inmany areas, including modeling longitudinal data. For example, Choleskydecompositions are frequently used for modeling longitudinal data, bothwith Bayesian and random effect models [e.g., see Daniels and Kass (2001); ON-EUCLIDEAN STATISTICS FOR COVARIANCE MATRICES Chen and Dunson (2003); Pourahmadi (2007)]. The Procrustes size-and-shape metric and matrix square root metric provide a further opportunityfor modeling, and may have advantages in some applications, for example, incases where the covariance matrices are close to being deficient in rank. Fur-ther applications where deficient rank matrices occur are structure tensors incomputer vision. The Procrustes approach is particularly well suited to suchdeficient rank applications, for example, with structure tensors associatedwith surfaces in an image. Other application areas include the averagingof affine transformations [Alexa (2002); Aljabar et al. (2008)] in computergraphics and medical imaging. Also the methodology could be useful in com-putational Bayesian inference for covariance matrices using Markov chainMonte Carlo output. One wishes to estimate the posterior mean and othersummary statistics from the output, and that the methods of this paper willoften be more appropriate than the usual Euclidean distance calculations.
Acknowledgments.
We wish to thank the anonymous reviewers and Huil-ing Le for their helpful comments. We are grateful to Paul Morgan (MedicalUniversity of South Carolina) for providing the brain data, and to Bai Li,Dorothee Auer, Christopher Tench and Stamatis Sotiropoulos, from the EUfunded CMIAG Centre at the University of Nottingham, for discussions re-lated to this work. REFERENCES
Alexa, M. (2002). Linear combination of transformations.
ACM Trans. Graph. Alexander, D. C. (2005). Multiple-fiber reconstruction algorithms for diffusion MRI.
Ann. N. Y. Acad. Sci.
Aljabar, P., Bhatia, K. K., Murgasova, M., Hajnal, J. V., Boardman, J. P.,Srinivasan, L., Rutherford, M. A., Dyet, L. E., Edwards, A. D. and
Rueckert,D. (2008). Assessment of brain growth in early childhood using deformation-basedmorphometry.
Neuroimage Amaral, G. J. A., Dryden, I. L. and
Wood, A. T. A. (2007). Pivotal bootstrapmethods for k -sample problems in directional statistics and shape analysis. J. Amer.Statist. Assoc.
Arsigny, V., Fillard, P., Pennec, X. and
Ayache, N. (2007). Geometric means in anovel vector space structure on symmetric positive-definite matrices.
SIAM J. MatrixAnal. Appl. Basser, P. J., Mattiello, J. and
Le Bihan, D. (1994). Estimation of the effectiveself-diffusion tensor from the NMR spin echo.
J. Magn. Reson. B
Basu, S., Fletcher, P. T. and
Whitaker, R. T. (2006). Rician noise removal in diffu-sion tensor MRI. In
MICCAI (1) (R. Larsen, M. Nielsen and J. Sporring, eds.)
LectureNotes in Computer Science
Bhattacharya, R. and
Patrangenaru, V. (2003). Large sample theory of intrinsic andextrinsic sample means on manifolds. I.
Ann. Statist. Bhattacharya, R. and
Patrangenaru, V. (2005). Large sample theory of intrinsic andextrinsic sample means on manifolds. II.
Ann. Statist. I. L. DRYDEN, A. KOLOYDENKO AND D. ZHOU
Bickel, P. J. and
Levina, E. (2008). Regularized estimation of large covariance matrices.
Ann. Statist. Bookstein, F. L. (1986). Size and shape spaces for landmark data in two dimensions(with discussion).
Statist. Sci. Chen, Z. and
Dunson, D. B. (2003). Random effects selection in linear mixed models.
Biometrics Daniels, M. J. and
Kass, R. E. (2001). Shrinkage estimators for covariance matrices.
Biometrics Daniels, M. J. and
Pourahmadi, M. (2002). Bayesian analysis of covariance matricesand dynamic models for longitudinal data.
Biometrika Dryden, I. L. (2005). Statistical analysis on high-dimensional spheres and shape spaces.
Ann. Statist. Dryden, I. L. and
Mardia, K. V. (1998).
Statistical Shape Analysis . Wiley, Chichester.MR1646114
Fillard, P., Arsigny, V., Pennec, X. and
Ayache, N. (2007). Clinical DT-MRI es-timation, smoothing and fiber tracking with log-Euclidean metrics.
IEEE Trans. Med.Imaging Fletcher, P. T. and
Joshi, S. (2007). Riemannian geometry for the statistical analysisof diffusion tensor data.
Signal Process. Fr´echet, M. (1948). Les ´el´ements al´eatoires de nature quelconque dans un espace dis-tanci´e.
Ann. Inst. H. Poincar´e Goodall, C. R. and
Mardia, K. V. (1992). The noncentral Bartlett decompositionsand shape densities.
J. Multivariate Anal. Gower, J. C. (1975). Generalized Procrustes analysis.
Psychometrika James, W. and
Stein, C. (1961). Estimation with quadratic loss. In
Proc. 4th BerkeleySympos. Math. Statist. and Prob., Vol. I
Karcher, H. (1977). Riemannian center of mass and mollifier smoothing.
Comm. PureAppl. Math. Kendall, D. G. (1984). Shape manifolds, Procrustean metrics and complex projectivespaces.
Bull. London Math. Soc. Kendall, D. G. (1989). A survey of the statistical theory of shape.
Statist. Sci. Kendall, D. G., Barden, D., Carne, T. K. and
Le, H. (1999).
Shape and ShapeTheory . Wiley, Chichester. MR1891212
Kendall, W. S. (1990). Probability, convexity, and harmonic maps with small image. I.Uniqueness and fine existence.
Proc. London Math. Soc. (3) Kent, J. T. (1994). The complex Bingham distribution and shape analysis.
J. Roy. Statist.Soc. Ser. B Le, H. (2001). Locating Fr´echet means with application to shape spaces.
Adv. in Appl.Probab. Le, H. and
Small, C. G. (1999). Multidimensional scaling of simplex shapes.
PatternRecognition Le, H.-L. (1988). Shape theory in flat and curved spaces, and shape densities with uniformgenerators. Ph.D. thesis, Univ. Cambridge.
Le, H.-L. (1992). The shapes of non-generic figures, and applications to collinearity test-ing.
Proc. Roy. Soc. London Ser. A
Le, H.-L. (1995). Mean size-and-shapes and mean shapes: A geometric point of view.
Adv. in Appl. Probab. Lenglet, C., Rousson, M. and
Deriche, R. (2006). DTI segmentation by statisticalsurface evolution.
IEEE Trans. Med. Imaging Mardia, K. V., Kent, J. T. and
Bibby, J. M. (1979).
Multivariate Analysis . AcademicPress, London. MR0560319
Moakher, M. (2005). A differential geometric approach to the geometric mean of sym-metric positive-definite matrices.
SIAM J. Matrix Anal. Appl. Pennec, X. (1999). Probabilities and statistics on Riemannian manifolds: Basic toolsfor geometric measurements. In
Proceedings of IEEE Workshop on Nonlinear Signaland Image Processing (NSIP99) (A. Cetin, L. Akarun, A. Ertuzun, M. Gurcan and Y.Yardimci, eds.) Pennec, X., Fillard, P. and
Ayache, N. (2006). A Riemannian framework for tensorcomputing.
Int. J. Comput. Vision Pourahmadi, M. (2007). Cholesky decompositions and estimation of a covariance ma-trix: Orthogonality of variance correlation parameters.
Biometrika R Development Core Team (2007).
R: A Language and Environment for StatisticalComputing . R Foundation for Statistical Computing, Vienna, Austria.
Schwartzman, A. (2006). Random ellipsoids and false discovery rates: Statistics for dif-fusion tensor imaging data. Ph.D. thesis, Stanford Univ.
Schwartzman, A., Dougherty, R. F. and
Taylor, J. E. (2008). False discovery rateanalysis of brain diffusion direction maps.
Ann. Appl. Statist. Wang, Z., Vemuri, B., Chen, Y. and
Mareci, T. (2004). A constrained variationalprinciple for direct estimation and smoothing of the diffusion tensor field from complexDWI.
IEEE Trans. Med. Imaging Zhou, D., Dryden, I. L., Koloydenko, A. and
Bai, L. (2008). A Bayesian methodwith reparameterisation for diffusion tensor imaging. In
Proceedings, SPIE conference.Medical Imaging 2008: Image Processing (J. M. Reinhardt and J. P. W. Pluim, eds.)69142J. SPIE, Bellingham, WA.
I. L. DrydenDepartment of StatisticsLeConte CollegeUniversity of South CarolinaColumbia, South Carolina 29208USAandSchool of Mathematical SciencesUniversity of NottinghamUniversity ParkNottingham, NG7 2RDUKE-mail: [email protected]
A. KoloydenkoDepartment of MathematicsRoyal Holloway, University of LondonEgham, TW20 0EXUKE-mail: [email protected]