Variance Measures for Symmetric Positive (Semi-) Definite Tensors in Two Dimensions

Calculating the variance of a family of tensors, each represented by a symmetric positive semi-definite second order tensor/matrix, involves the formation of a fourth order tensor Rabcd. To form this tensor, the tensor product of each second order tensor with itself is formed, and these products are then summed, giving the tensor Rabcd the same symmetry properties as the elasticity tensor in continuum mechanics. This tensor has been studied with respect to many properties: representations, invariants, decomposition, the equivalence problem et cetera. In this paper we focus on the two-dimensional case where we give a set of invariants which ensures equivalence of two such fourth order tensors Rabcd and R~abcd. In terms of components, such an equivalence means that components Rijkl of the first tensor will transform into the components R~ijkl of the second tensor for some change of the coordinate system.

All such second order tensors share the same mathematical properties, namely, they are realvalued, symmetric, and positive semi-definite. Moreover, in these disciplines, one encounters a collection of such tensors, e.g., at different locations of the image. Populations of such tensors have also been key to some studies aiming to model the underlying structure of the medium under investigation [8,12,18].
Irrespective of the particular application, let R ab denote such tensors, 1 and we shall refer to the set of n tensors as {R ab (i) } i . Our desire is to find relevant descriptors or models of such a family. One relevant statistical measure of this family is the (population) variance where R ab = 1 n ∑ i = 1 n R ab (i) is the mean. (For another approach, see e.g., [8]). In this paper, we are interested in the first term, i.e., we study the fourth order tensor (skipping the normalization) where R ab (i) ≥ 0 stands for R ab (i) being positive semi-definite. It is obvious that R abcd has the symmetries R abcd = R bacd = R abdc and R abcd = R cdab , i.e., R abcd has the same symmetries as the elasticity tensor [14] from continuum mechanics. The elasticity tensor is well studied [13], e.g. with respect to classification, decompositions, and invariants. In most cases this is done in three dimensions. The same (w.r.t. symmetries) tensor has also been studied in the context of diffusion MR [2].
In this paper we will focus on the corresponding tensor R abcd in two dimensions. First, there are direct applications in image processing, and secondly, the problems posed will be more accessible in two dimensions than in three. In particular we study the equivalence problem, namely, we ask the question: given the components R ijkl and R ijkl of two such tensors do they represent the same tensor in different coordinate systems (see Sects. 2.1.2 and 4)?

Outline
Section 2 contains tensorial matters. We will assume some basic knowledge of tensors, although some definitions are given for completeness. The notation(s) used is commented on and in particular the three-dimensional Euclidean vector space V (ab) is introduced.
In Sect. 2.1.2, we make some general remarks concerning the tensor R abcd and specify the problem we focus on. Section 2.1 is concluded with some remarks on the Voigt/Kelvin notation and the corresponding visualisation in ℝ 3 .
Section 2.2 gives examples of invariants, especially invariants which are easily accessible from R abcd . Also, more general invariant/canonical decompositions of R abcd are given.
In Sect. 3, we discuss how the tensor R abcd can (given a careful choice of basis) be expressed in terms of a 3 × 3 matrix, and how this matrix is affected by a rotation of the coordinate system in the underlying two-dimensional space on which R abcd is defined.
In Sect. 4 we return to the equivalence problem and give the main result of this work. In Sect. 4.1.1 we provide a geometric condition for equivalence, while in Sect. 4.1.2, we present the equivalence in terms of a 3 × 3 matrix. Both these characterisations rely on the choice of particular basis elements for the vector spaces employed. In Sect. 4.1.3 the same equivalence conditions are given in a form which does not assume a particular basis.

Preliminaries
In this section we clarify the notation and some concepts which we need. Section 2.1 deals with the (alternatives of) tensor notation and some representations. The equivalence (and related) problems are also briefly addressed. Section 2.2 accounts for some natural invariants, traces and decompositions of R abcd .
We will assume some familiarity with tensors, but to clarify the view on tensors we recall some facts. We start with a (finite dimensional) vector space V with dual V*. A tensor of order (p,q) is then a multi-linear mapping (non-degenerate) metric/scalar product g : V × V ℝ gives an isomorphism from V to V* through v → g(v, ·), and it is this isomorphism which is used to 'raise and lower indices', see below. Indeed, for a fixed v ∈ V, g(v, ·) is a linear mapping V ℝ, i.e., an element of V*.

Tensor Notation and Representations
There is a plethora of notations for tensors. Here, we follow the well-adopted convention [16] that early lower case Latin letters (T a bc ) refer to the tensor as a geometric object, its type being inferred from the indices and their positions (the abstract index notation). g ab denotes the metric tensor. When the indices are lower case Latin letters from the middle of the alphabet, T i jk , they refer to components of T a bc in a certain frame. The super-index i denotes a contravariant index while the sub-indices j, k are covariant. For instance, a typical vector (tensor of type (1, 0)) will be written v a with components v i , while the metric g ab (tensor of type (0, 2)) has components g ij . At a number of occasions, it will also be useful to express quantities in terms of components with respect to orthonormal frames, i.e., Cartesian coordinates. This is sometimes referred to as 'Cartesian tensors', and the distinction between contra-and covariant indices is obscured. In these situations, it is possible (but not necessary) to write all indices as sub-indices, and sometimes the symbol ≐ is used to indicate that an equation is only valid in Cartesian coordinates. For example T i ≐ T ijk δ jk instead of T i = T i jk g jk = T ik k . Often this is clear form the context, but we will sometimes use ≐ to remind the reader that a Cartesian assumption is made. Here, the Einstein summation convention is implied, i.e., repeated indices are to be summed over, so that for instance T i = T jk i g jk = T k ik = ∑ j = 1 n ∑ k = 1 n T jk i g jk = ∑ k = 1 n from 1 to n. We have also used the metric g ij and its inverse g ij to raise and lower indices. For instance, since g ij v i is an element of V*, we write g ij v i = v j .
We also remind of the notation for symmetrisation. For a two-tensor T (ab) = 1 2 (T ab + T ba ), while more generally for a tensor T a 1 a 2 ⋯a n of order (0, n) we have T (a 1 a 2 ⋯a n ) = 1 n! ∑ π T a π(1) a π(2) ⋯a π(n) where the sum is taken over all permutations π of 1, 2, …, n. Naturally, this convention can also be applied to subsets of indices. For instance, H a(bc) = 1 2 (H abc + H acb ).

The Vector
Space of Symmetric Two-Tensors-In any coordinate frame a symmetric tensor R ab (i.e., R ab = R ba ) is represented by a symmetric matrix R ij (2 × 2 or 3 × 3 depending on the dimension of the underlying space). In the two-dimensional case, with the underlying vector space V a ∼ ℝ 2 , this means that R ab lives in a three-dimensional vector space, which we denote by V (ab) . V (ab) is equipped with a natural scalar product: < A ab , B ab >= A ab B ab , making it into a three-dimensional Euclidean space. Here A ab B ab = A ab B cd g ac g bd , i.e, the contraction of A ab B cd over the indices a, c and b, d, and the tensor product A ab B cd itself is the tensor of order (0, 4) given by (A ab B cd )v a u b w c m d = (A ab v a u b ) (B cd w c m d ) together with multi-linearity.

The
Tensor R abcd and the Equivalence Problem-As noted above, R abcd given by (1) has the symmetries R abcd = R (ab)cd = R ab(cd) and R abcd = R cdab , and it is not hard to see that this gives R abcd six degrees of freedom in two dimensions. (See also Sect. 2.1.3.) It is also interesting to note that R abcd provides a mapping V (ab) → V (ab) through R ab R abcd R cd , and that this mapping is symmetric (due to the symmetry R abcd = R cdab ). Given R abcd there are a number of questions one can ask, e.g., • Feasibility-given a tensor R abcd with the correct symmetries, can it be written in the form (1)?
• Canonical decomposition-given R abcd of the form (1), can you write R abcd as a canonical sum of the form (1), but with a fixed number of terms (cf. eigenvector decomposition of symmetric matrices)?
• Visualisation-since fourth order tensors are a bit involved, how can one visualise them in ordinary space?
• Characterisation/relevant sets of invariants-what invariants are relevant from an application point of view?
• The equivalence problem-in terms of components, how do we know if R ijkl and R ijkl represent the same tensor when they are in different coordinate systems?
We will now focus on the equivalence problem in two dimensions. This problem can be formulated as above: given, in terms of components, two tensors (with the symmetries we consider) R ijkl and R ijkl , do they represent the same tensor in the sense that there is a coordinate transformation taking the components R ijkl into the components R ijkl ? In other words, does there exist an (invertible) matrix P m i so that This problem can also be formulated when R ijkl and R ijkl are expressed in Cartesian frames. Then the coordinate transformation must be a rotation, i.e., given by a rotation matrix Q i j ∈ SO(2). Hence, the problem of (unitary) equivalence is: Given R ijkl and R ijkl , both expressed in Cartesian frames, is there a matrix (applying the 'Cartesian convention') Q ij ∈ SO(2) so that R ijkl = R mnop Q mi Q nj Q ok Q pl ?

The Voigt/Kelvin Notation-Since
and use vector algebra on ℝ 3 . This is used in the Voigt notation [15] and the related Kelvin notation [6]. As always, one must be careful to specify with respect to which basis in V (ab) the coordinates (ξ a ξ b + η a η b ), i.e., in the induced basis we have In this basis, we write an arbitrary element M ab ∈ V (ab) as M ij = z + x y y z − x , which means that M ab gets the coordinates M i = 2 x y z . Note that M ij is positive definite if z 2 − x 2 − y 2 ≥ 0 and z ≥ 0. In terms of the coordinates of the Voigt notation, the tensor R abcd corresponds to a symmetric mapping ℝ 3 ℝ 3 , given by a symmetric 3 × 3 matrix, which also shows that the degrees of freedom for R abcd is six.

Visualization in
ℝ 3 -Through the Voigt notation, any symmetric two-tensor (in two dimensions) can be visualised as a vector in ℝ 3 . Using the basis vector given by (2), we note that e ij (1) and e ij (2) correspond to indefinite quadratic forms, while e ij (3) is positive definite.
In Fig. 1 (left) these matrices are illustrated as vectors in ℝ 3 . The set of positive semidefinite matrices corresponds to a cone, cf. [4], indicated in blue. When the symmetric 2 × 2 matrices are viewed as vectors in ℝ 3 , the outer product of such a vector with itself gives a symmetric 3 × 3 matrix. Hence we get a positive semi-definite quadratic form on ℝ 3 , which can be illustrated by an (degenerate) ellipsoid in ℝ 3 . In Fig. 1 (right) (e ab (1) + e ab (3) )(e cd (1) + e cd (3) ), (e ab (2) + e ab (3) )(e cd (2) + e cd (3) ) and e ab

Invariants, Traces and Decompositions
By an invariant, we mean a quantity that can be calculated from measurements, and which is independent of the frame/coordinate system with respect to which the measurements are performed, despite the fact that components, e.g., T i jk themselves depend on the coordinate system. It is this property that makes invariants important, and typically they are formed via tensor products and contractions, e.g., T i jk T k il g jl . Sometimes, the invariants have a direct geometrical meaning. For instance, for a vector v i , the most natural invariant is its squared length v i v i . For a tensor T i j of order (1,1) in three dimensions, viewed as a linear mapping ℝ 3 ℝ 3 , the most well known invariants are perhaps the trace T i i and the determinant det(T i j ). The modulus of the determinant gives the volume scaling under the mapping given by T i j , while the trace equals the sum of the eigenvalues. If T i j represents a rotation matrix, then its trace is 1 + 2 cos ϕ, where ϕ is the rotation angle. In general, however, the interpretation of a given invariant may be obscure. (For an account relevant to image processing, see e.g., [9]. A different, but relevant, approach in the field of diffusion MRI is found in [20].) (1), and considering the symmetries of R abcd , two (and only two) natural traces arise. For a tensor of order (1, 1), e.g., R i j , it is natural to consider this as an ordinary matrix, and consequently use stem letters without any indices at all. To indicate this slight deviation from the standard tensor notation, we denote e.g., R i j by R . Using [·] for the trace, so that [R ] = Tr(R ) = R a a , we then have

Natural Traces and Invariants-From
and Hence, in a Cartesian frame, where the index position is unimportant, we have for the To proceed there are two double traces (i.e., contracting R abcd twice): and S = S a a = R ac In two dimensions, the difference T ab −S ab is proportional to the metric g ab . Namely, Lemma 1 With T ab and S ab given by (3) and (4), it holds that (in two dimensions) T ab − S ab = ∑ i = 1 n det(R (i) )g ab .
Proof By linearity, it is enough to prove the statement when n = 1, i.e., when the sum has just one term. Raising the second index, and using components, the statement then The first inequality becomes equality when ac = 0, i.e., when A has rank one. The second inequality becomes equality when a = c, i.e., when A is isotropic. □

Definition 1
We define the mean rank, r m , by r m = T/S, with T and S as above. Hence, in two dimensions, 1 ≤ r m ≤ 2.

A Canonical
Decomposition-It is customary [3,7] to decompose a tensor with the symmetries of R abcd into a sum where one term is the completely symmetric part:

Lemma 3
Suppose that W abcd is given by (7), and put W ab = W abcd g cd , W = W ab g ab . Then (in two dimensions) W abcd = W 4 (2g ab g cd − g ac g bd − g ad g bc ) Proof By linearity, it is enough to consider the case when R abcd = A ab A cd for some (symmetric) A ab . In terms of eigenvectors (to A a b ) we can write A ab = αx a x b + βy a y b , where x a x a = y a y a = 1, x a y a = 0. In particular g ab = x a x b + y a y b . From (7) we then get W abcd = 1 3 (2R abcd − R acbd − R adbc ) = 1 3 (2A ab A cd − A ac A bd − A ad A bc ) = 1 3 (2(αx a x b + βy a y b )(αx c x d + βy c y d ) − (αx a x c + βy a y c )(αx b x d + βy b y d ) − (αx a x d + βy a y d )(αx b x c + βy b y c )) . = αβ 3 (2g ab g cd − g ac g bd − g ad g bc ) , (9) where the last equality can be seen by inserting g ab = x a x b + y a y b (for all indices) and expanding. Taking one trace, i.e., contracting with g cd gives W ab = 2αβ 3 g ab , and another trace gives W = 4αβ 3 , which proves the lemma. □

R abcd as a Quadratic Form on ℝ 3
Through the orthonormal basis for the space of symmetric two-tensors (in two dimensions) given by (2), the tensor R abcd viewed as a quadratic form can be represented by a 3 × 3-matrix. Here, we will restrict ourselves to an orthonormal basis for V (ab) , namely the basis It is instructive to see how the various derived tensors show up in M ij . In terms of the basis (2) it is natural to look at the various parts of M ij as follows This splitting is natural for reasons which will become apparent in the next sections. Note, however, that with this representation it is tempting to consider coordinate changes in ℝ 3 , which is not natural in this case. Rather, of interest is the change of basis in V a and the related induced change of coordinates in the representation (10). See Sect. 3.2.

Representation of the Canonically Derived Parts of R abcd
It is helpful to see how the components of the various tensors T ab , S ab , T, S, H , (12) where the last equality follows form the trace-freeness of e ab (1) and e ab (2) . This means that the components of T ∘ ab (properly rescaled) goes into M ij as the components of v (and v t ) in (10 . When M ij is constructed from R abcd which is an outer product R ab R cd the trace is given by In addition, the relations below Eq. (7) show that The two degres of freedom in Å corresponds to the two degrees of freedom in H ∘ abcd .

The Behaviour of M ij Under a Rotation of the Coordinate System in V a
The components of M ij are expressed in terms of the orthonormal basis tensors given by (2) This means that for a vector u = ξ η ξ η = (ξ η) ξ η , the coordinates transform as For the components of the basis vectors e ab (1) , e ab (2) , e ab (3) we find (omitting the factor 1 ∕ 2)  (14) and this means that the components M ij transform as But this latter expression is just hence we have the following important remark/observation:

Remark 1
Viewing the matrix M ij as an ellipsoid in ℝ 3 , the effect of a rotation by an angle v in V a corresponds to a rotation of the ellipsoid by an angle 2v around the z-axis in ℝ 3 (where the z-axis corresponds to the 'isotropic direction' given by g ab ).

The Equivalence Problem for R abcd
The equivalence problem for R abcd can be formulated in different ways (for an account in three dimensions, we refer to [3]). Given two tensors R abcd and R abcd , both with the symmetries implied by (1), the question whether they are the same or not is straightforward as one can compare the components in any basis. However, R abcd and R abcd could live in different (but isomorphic) vector spaces, e.g. two tangent spaces at different points, and the concept of equality becomes less clear. On the other hand, in terms of components R ijkl and R ijkl , one could ask whether there is a change of coordinates which takes one set of components into the other. If so, one can find a (invertible) matrix P i j so that R ijkl = R mnop P m i P n j P o k P p l , and the tensors are then said to be equivalent. As already mentioned, it is convenient to restrict the coordinate systems to orthonormal coordinates. This means that two different coordinate systems differ only by their orientation, i.e., the change of coordinates are given by a rotation matrix Q ∈ SO (2). Under the 'Cartesian convention' that all indices are written as subscripts, R abcd and R abcd are equivalent if there is a matrix Q ∈ SO(2) so that (their Cartesian components satisfy) R ijkl = R mnop Q mi Q nj Q ok Q pl .

Different Ways to Characterize the Equivalence of R abcd and R abcd
In this section, we will discuss three ways to determine whether two tensors R abcd and R abcd are equivalent or not. In Sects. 4.1.1 and 4.1.2 we present two such methods briefly, while Sect. 4.1.3, which is more complete, contains the main result of this work.
As mentioned in Sect. 1.1, the results of Sects. 4.1.1 and 4.1.2, which may be used in their own rights, rely on particular choices of basis matrices for V (ab) . The formulation in Sect. 4.1.3 on the other hand, is expressed in the components of R abcd (in any coordinate system) directly.

Orientation of the Ellipsoid in ℝ 3 -A necessary condition for R abcd and R abcd
to be equivalent is that their corresponding 3 × 3-matrices M ij and M ij have the same eigenvalues. On the other hand, this is not sufficient since the representation in ℝ 3 should reflect the freedom in rotating the coordinate system in V a ∼ ℝ 2 . With the coordinates adopted, this corresponds to a rotation of the associated ellipsoid around the z-axis in . This is illustrated in Fig. 2 where three ellipsoids, all representing positive definite symmetric mappings having identical eigenvalues, are shown. The two first ellipsoids can be rotated into each other by a rotation around the z-axis. This implies that the corresponding tensors R abcd and R abcd are equivalent. The third ellipsoid can also be rotated into the two others, but these rotations are around directions other than the z-axis, which means that this ellipsoid represents a different tensor.
In the generic case, with all eigenvalues different, it is easy to test whether two different ellipsoids can be transfered into each other through a rotation around the z-axis. This will be the case if the corresponding eigenvectors (of M ij and M ij ) have the same angle with the z-axis. Hence it is just a matter of checking the z-components of the three normalized eigenvectors and see if they are equal up to sign.

Components in a Canonical
Coordinate System-In a sense, this is the most straightforward method. In a coordinate system which respects e ab (3) as the z-axis in V (ab) ∼ ℝ 3 , two tensors R abcd and R abcd are equivalent if there is a rotation matrix (in two Hence, equivalence can be easily tested by first checking that T = T and that ‖ T If this is the case, (and if ‖ T ∘ ‖ > 0) one determines the rotation matrix Q which gives Examples of invariants are T = R abcd g ab g cd , S = R abcd g ac g bd and the invariants H = H ab g ab , W = W ab g ab . To produce the invariants, we use the tensor R abcd and the metric g ab . However, if we regard V a ∼ ℝ 2 as oriented, so that the orthonormal basis {ξ, η} for V a also is oriented, then invariants can also be formed in another way. Namely, since the space of symmetric 2 × 2 matrices is 3-dimensional, and since the metric g ab singles out a 1-dimensional subspace, it also determines a 2-dimensional subspace L; all elements orthogonal to g ab . This subspace is the set of all symmetric 2 × 2 matrices which are also trace-free. L can be given an orientation through an area form, which in turn inherits the orientation from V a .

Discussion
In this work, we started with a family of symmetric positive (semi-)definite tensors in two dimensions and considered its variance. This lead us to a fourth order tensor R abcd with the same symmetries as the elasticity tensor in continuum mechanics. After listing a number of possible issues to address, we focused on the equivalence problem. Namely, given the components of two such tensors R abcd and R abcd , how can one determine if they represent the same tensor (but in different coordinate systems) or not? In Sect. 4, we saw that this could be investigated in different ways. The result of Theorem 1 is most satisfactory in the sense that it is expressible in terms of the components of the fourth order tensors directly.
There are two natural extensions and/or ways to continue this work. The first is to apply the result to realistic families of e.g., diffusion tensors in two dimensions. The objective is then, apart from establishing possible equivalences, to investigate the geometric meaning of the invariants. The other natural continuation is to investigate the corresponding problem in three dimensions. The degrees of freedom of R abcd will then increase from 6 to 21, leaving us with a substantially harder, but also perhaps more interesting, problem. Left: the symmetric matrices e ab (1) , e ab (2) , e ab (3) (red) and e ab (1) + e ab (3) , e ab (2) + e ab (3) (blue) as vectors in ℝ 3 . The positive semi-definite matrices correspond to vectors which are inside/above the indicated cone (including the boundary). Right: the fourth order tensors (e ab (1) + e ab (3) )(e cd (1) + e cd (3) ) and (e ab (2) + e ab (3) )(e cd (2) + e cd (3) ) depicted in blue, and e ab (3) e cd (3) shown in red are viewed as quadratic forms and illustrated as ellipsoids (made a bit 'fatter' than they should be for the sake of clarity)

Fig. 2.
Three identical (truncated) ellipsoids in ℝ 3 with different orientations. The two leftmost ellipsoids can be carried over to each other through a rotation around the (vertical in the figure) z-axis, which implies that they represent the same tensor R abcd (up to the meaning discussed). The right ellipsoid, despite identical eigenvalues with the two others, represent a different tensor since the rotation which carries this ellipsoid to any of the others is not around the z-axis