Stochastic Modelling of Symmetric Positive Definite Material Tensors

Spatial symmetries and invariances play an important role in the behaviour of materials and should be respected in the description and modelling of material properties. The focus here is the class of physically symmetric and positive definite tensors, as they appear often in the description of materials, and one wants to be able to prescribe certain classes of spatial symmetries and invariances for each member of the whole ensemble, while at the same time demanding that the mean or expected value of the ensemble be subject to a possibly ‘higher’ spatial invariance class. We formulate a modelling framework which not only respects these two requirements—positive definiteness and invariance—but also allows a fine control over orientation on one hand, and strength / size on the other. As the set of positive definite tensors is not a linear space, but rather an open convex cone in the linear space of physically symmetric tensors, we consider it advantageous to widen the notion of mean to the so-called Fréchet mean on a metric space, which is based on distance measures or metrics between positive definite tensors other than the usual Euclidean one. It is shown how the random ensemble can be modelled and generated, independently in its scaling and orientational or directional aspects, with a Lie algebra representation via a memoryless transformation. The parameters which describe the elements in this Lie algebra are then to be considered as random fields on the domain of interest. As an example, a 2D and a 3D model of steady-state heat conduction in a human proximal femur, a bone with high material anisotropy, is modelled with a random thermal conductivity tensor, and the numerical res-ults show the distinct impact of incorporating into the constitutive model different material uncertainties—scaling, orientation, and prescribed material symmetry—on the desired quantities of interest.


Introduction
Heterogeneous materials and their statistical description are of great interest in many fields.The properties of such materials often vary considerably on the spatial scales of interest, and additionally in many cases the material properties are only described statistically, reflecting some underlying uncertainty as to the exact values.This leads to the idea of a stochastic/random representation of these properties.The material properties are naturally collected in tensor quantities, and these tensors are often of even order and physically symmetric-in elasticity this is often termed the major symmetry.Furthermore, when these tensorial properties appear in the definitions of stored energy or entropy production of stable systems, they are additionally positive definite.It is this class of even-order physically symmetric and positive definite (SPD) tensors which we are interested in describing.Well known examples of such tensors are the second-order tensors describing diffusion phenomena, mapping the gradient of some quantity to the flux of some related quantity, e.g. the thermal conductivity tensor mapping temperature gradients to fluxes of thermal energy.An example of a fourth-order tensor is the elasticity tensor of a linearly elastic material, which maps the strain tensor to the stress tensor.Any even-order tensor may be viewed as a linear mapping on the space of tensors of half the order by contracting over half the indices.By choosing a basis in the space of tensors of half the order, each such even-order tensor can hence be represented by a matrix; so the class we are interested in can thus be represented by SPD matrices.We shall switch to this matrix representation whenever it is convenient to do so.
Even though frequently the exact values of the tensorial entries may not be known, some properties resulting from general principles often are known; namely, often it is known what kind of symmetries or invariances are to be expected [50,36].These invariances-we prefer the term invariance over symmetry here to avoid confusion with the afore mentioned physical symmetry, resp.the symmetry of the tensor as a linear mapping-under the operation of some symmetry group are well known, and here we concentrate on spatial point operations such as those defining isotropy.These symmetry groups are represented as linear transformation groups on the relevant space of even order tensors, or, equivalently, on the space of representing symmetric matrices mentioned above.Such transformation groups then define invariant subspaces which represent the tensors with the given invariance property.As the number of such point groups is larger in the case of fourth-order spatial tensors (e.g.[9,7]) than in the case of second-order spatial tensors [41], in our explicit treatment here we shall confine ourselves to second-order spatial tensors for the sake of simplicity and clarity.But the general idea applies to any even-order SPD tensor field.
So, after a choice of bases, what will be considered here are random SPD matrices, where certain invariances are known for the whole population, as well as for the mean.There is already quite a bit of previous work in this field, drawing on several sources, which will be reviewed only very cursory and briefly here.Random SPD matrices (e.g.[30]) arise in a number of fields, e.g.[54,57,29], to name just three.
As each SPD matrix can be spectrally decomposed with positive eigenvalues and has orthogonal eigenvectors, and as this separates size / strength as encoded in the eigenvalues and orientation as encoded in the eigenvectors resp.the orthogonal transformation into the eigenvector basis.We see this as vital information about how the tensor / matrix acts, and we shall base our treatment on this fundamental decomposition.On the other hand, for the purpose of modelling or representing a tensor, the most comfortable and desirable is to find a representation on a linear space.It is one of the aims to show here how this may be achieved for the spectral ecomposition.Positive definite symmetric even-order tensors resp.SPD matrices are an open convex cone in the vector space of symmetric even-order tensors resp.symmetric matrices, and are thus a differentiable manifold but not a vector space.
Differentiable manifolds can be represented-mostly only locally-on their tangent spaces.As with any manifold, the analysis of matrix manifolds is typically facilitated if they can be adorned with a Riemannian structure [52,61,63], and even more so if the manifold carries also a group structure, more specifically in case it is a Lie group, as then everything can be played down to the tangent space at the neutral element i.e. the associated Lie algebra.Unfortunately, the manifold of SPD matrices is not a Lie group under normal matrix multiplication, although it is possible [3] to give the SPD manifold a new commutative or Abelian multiplicative group structure, which coincides with matrix multiplication on commuting matrices, and is lifted from its Lie algebra to the SPD manifold via exponentiation.One may use this map to also carry any Euclidean structure on the Lie algebra to a Riemannian one on the Lie group.
For second-order spatial tensors the spatial group reduction completely determines the eigenvector transformation of the associated SPD matrix, whereas for e.g. for fourth-order tensors there is an iso-spectral sub-group in the relevant Lie group of orthogonal transformations, cf.e.g.[38,39,40,41], which determines the eigen-strain distribution in the relevant elasticity class [50,9,10,36,7].This is one more reason to confine ourselves here initially to second-order spatial tensors for the sake of brevity and simplicity of exposition, the corresponding more involved results for fourth-order spatial tensors (i.e.elasticity tensors) will be published elsewhere.The spectral decomposition of a random SPD matrix has thus two components, the SPD diagonal matrix of eigenvalues and the orthogonal orientation matrix.In particular-as will be shown-this allows the independent control of the strength (or size or norm) of the tensor and of its directional properties.Both components may be random, and thus random orthogonal matrices have some relevance here as well [42,32,45].Thus the spectral decomposition will lead to a representation on a product of Lie algebras-the orthogonal matrices and the positive definite diagonal ones.The Euclidean structure on their joint Lie algebra-the direct sum of the skewsymmetric and diagonal matrices-may be mapped onto a product Riemannian structure on the product manifold of the two Lie groups which represents the SPD matrices.
A Riemannian structure on a manifold as just alluded to allows one to define the length of paths, and thus shortest paths resp.geodesics, and hence to measure distances on the manifold.As the variational definition of mean, resp.average, or expectationwe will use these terms interchangeably-depends on the distance measure (e.g.[48], cf. also Section 2), this leads to the question as to what is an adequate definition of the mean for SPD tensors.The usual arithmetic mean is wedded to a flat (vector-)space and a Euclidean distance, as is well known and will be explained in more detail later in Section 2.3.On a metric space the mean of two points is the point half-way along the shortest path between them.This generalisation is known as the Fréchet mean.This connection between the mean (and variance to be defined later) and a distance metric may be used to switch the discussion of what is desired from the mean to the underlying distance metric, as on a non-flat manifold the Euclidean metric of the ambient space is not advantageous.For example on the Lie group of orthogonal matrices it is completely common to use such a Riemannian metric [42,32,45].
On the manifold of SPD matrices the usual Euclidean metric of the ambient space of symmetric matrices is derived from the Frobenius resp.Hilbert-Schmidt inner product and norm.It has long been recognised that this metric has some undesirable properties [52,1,47,2,3,15,16,20], in particular the so-called swelling, fattening, and shrinking effects in interpolation resp.averaging [54,33,55,22,21].This refers to the ellipsoid representing the SPD matrix, and is connected to the non-monotonic interpolation of the eigenvalues or some of their functions (e.g. the determinant) along Euclidean geodesics, i.e. straight lines.In the context of stochastic material modelling this means a loss of anisotropy for the whole ensemble, which may be undesirable in some situations.In the literature just cited, there is an intensive discussion-particularly coming from the medical field of diffusion tensor imaging, where these undesirable effects are not acceptable-on how to make the SPD tensors into a Riemannian manifold with different metrics [53].
From this one obtains variatonally the Fréchet or Karcher mean-the generalisation of the arithmetic resp.Euclidean mean to general metric spaces [48].This is an area of active research [61,62,63], the 'best' metric and mean seem to be application dependent, and here we shall formulate some desiderata for a metric-and thus mean-for SPD tensors of material properties, and propose one which in our view is most suited for this kind of field.
As to the modelling of random tensor fields with given invariance properties, the culmination of a long series of works [39,51,38,37,40] is reported in [41], which looks at stochastically homogeneous and isotropic random tensor fields of any order in a three dimensional domain.The results reported there allow a fine control about what invariance the mean of such a tensor field has, and what kind of invariance may be required of each realisation.Using the spectral theorem for the covariance of such homogeneous and isotropic fields, as well as the representation theory of the appropriate groups, here one may find the spectral resolution of the covariance function and from it a Karhunen-Loève like representation (e.g.[43]) of the random tensor field with the desired invariance properties.In these deliberations, no special attention is payed to the topic of positive definiteness, as the methods used are purely vector space like linear combinations or generalisations such as series and integrals.And although the limiting expression may well represent a positive definite tensor, in a practical computational situations these series and integrals have to be truncated resp.numerically approximated, and this step may lead to a numerical loss of positive definiteness.
A reduced non-parametric approach to generate random SPD matrices, which takes explicit account of the topic of positive definiteness, was presented in [57,58].Here an algebraic property, namely that positive elements in the algebra are squares of other elements in the algebra, is used to ascertain that each generated tensor field is indeed SPD, and the theory of C * -algebras [56,64] ensures that each SPD tensor(-field) has a square root.Based on this, the generation of random elasticity tensors with specified invariance in the mean and fully anisotropic invariance in each realisation, resp.controlled elasticity class invariance with constant spatial orientation of the symmetry axes in each realisation, is shown in a series of papers [23,24,25,26,27,28,49,60,59].
Another example for a frequent approach to ensure the SPD property, based on spectral and C * -algebra properties, is to use the exponential function, resp.the logarithm in the opposite direction.This is employed in e.g.[27,49,33,55,22,21], and it can even be combined with the previous squaring approach [27,49] to ensure compliance with a particular invariance class.This technique is used quite widely when describing measured diffusion and conductivity tensors, namely instead of addressing the SPD tensor itself, one focuses with the description on its logarithm-and C * -algebra theory assures us that a SPD tensor has a unique main-branch logarithm.To obtain the tensor, one finally takes the exponential, thus ensuring positive definiteness.This approach can take advantage of the fact that the invariance classes are the same for a tensor and its exponential.We shall follow a similar path here and essentially use the exponential function to ensure positive definiteness numerically in each realisation.
The plan of the paper is as follows.In Section 2 we formulate the overall problem more precisely and spell out the desiderata for the modelling and generation of symmetric positive definite (SPD) matrices, as well as for the mean, based on requirements for the metric.The desiderata for the mean or expectation will be set forth in Section 2.3.As we shall use the spectral decomposition in Section 3 to separate orientation and strength / size, and an appropriate definition of mean and corresponding modelling approach, one has also to be able to generate random orthogonal transformations.In the end the SPD matrix at each point in the domain in for each realisation will be a function of a few new numerical parameters, here the logarithms of the eigenvalues as well as a description for the orthogonal eigenvector matrix, which we shall base on the Rodrigues formula via the so called rotation vector.
It is these eigenvalue parameters as well as the rotation vector which vary randomly and in the domain.On this level, the step to spatially random fields is relatively simple, and one may borrow the idea formulated in [23,24,25,26,27,28] and regard the representation parameters as random fields without any positivity or invariance constraints.This then automatically produces random SPD tensor fields when played back via the spatially point-wise exponential map on the product Lie group, and from these to the SPD tensors.
As the spatial variation of the random tensor fields is well known from the corresponding correlation functions-no matter what kind of mean was used for the computation of the correlation function, this can be converted to a joint correlation function for the log eigenvalues and the rotation vector.Spatially homogeneous and invariant -wise isotropic tensor fields, which is the subject of [41], are a special case of this more general approach.As may be gleaned from this short overview, ensuring positive definiteness and the proper invariance at each point and realisation is a main concern, so we focus mainly on this part.Random fields of the parameters can then be generated through well known synthesis techniques connected to the Karhunen-Loève representation (e.g.[43]), or even more general representations with tight frames [35], where no nonlinear constraints like positivity have to be observed, yielding the SPD matrix through a memoryless transformation.
We also have an application in mind, which will be used in numerical examples to demonstrate the workings of the different modelling assumptions.This is the thermal conductivity of bones [65,11,12] -a highly anisotropic material -which is important in some surgical procedures.It is described in Section 4. The results of the computations are shown in Sections 4.1 and 4.2, which then displays the effect on actual computational results of different modelling assumptions about the anisotropy.While Section 5 concludes, some additional material on distance metrics for symmetric positive definite matrices is collected in Appendix A.

Problem description
Let G ⊂ R d be a bounded domain in a d-dimensional Euclidean space R d (here d = 2 or d = 3), such that the behaviour of a physical system in the domain G is governed by an abstract equation Here u ∈ U describes the state of the system lying in a Hilbert space U (for the sake of simplicity), A is a-possibly non-linear-operator modelling the physics of the system, and f ∈ U * is some external influence (action/excitation/loading).Furthermore, the coefficient C is regarded as a material tensor which represents intrinsic physics-based properties of materials, such as thermal conductivity, magnetic permeability, chemical diffusivity, and similar [50,36].It is modelled as a tensor-valued field, at each point in the domain G taking values in the open convex cone of real-valued second order positive-definite tensors: (2) As a simple example of such an equation, the reader may think of the well known thermal conduction equation -the example to be used in the numerical experiments: plus some appropriate boundary conditions on ∂G.Here the system state is described by the temperature T (x) at each point x ∈ G, the 'loading' f (x) are thermal sinks and sources, and C(x) = κ(x) ∈ Sym + (d) is the second order tensor field of thermal conductivities.

Stochastic model
Coming back to the general model Eq.( 1), due to uncertainties in the exact value of C, which may be due to a highly heterogeneous material, or simply a lack of knowledge about the exact value, we assume that a probabilistic model is adopted, so that it is modelled as a second-order tensor valued random field, i.e. a mapping on a probability space (Ω, F, P), where Ω represents the sample space containing all outcomes ω ∈ Ω, F is the σ-algebra of measurable subsets of Ω, and P is a probability measure.Similarly, one may assume the boundary or initial conditions, as well as external actions to be uncertain.However, this is secondary to the purpose of this paper and will not be covered further.To avoid unnecessary clutter in the notation the argument x ∈ G will be omitted, unless it is necessary for an understanding of the expressions.The same will often be true of the argument ω ∈ Ω.
Evidently, this change transforms Eq.( 1) to its stochastic counterpart, taking the form in which u(ω) is a stochastic response, due to the randomness in C(ω).Here we shall not be concerned on how this stochastic response may be computed (see e.g.[43]), but primarily on how to represent the symmetric positive definite (SPD) tensor C(ω) ∈ Sym + (d).

Representation requirements
One important point in the description of material properties of C(ω) ∈ Sym + (d) are symmetries in the sense of invariance against certain transformations.It is well known that the collection of all such transformations form a group G [50,41], and here we shall only be concerned with so-called point groups.The representation of any such subgroup G ⊆ O(d) of the spatial orthogonal group O(d) of R d is fairly direct, and invariance here means that for all R ∈ G one has R T C(ω)R = C(ω), a linear constraint on C(ω).
The set of tensors C ∈ Sym(d) which satisfy such an invariance are easily seen to form a subspace in Sym(d) (e.g.[41]), the invariant subspace Sym(d) G ⊆ Sym(d) of the action of G. Hence the distribution of the C(ω) lives only on this subspace, or more precisely on Sym + (d) G := Sym + (d) ∩ Sym(d) G , as there will be no realisations outside of it.For the mean C ϑ := E ϑ (C(•)), which is an averaging operation, it is possible (e.g.[41]) that it is invariant w.r.t. a larger group G m ⊇ G.
As was noted already in Section 1, it is not all obvious which kind of mean to take here, as it depends on how Sym + (d) is seen as a metric manifold with metric ϑ; more about that will follow further below in Section 2.3 and in Appendix A. This implies that as due to G m ⊇ G the G m -invariant subspace Sym(d) Gm must be a subspace of the G-invariant subspace Sym(d) G .All this leads to the following requirements for the parametrisation of C(ω): The first, third, and fourth item are not difficult to satisfy, as they are a linear constraints in the vector space of all d × d matrices, all one has to do is to pick an appropriate basis.The first item is almost trivial, and the second item is elaborated on in Section 3.1.The third and fourth point are discussed in Section 3.2, whereas the question of the mean in the fourth item will be addressed in the following Section 2.3.

Means for symmetric positive definite matrices
The second and third items on the above requirements list will be discussed in more detail later, but to discuss the fourth item, one has to first define the kind of mean or average resp.expected value one wants to use.As was already alluded to in Section 1, this is very much connected to the way in which distances on Sym + (d) are measured.And computing a weighted mean between two items is actually a way of interpolating.The next few paragraphs are meant to illustrate and remind the reader of this connection.
It has been well known since antiquity that there are different ways to compute a mean or an average, one just has to look at the classical Pythagorean means.The arithmetic or Euclidean mean x of a set of vectors x 1 , . . ., x n ∈ R k -best known by almost anyone is the case k = 1 -with standard Euclidean squared norm x 2 = x T x and corresponding Euclidean distance ϑ E (y, x) = y − x , is given by which displays right away the variational characterisation of the mean x by elementary calculus.The minimal value Ψ (x) of the function Ψ is called the variance of the set 2 in Eq.( 6) are the squares of the Euclidean distances ϑ E (x j , x) between x j and x.One may also note that for n = 2 one gets from Eq.( 6) that x = is the half-way interpolation on the 'shortest path'-the straight line-between x 1 and x 2 , which is an elementary example of the connection between averaging and interpolation.The Eq.( 6) is easily generalised to not necessarily equal (probability) weights w j > 0, n j=1 w j = 1, instead of the constant w j ≡ 1/n, and vectors x j ∈ V in any Euclidean vector space V with inner product x, y V and corresponding norm, giving xw = n j=1 w j x j .Hereby it was assumed that-as is the case for an isotropic Euclidean distance-the gradient ∇ x x, x V of the inner product is 2x.
The generalisation of Eq.( 6) corresponding to a general probability distribution P of a V-valued random variable Y would be and again the minimum Ψ P (Y ) is the variance of the distribution P.
In case the underlying set does not have a vector space structure but only that of a metric space, i.e. when only a distance metric ϑ is defined, or, like in our case, another distance than the Euclidean is to be used, one may define the so-called Fréchet mean or average attached to that metric [48,53] through a variational characterisation similar to Eqs.( 6) and (7), by replacing the Euclidean metric ϑ E in those expressions by ϑ.Thus one may define the Fréchet-Karcher mean on such a metric space S with a distance function ϑ which is at the same time a measure space with a general probability measure P for an S-valued random variable R as and again the variance is defined by the minimum value Ψ P (R).If a unique minimiser exists, for general metric spaces this is called the Fréchet mean and variance, and for Riemannian manifolds a local minimiser of the variance functionals in Eq.( 8) is also sometimes called the Karcher mean [48,53].These foregoing observations show the intimate connection between distance metric, mean or average and expectation, and interpolation.The properties of the Fréchet mean or average depend completely on the underlying metric, see Eqs.( 6) and (8).Turning now to positive definite tensors C ∈ Sym + (d), as mentioned already in Section 1, the Euclidean or Frobenius resp.Hilbert-Schmidt norm Eq.( 47) in Appendix A) has long been recognised to result in some undesirable properties [52,1,47,2,3,15,16,20], especially the swelling, fattening, and shrinking effects in interpolation resp.averaging [54,33,55,22,21].It is therefore of considerable interest to consider different metrics on Sym + (d) [53,48].In order not to distract the reader's attention too much by reviewing different metrics used on Sym + (d), we relegate this material to Appendix A. Having chosen some distance metric ϑ D with desirable properties on Sym + (d), one defines the B-based Fréchet variance analogous to Eqs.( 6) and ( 8): From this one obtains the definition of the corresponding Fréchet or Karcher mean resp.expectation Inserting this minimiser-if it exists-into the functional Eq.( 9), one obtains the variance of the Sym + (d)-valued random variable C(ω) as the minimum value Let us collect here some, in our view desirable, properties of a possible metric for material tensors C(ω) ∈ Sym + (d), which then may point to a certain kind of mean to use.As an extreme case one could argue that the mean C ϑ D should be such, that if used in Eq.( 1), one obtains the mean ū of the stochastic solution of Eq.( 5): This is in fact the goal of homogenisation techniques, see e.g.[34,31], but it is probably too much to require from a general kind of average, as it would depend on the model Eqs.( 1) and ( 5) and thus be highly problem dependent.Hence we try and settle with some more general requirements based on invariance considerations.The invariance requirements for the metric translate immediately into invariance requirements for the corresponding Fréchet mean.
From the point of physical modelling, if C(ω) appears in a constitutive model, there is typically an inverse formulation involving C −1 (ω).This means that one normally knows as much about the distribution of C(ω) as about the distribution of C −1 (ω), as the set Sym + (d) is stable under inversion.So one property of a desired metric ϑ D on Sym + (d) would be invariance under inversion.Changing physical units would correspond to a scaling by a positive factor, so the desired metric should be invariant under scaling.And as a change of the coordinate frame by an orthogonal transformation R ∈ O(d) changes the constitutive tensor C to R T CR, one would like the desired metric ϑ D to be invariant under such transformations in order to respect the principle of frame indifference.Collecting everything, the desiderata for a metric ϑ D , and hence for a Fréchet mean based on it via Eq.( 10), then are to satisfy for all C 1 , C 2 ∈ Sym + (d), α > 0, and all R ∈ O(d): • invariance under inversion: • invariance under frame transformations: By choosing a metric ϑ D on Sym + (d) which satisfies those requirements, the Fréchet-Karcher mean E ϑ D [C] based on that metric will automatically inherit all those properties, as can be easily gleaned from the very definition above in Eqs.( 9) and (10).
The usual Frobenius resp.Euclidean distance which leads to the normal arithmetic mean Eq.( 7) satisfies the last requirement, but not the first two.Additionally, the swelling, fattening, and shrinking effects that result from the straight line interpolation or averaging in the vector space Sym(d) of two elements of Sym + (d) are highly undesirable [2,54,33,55,22,21].Any Fréchet-Karcher mean based on a distance function which exhibits these undesirable effects in interpolation would also exhibit these effects in the averaging procedure, i.e. the thus defined mean would be more swollen, fatter, or shrunk than the population it is derived from.This short consideration shows that the normal arithmetic or Euclidean mean based on the Frobenius distance ϑ F (Eq.( 48) in Appendix A) is not ideal for SPD material tensors.
What comes out of from the more detailed considerations in Appendix A regarding possible metrics on Sym + (d), is that all the other distance metrics considered there, namely the geometric or affine-invariant distance ϑ G (Eq.( 49) in Appendix A), the logarithmic or log-Euclidean distance ϑ L (Eq.( 52) in Appendix A), and the scaling-rotation distance-which we will often refer to as just the scaling distanceϑ S (Eq.( 55) in Appendix A), and which in turn is based on the product Lie group distance ϑ P (Eq.( 54) in Appendix A), all satisfy all of the above three requirements.This means that for a choice between them additional considerations are necessary.These are how to ensure the requirement listed above in Section 2.2, namely how to ensure positive definiteness through an appropriate representation, and at the same time ensure the symmetry class of the corresponding mean and each single realisation.We shall employ the afore mentioned scaling distanceϑ S (Eq.( 55) in Appendix A), based on the product Lie group distance ϑ P (Eq.( 54) in Appendix A) -as it allows the desired fine control of size / strength and orientation, as well as making it easy to ensure positive definiteness.

Stochastic modelling
Methods to ensure positive definiteness will be briefly considered in Section 3.1.We shall settle for an approach based on the exponential map, which is most easily handled in the spectral decomposition representation.This kind of representation is also well suited to address the problem of symmetry invariance, which will be done in Section 3.2.The outcome of this approach is that the tensor C(ω) at each spatial location may be described by a number-depending on the symmetry class-of parameters, which can be clearly separated into scaling or 'size' parameters, and others which determine the orientation of the symmetry axes.Although this is not the main topic of this paper, in Section 3.3 it is sketched how, borrowing an idea from [23,24,25,26,27,28], these parameters can be modelled as random fields without any positivity or invariance constraints, thus giving a random tensor field C(x, ω) through a memoryless point-wise transformation.

Ensuring positive definiteness
The second item on the list of requirements for representation above in Section 2.2, that each C(ω) is SPD, can pose some problems, as Sym + (d) ⊂ Sym(d) is not a subspace, but geometrically an open convex cone.Therefore we shall look only at representations which can ensure that the resulting C(ω) is SPD, even under numerical approximations.Some of these representations are actually closely connected to some of the metrics described in Appendix A; and although the kind of metric and associated mean and the representation are formally independent, some representations make it easier to compute certain means.
The first representation we shall look at is the exponential map, e.g.see [2,3,27,49,33,55,22,21] for some recent examples, which is analytic and one-to-one and onto, i.e. it can be taken as a global chart.For a H ∈ Sym(d) we set C = exp(H) ∈ Sym + (d).This can be computed either via spectral calculus, or simply through the everywhere convergent power series.One may also note that thanks to the Cayley-Hamilton theorem the series can actually be written as a finite sum.This representation also makes it very easy to compute the logarithmic or log-Euclidean distance ϑ L , Eq.( 52) in Appendix A.
The main result is that whatever errors may be committed by approximating H(ω), by computing the exponential exp(H(ω)) in the end the result is guaranteed to be SPD.One advantage of the representation is that Sym(d) is a linear vector space, and H ∈ Sym(d) is completely free in that vector space, which makes stochastic modelling a lot easier.One should also note that ensuring any kind of group invariance, the subject of the next Section 3.2, is very easy, as C and its logarithm H have exactly the same group invariance.For a fixed orientation, this then manifests itself as requiring H to be in some linear subspace of Sym(d).
One may note that by wanting to model H = log(C), the tensor C has physical units, so that it does not directly make sense to take its logarithm.This is usually dealt with by some kind of scaling with a representative inverse tensor C −1 which has the right invariance properties-this could be the mean C ϑ one wants to use-as a scaling, and using its spectral decomposition (cf.also Eq.( 13) so one may take its logarithm H(ω) = log( C(ω)).Usually it is tacitly assumed that such a transformation or scaling to a dimensionless form was performed before computing functions of C.
The second kind of representation we shall consider, is the approach (see e.g.[57,58,23,24]) where the problem to ensure that each C(ω) is SPD is solved by representing the tensor as a generalised square.The very definition of being positive in an algebra is that the element is a generalised square, i.e.C(ω) = G T (ω)G(ω) with G(ω) ∈ GL(d) nonsingular.In this approach as described in the cited literature, the matrix G(ω) is chosen as upper triangular, and thus one is certain to use only the right number of parameters.In [27,28] this representation as a generalised square is combined with the exponential representation described above.Now GL(d) is not a vector space, and any L ∈ GL(d) has to satisfy the non-linear constraint det L = 0.But by controlling the diagonal elements of the upper triangular matrix G(ω) ∈ GL(d) such that det G(ω) > 0 -the matrix G(ω) may thus be seen as a Cholesky-like factor of C(ω) -one can indeed make sure that it is invertible, and this otherwise troublesome nonlinear constraint of non-vanishing determinant becomes simpler.
The third kind of representation, and this will be the one we will eventually use in our computations, allows a fine control of orientation and 'size'.This is achieved by looking at the spectral decomposition of a C ∈ Sym + (d) and its logarithm H = log(C): with Q ∈ SO(d) and the diagonal positive definite matrix Here the eigenvector matrix Q controls the orientation, whereas the diagonal matrix Λ (resp.Y = diag(y (i) ) = diag(log(λ (i) )) = log(Λ)) controls the scale, size, or strength of the tensor C.
In Appendix A the spectral decomposition Eq.( 13) is also used to investigate the scaling-rotation distance ϑ S on Sym + (d).Informally, Eq.( 13) says that it induces the representation One should pointed out though [54,33,55,22,21] that the representation map Eq.( 14) is onto, but not necessarily one-to-one; an example where this can be seen is the occurrence of multiple eigenvalues.
While it was pointed out already that Sym + (d) is not a Lie group under matrix multiplication, the subset Diag + (d) ⊂ Sym + (d) of positive definite diagonal matrices obviously is a commutative Lie group with matrix multiplication.And as SO(d) is well known to be a Lie group, the product on the left hand side of Eq.( 14) is a product of Lie groups, and thus is again a Lie group.And the advantage of having a representation on a Lie group is that the Lie group can itself be easily represented on its Lie algebrathe tangent space at the identity and thus importantly an unconstrained linear vector space-via the exponential map.In the vicinity of the origin the exponential map is an analytic diffeomorphism, and for the Lie group Diag + (d) it is actually a globally valid chart, similarly as for Sym + (d) in C = exp(H) before.
The Lie algebra of the commutative Lie group Diag + (d) ⊂ Sym + (d) is easily seen to be the subspace diag(d) ⊂ Sym(d) of diagonal (and thus symmetric) matrices.The Lie algebra of the special orthogonal group SO(d) is known to be the Lie algebra of skew symmetric matrices so(d) (W = −W T , see the following Section 3.2), and thus the spectral synthesis representation Eq.( 14) can be extended to the left with the exponential map, giving a two step modelling process: Obviously Λ = diag(λ (i) ) = diag(exp y (i) ) = Y ensures that C is SPD, and the exponential of a diagonal matrix is easy to compute.The two steps of exponential map and spectral synthesis can actually be interwoven, explicitly showing that the exponential map is used in the last step: The important advantage we see in this last representation, vis-a-vis the first two considered before, is that it separates out the eigenvalues λ (i) resp.their logarithms y (i) = log λ (i) -which physically may be conductivities or resistances-from the anglelike parameters defining W ∈ so(d), which one may use for fine-tuning in the stochastic modelling.

Group invariance
As was already indicated earlier, and as easily seen from Eq.( 13), the group invariance properties of of a symmetric tensor H and its exponential C = exp H are the same, a fact which was also used in [27,28].This means that all the invariance constraints in the Lie group resp.Lie algebra representation Eq.( 15) to be proposed here may be applied on the Lie algebra level.Thus, building on the spectral decomposition Eq.( 13), one uses the inverse of the scaling equation Eq.( 12) above together with the spectral decomposition Eq.( 13).
As regards the symmetry groups, it was mentioned already that for second order symmetric tensors in 2D or 3D the situation is relatively simple [40,41], as by reduction through the symmetry group [18] the matrix representation can be fully diagonalised.Thus in 3D there are only three symmetry classes [41] depending on the number of distinct eigenvalues, the symmetry subgroup relations are Here "orthotropic: 3 dEVs" means that there are 3 distinct eigenvalues, and the symmetry group for the orthotropic class is a subgroup of the symmetry group for the plan isotropic class (which has only 2 distinct eigenvalues), and similarly for the other cases.This means that it should be possible e.g. that each realisation is orthotropic, but the mean is isotropic, or any other such combination up the inclusion chain in Eq.( 17), such that the mean has the same or a larger symmetry group.The case in 2D is even simpler, the symmetry subgroup relations are To show the different possibilities, this will be done in three steps: first the spatial orientation will be kept fixed and only the eigenvalues (i.e. the scaling) will be random.In the case of second order tensors which is now considered here, this determines the group invariance class by the number of coinciding eigenvalues.In a second step, we shall keep the eigenvalues fixed but allow the spatial orientation of the symmetry axes to vary randomly; and in the last step both of those sources of randomness will be modelled together.

Fixed spatial orientation:
The modelling approach is easiest to see by modelling only the anisotropic scaling as random, and keep a fixed spatial orientation of the anisotropy.
With this scaling Λ, which together with the fixed orientation Q and the random element Λ(ω) = exp( Y (ω)) determine the symmetry class for each realisation corresponding to the group G, it can be arranged that for both for the log-Euclidean resp.logarithmic mean based on the log-Euclidean distance metric ϑ L (Eq.( 52) in Appendix A), as well as for the scaling mean based on the scalingrotation resp.scaling distance ϑ S (Eq.( 55) in Appendix A), which in turn is based on the product Lie group distance ϑ P (Eq.( 54) in Appendix A).Hence one has from Eq.( 19) (20) which one can arrange to have the invariance properties of the group G m ⊇ G.
The representation equation Eq.( 19) may just as well be written for H(x, ω): as diagonal matrices commute, with Y (x) = log Λ(x).Thus one establishes a link to the exponential representation.From the last expression in Eq.( 21) one may observe that in the unconstrained linear vector space Sym(d) of the log the scaling is merely a shift of origin, and as the non-dimensional form in Eq.( 12) typically just makes the formulas look more complicated, we shall drop it in the notation and simply assume that this has been carried out in the background.The Lie algebra representation Eqs.( 15) and ( 16) thus can be simplified from Eq. (19) to where an index "s" has been attached to signal that only random anisotropic scaling, i.e. a random Y s was used.Additionally, one has to require that -here we are using the Euclidean mean in the unconstrained vector space Sym(d), which corresponds to the log-Euclidean and scaling means of Λ s (x, ω), and thus Obviously there is an analogous equation to Eq.( 22) involving the reference tensor C(x), i.e. from Eq.( 20) one has: Let us first consider the simpler case of 2D in Eq.( 18), where the tensor can be isotropic with two equal eigenvalues y (1)   s = y (2)   s -a multiple of the identity-or orthotropic with distinct eigenvalues y (1)   s = y (2)  s , in which case the orthogonal matrix Q s determines the directions of the major axes of the tensor.If in the first case one demands that each realisation is isotropic-and obviously also the reference C used in Eq.( 12) -then only one random variable is needed to model both eigenvalues.On the other hand, if only the mean has to be isotropic but individual realisations may be orthotropic, then two identically distributed random variables could be used for the eigenvalues.Turning now to the case where the reference C is orthotropic, obviously one eigenvalue logarithm has to be larger than the other, say ȳ (1) s > ȳ(2) s .For the individual realisations, one may now choose to require that this relation holds also for each realisation, or only in the mean, together with some possible correlation of the eigenvalues.In case one wants to require for each realisation y (1)   s > y (2)  s , one trick to achieve this is to take two unconstrained random variables, say ξ 1 and ξ 2 , and set y (2)  s = ξ 1 and y (1)  s = ξ 1 + exp(ξ 2 ).Hence, in case individual realisations may be orthotropic, two log-stiffness like parameters y (1)  s , y (2)   s are needed, and in case the individual realisations are always isotropic, just one log-stiffness like parameters y (1)  s is sufficient.After this somewhat broad discussion of the 2D case, the 3D in Eq.( 17) case may be treated in a similar way.The three invariance classes are determined again by the multiplicity of the eigenvalues, and the three cases are isotropy-all eigenvalues equal y (1)  s = y (2)  s = y (3)  s -or plan isotropy, where two eigenvalues are equal (defining the plane of isotropy) and the third eigenvalue differs from these, or finally orthotropy, where all three eigenvalues are different.So it is possible to require that each realisation is planisotropic (or orthotropic), but that the mean is isotropic, or that each realisation may be orthotropic but the mean plan isotropic (or isotropic).Thus, depending on which case one wants to model, one to three log-stiffness like parameters y (j)  s are needed.
Random spatial orientation with fixed scaling: Nevertheless, it is a well known fact that anisotropic materials often also exhibit uncertainty in orientations due to the randomly oriented fibres/grains/crystals.The general procedure for this was pointed out before in Eq.( 31).The directional uncertainty in H(ω)-by keeping the material symmetry constant-is incorporated by subjecting the eigenvectors in Q in as in Eq.( 31) to random rotations.So for the start we keep the eigenvalues constant, and only subject the eigenvectors to random fluctuations; we use a subscript "r" to signify that only the eigenvectors are random.The stochastic tensor H r (ω r ) ∈ Sym(d)-with random orientations only-is then analogously to Eq.( 31) defined in the form: where R(ω r ) ∈ SO(d) is a random rotation matrix chosen such that each realisation of H r (ω r ) resp.C r (ω r ) satisfies the appropriate invariance requirements.We also want which requires It now becomes necessary to represent random rotations like R(ω r ) ∈ SO(d) in Eqs.( 24) and (25).There are very many approaches to represent rotations, as well as a large amount of easily accessible literature; this will not be reviewed presently as this would lead the discussion astray.In order to arrive at a vector space setting preferred by us, we resort to the well known correspondence Eq.( 27) between the Lie group SO(d) of orthogonal matrices with unit determinant and its Lie algebra so(d) of skew-symmetric matrices [45,46,8], also alluded to in Appendix A. R = exp(W ) ∈ SO(d), and W = log(R) ∈ so(d). ( As one wants that where the Lie algebra so(d) is an unconstrained linear vector space.Any distribution of W (x, ω) symmetric to the origin 0 ∈ so(d) will produce the desired result.Of course, it will not be possible to establish a one-to-one relation between the Lie group SO(d) and some real vector space, as SO(d) is compact-this is easily seen in 2D where SO(2) is homeomorphic to the one one-sphere, i.e. the unit circle in R 2 , or in 3D where the representation through unit quaternions shows that SO(3) is homeomorphic to the three-sphere, i.e. the unit sphere in R 4 -and a real vector space is not compact.But the relations in Eq.( 27) come pretty close to this desideratum, which makes the probabilistic modelling quite a bit easier.
Especially in 3D one may use an explicit version of Eq.( 27), the well-known Rodrigues rotation formula, e.g.[8].There is a familiar correspondence between a skew-symmetric matrix W ∈ so(3) and the defining Euler vector w = [w 1 , w 2 , w 3 ] T ∈ R 3 along its unit rotation axis vector w/ w which is the unique rotation axis of a non-trivial rotation R with rotation angle φ = w .Rodrigues's rotation formula may then, thanks to the Cayley-Hamilton theorem, be written as and inversely [17] The defining vector w(ω r ) is a R 3 -valued random vector with values in the ball with radius π, and the random angle of rotation is φ(ω r ) = w(ω r ) , defined on a probability space (Ω r , F r , P r ).But in case each realisation is plan isotropic (this assumed to be the 1-2 plane), there is no point in rotations with a w 3 component -a rotation in the 1-2 plane.In this instance one only needs two angle-like parameters w 1 , w 2 .Note that in the 2D case the axis of rotation w = (0, 0, φ) is fixed perpendicular to the plane; only the rotational angle φ(ω r ) resp. the parameter w 3 is modelled as a circular/angular random variable [42,32].And in case each realisation is isotropic, no rotation at all is required, as it will have no effect.
In any case one obtains a random R(ω r ) ∈ SO(d), determined by 0-3 angle-like parameters w j , defining the axis of rotation w.Again, as long as the distribution of w(ω r ) resp.φ(ω r ) is symmetric around the origin, the expected value of R(ω r ) of the here adequate Fréchet or Karcher mean with the metric ϑ R (defined in Appendix A) is the identity, I = E ϑ R (R).Hence Eq.( 26) holds in this case.

Random spatial orientation and scaling:
To account for both random orientation and scaling attributes in the tensor H, one may combine the two approaches.If these two aspects are considered independent, one may define the probability space as a product space (Ω, F, P) with Ω := Ω s ×Ω r , F := F s ⊗F r , and a product measure P := P s P r .In any case, denoting a combined rotational-scaling uncertainty with a subscript "rs", and reusing the already defined subscripts "s" and "r", let C s (x, ω s ) denote a representation of an SPD tensor field with given group invariance and a fixed orientation of anisotropy axes, and E ϑ L [C s (x, •)] = C(x), like the one just mentioned in Eq.( 22), and let R(x, ω r ) ∈ SO(d) be a random orientation as in Eqs.( 24)- (27).
is then a representation of an SPD tensor field C(x, ω rs ) with given group invariance of C s (x, ω s ) and random anistropic scaling, but in addition with a random orientation of anisotropy axes-with the random orientation described by the SO(d)-valued random field R(x, ω r ).One then has for the scaling mean, based on the scaling distance ϑ S (Eq.( 55) in Appendix A):
The joint correlation structure of the random parameters z(x, ω) is described through a joint correlation matrix J ∈ R m×m : Aiming at a Karhunen-Loève expansion of a corresponding spatial field (see e.g.[43]), one has to solve the eigenproblem for eigenvalues µ and eigenvector fields η (x) ∈ R m .It should be noted that spatial homogeneity would be evident as usual from the fact that J (x 1 , x 2 ) would be a function of only the difference of the arguments, x 1 − x 2 .In that case the eigenproblem Eq.( 34) is a convolution equation, and it is well known how to solve those by Fourier methods, as indeed the Fourier functions are (possibly generalised) eigenfunction of the operator J represented by Eq.( 34).In any case, the orthonormal eigenfunctions η (x) determine uncorrelated mean zero unit variance random variables ξ (ω) through the projection of z(x, ω): But normally we want to perform the reverse task, a spectral synthesis, and thus the random variables ξ (ω) have to be chosen such that z(x, ω) has the proper statistics at each x ∈ G.This problem is outside the scope of this work and has been addressed elsewhere, see e.g.[43] and the literature therein.
Assuming the proper mean zero unit variance random variables ξ (ω) have been chosen, for the spectral synthesis resp.Karhunen-Loève expansion one typically has an infinite series-or an integral in case of a continuous spectrum-over the spectrum µ ∈ σ(J) of the integral covariance operator in Eq. (34) (see e.g.[43]): It is of course also possible to use more general expansion functions η , e.g.tight frames [35], an approach which may have advantages especially in the lognormal case [4,5].Note that through the use of a Lie group mapping onto Sym + (d), we are able to essentially do all the modelling in the Lie algebra, which is a free vector space, and the expression Eq.( 36) is an example of this.This brings the modelling of random spatial fields or temporal stochastic processes of Sym + (d) tensors to the same level of difficulty as the usual modelling of random fields resp.stochastic processes.

Application: heat conduction in a proximal femur
The application we have in mind is the steady-state heat conduction in bones, a highly anisotropic material.More specifically, we shall look at both a 2D and 3D model of the proximal femur, this is the upper part of the thigh bone close to the hip.Here we shall investigate the effects of different models of the anisotropy of the heat conduction tensor on the heat conduction and temperature distribution in the proximal femur.
Referring to the deterministic model of heat conduction defined in Eq.( 3), let the conductivity tensor κ ∈ Sym + (d) be modelled as a random tensor as described in Eq. (31).The scaling terms Λ s = diag(λ (i) s ) ∈ Sym + (d) of the tensor κ are modelled as positive log-normal random variables: i.e. a non-linear transformation of the Gaussian random variables where the parameters {µ i } d i=1 ∈ R and {σ i } d i=1 ∈ R + denote the mean and standard deviation corresponding to the Gaussian distribution.
Further, the directional uncertainty of κ is accounted by modelling the eigenvectors in Q r as von Mises-Fisher distributions defined on the unit sphere S d−1 [42].If v is one of the eigenvectors of the matrix Q r , then a von Mises-Fisher unit vector v(ω r ) with mean direction µ circ ∈ S d−1 has a probability density function (PDF) of the form: where η ∈ R + is the measure of concentration, and I p is the modified Bessel function of order p.When d = 2, as only the rotational angle φ(ω r ) is modelled as a random variable on the unit circle S 1 , the above PDF reduces to the von Mises distribution [6]: in which φ is the mean direction of the von Mises circular random variable φ(ω r ).
By rewriting the deterministic model in Eq.( 3) to a stochastic one which has to be understood in a weak sense, both as regards the spatial as well as the stochastic dependence.The goal is to determine the random temperature field T (x, ω) : G × Ω → R, assuming deterministic boundary conditions and heat source.
Given the weak form of Eq.( 41), we use the finite element method to seek for an approximate solution T h (x, ω) on a finite-dimensional subspace U h ⊂ U of the solution space U of Eq.( 3), where as usual h is the discretisation parameter-an indicator of element size.Subsequently, the objective of this study is to determine the statistics like the mean and standard deviation of the solution T h (x, ω) using the classical Monte Carlo method (MC) [44,19].Given, i.i.d.{κ(ω i )} N i=1 samples, where N is the sample size, one may compute the statistics of discretised solution T h,N := {T h (x, ω i )} N i=1 .The unbiased sample mean [13] is then calculated as and the corrected sample standard deviation as Along with the temperature field, we also evaluate the heat flux vector field q(x, ω) : G × Ω → R d , whose approximate solution q h (x, ω) is determined on a finite-dimensional subspace V h ⊂ V, where V is the solution space of heat flux, q(x) = κ∇T (x).Further, we determine the second order statistics of the magnitude of approximate heat flux vector field in the Euclidean norm i.e. q (t) h (x, ω) = q h (x, ω) .More specifically, the statistics such as the mean and standard deviation are obtained by sampling based MC estimators µ (MC) (q (t) h,N ) and Std (MC) (q (t) h,N ), respectively.Apart from evaluating the statistics of magnitude of heat flux field q h (x, ω), one may also be interested in understanding the statistics of direction of heat flow.To this, one may compute the second order circular statistics of the normalized heat flux, qh (x, ω) = q h (x, ω)/q (t) h (x, ω).Subsequently, the MC mean estimate of sampled unit vector field qh,N := {q h (x, ω i )} N i=1 can be estimated accordingly [42,32]: where L = µ (MC) (q h,N ) ≤ 1 represents resultant Euclidean length of the sample mean µ (MC) (q h,N ) ∈ R d , and µ (MC) circ (q h,N ) = µ (MC) (q h,N )/ µ (MC) (q h,N ) is a unit vector defining the sample mean direction of qh,N .With the help of the former one, the sample circular standard deviation of qh,N can further be evaluated as A two-dimensional proximal femur with a body of approximately 7 cm in width and 21.7 cm in height is considered, see Fig. 1.A uniformly distributed surface heat flux with a resultant heating of 0.1 W and a fixed temperature of 0 • C are applied at locations as shown in Fig.

Results on 2D proximal femur
1 [11].The computational domain is discretized by the finite element discretization method with 171 four-noded plane-stress elements, giving a mesh with 206 nodes and 390 degrees of freedom.In this study, each deterministic simulation obtained by sampling the stochastic parametric space is solved by the Calculix solver [14].For the MC simulation, the Pearson coefficient of dispersion δ = 0.1 (for the scaling uncertainty) and N = 10 5 samples are considered for all the cases in this study.
For simplification, we use abbreviations defined in Table 1 which are used to refer different modelling scenarios considered in this study.Each model will have a unique name in the form: reference tensor symmetry-realisation symmetry-type of model.For instance, "iso-ortho-scl" represents a model with random scaling only, such that the reference tensor belongs to isotropic symmetry, whereas the realisations are of orthotropic symmetry.Three such possible scenarios are explored in this article-"iso-iso-scl" (isotropy-isotropy-random scaling), "iso-ortho-scl" (isotropy-orthotropy-random scaling), and "ortho-ortho-dir" (orthotropy-orthotropy-random direction).Consequently, using the Mean/realisation symmetry Type of model isotropy orthotropy random scaling random direction iso ortho scl dir Table 1: Model abbreviations subscripts defined in Section 3.2 for random scaling and random orientation, two reference conductivity tensorsκ iso s ∈ Sym + (2) belonging to the isotropic symmetry class (denoted by superscript iso) and κ ortho r ∈ Sym + (2) with orthotropic symmetry (with superscript ortho)-are considered, whose values are tabulated in Table 2 [12].Furthermore, we simulate the grain orientations of the femur in a spatially homogeneous sense.As a result, the eigenvectors of the tensor κ ortho r corresponding to the eigenvalues λ (1)   r and λ (2)   r are oriented at an angle of π/4 in a clockwise direction to the x and y-axis, respectively.
Table 2: Reference conductivity tensors with corresponding eigenvalues of 2D femur Random scaling with fixed symmetry: the first modelling scenario is the one in which we perform random scaling only, such that the reference tensor κ iso s , as well as the realisations, belong to isotropic symmetry (denoted by iso-iso-scl); see Fig. 2 for a schematic overview of the model.In such a case, the random tensor κ s (ω s ) ∈ Sym + (2) is modelled-as in Eq.( 22) -with random scaling only.Here the uncertain scaling parameters λ (1)  s (ω s ) = λ (2)  s (ω s ) are modelled as identical log-normal random variables.Their corresponding PDFs are shown in Fig. 2

(a). A geometrically equivalent (circular) representation of reference tensor and realisations {κ
, where the radius of 20 circle is defined by 1/ λ (1) s (ω s ) and 1/ λ (2) s (ω s ) [36], is presented in Fig. 2(b).As can be seen the preservation of isotropic symmetry is apparent given the constant circular shape among the random variable realisations.However, the varying size of circles signifies a variation in scaling parameters.
Figure 2: Stochastic tensor κ s (ω s ) with fixed symmetry (iso-iso-scl): (a) log-normal PDF of identical random scaling values λ (1)  s (ω s ) = λ (2)  s (ω s ), (b) geometric visualization of reference tensor κ iso s and realisations of random tensor κ s (ω s ); the geometries are enhanced by scaling the eigenvalues by a factor of 0.5 Random scaling with varying symmetry: in the next scenario, we consider a similar stochastic tensor κ s (ω s ) ∈ Sym + (2) as described previously with a similar reference tensor κ iso s .The difference being that the eigenvalues {λ (i) s (ω s )} 2 i=1 are modelled as independent identically distributed (i.i.d.) log-normal random variables.Due to which the realisations of model κ s (ω s ) belong to orthotropic symmetry, thus, this model is denoted as iso-orthoscl.Furthermore, in orthotropic realisations, the eigenvectors of κ s (ω s ) related to random eigenvalues λ (1)  s (ω s ) and λ (2)  s (ω s ) are constrained at an angle of π/4 in clockwise direction with respect to x and y-axis respectively.The schematic representation of this model with varying symmetry is displayed in Fig. 3, in which the corresponding PDFs are shown in Fig. 3(a).Closer inspection of the Fig. 3(b) shows the shift from circular geometry of isotropic reference tensor to orthotropic elliptical shape in realisations.One may also notice the variation in size of ellipses in the realisations as the semi-axes of ellipse are determined by 1/ λ  Here the model is designated as ortho-ortho-dir, since, a fixed orthotropic symmetry is preserved in mean and realisations.The random angle of rotation φ(ω r ) is modelled as a von Mises random variable (see Eq.( 40)) with zero mean φ = 0 with respect to x-axis and concentration parameter η = 75.In Fig. 4(a), the PDF f (φ(ω r )) is shown on a real line over the domain (−π, π], whereas the realisations of random variable φ(ω r ) and its sample mean vector (solid straight line) are plotted on a unit circle in Fig. 4(b).In this study, we model directional uncertainty on the x − y plane with z as the axis of rotation.Thereby, the corresponding random rotation matrix R(ω r ) in Eq.( 24) takes the form: Here we consider a right-handed Cartesian coordinate system, meaning that positive values of random variable φ(ω r ) indicate counterclockwise rotation and vice versa.Fig. 4(c) depicts the corresponding elliptical geometric visualization of the model.As can be seen that the orientation of elliptical semi-axes fluctuates in the realisations {κ r (ω r i )} N i=1 ≡ {κ r i } N i=1 , whereas its size and shape remain constant in the mean and realisations.
Uncertainty quantification: Fig. 5 shows the MC mean estimate of the nodal temperature (NT) T h,N for the three described stochastic models.For easier interpretation, in  this study, the results are displayed on a uniform scale.Hence, the values range from zero to a maximum value (corresponding to a given estimate of the desired output quantity).Clearly, the results in Figs.5(a) and 5(b) are similar, as the considered reference tensor κ iso s is identical in both the cases.Due to the assumption of an orthotropic mean tensor κ ortho r , in the third figure Fig. 5(c), we see a significant difference in the magnitude/contour pattern of the temperature.
Fig. 6 further depicts the standard deviation of temperature T h,N .Due to different modelling assumptions, it is apparent that all three results in this figure are different from each other.But Fig. 6(c) stands out the most, where the stochastic influence on nodal temperature is much lower as compared to the other two cases.
In the first two scenarios, only the scaling parameter is uncertain, whereas in the third example we consider only directional uncertainty.We know that under the assumption of deterministic boundary conditions, the random temperature field T (x, ω) varies inversely to the stochastic coefficeint κ r (ω r ).As the scaling parameters of the model κ r (ω r ) remain constant, the impact on the variation of temperature field T (x, ω) due to directional randomness is small.On the other hand, when the scaling parameters are modelled as varying, the standard deviation estimate of temperature T h,N becomes more sensitive as evident in Figs.6(a) and 6(b).Interestingly, when Fig. 6(b) is compared to Fig. 6(a), it is clear that, varying the material symmetry from higher (isotropy) to lower (orthotropy)  (a) iso-iso-scl (b) iso-ortho-scl (c) ortho-ortho-dir Furthermore, the sample mean estimates of total heat flux (THFL) q (t) h,N are presented in Fig. 7, where, Figs.7(a) and 7(b) showcase similar results as described previously (see Figs.5(a) and 5(b)).However, in comparison to these two figures, a visible difference in maximum total heat flux and contour pattern noticed in Fig. 7(c).In Fig. 8, the corresponding estimated standard deviation of THFL q (t) h,N is plotted.Here, the most interesting aspect is seen in Fig. 8(a), where, the standard deviation is close to zero.It turns out that, in the model iso-iso-scl, the stochastic coefficient κ s (ω s ) has almost perfect inverse correlation with temperature gradient field ∇T (x, ω).Thus, the randomness in the model κ s (ω s ) has almost no stochastic impact on THFL q (t) h,N .However, on the contrary, the estimate in Figs.8(b) and 8(c) is visible, signifying the stochastic influence of the considered input models on THFL q (t) h,N .One may conclude that directivity has more impact on the heat flux than the scaling parameter.Additionally, the directional mean µ M C circ (q h,N ) and standard deviation Std (MC) circ (q h,N ) of normalized heat flux (NHFL) qh,N -from Eq.( 44)-are displayed in Fig. 9.The mean orientation (vector representation) in Figs.9(a) and 9(b) look similar, however, with a closer inspection of Fig. 9(c), the difference is apparent.Similar to the results in Fig. 8, the circular standard deviation estimate in Fig. 9(a) is also close to zero, showing once again the insensitivity of directional attribute of NHFL qh,N to scaling uncertainty present in the model κ s (ω s ).Also, in comparison, the standard deviation estimate in Fig. 9(c) is more significant than in Fig. 9(b).

Results on 3D proximal femur
A three-dimensional proximal femur configuration of dimensions 45 mm in width and 154 mm in height is considered.The boundary conditions with identical values as used in the 2D example are imposed, shown in Fig. 10.The computational domain is discretized by a four-noded tetrahedral finite element mesh, comprising 12504 nodes, 3166 elements and 35115 degrees of freedom.We use similar modelling scenarios as in 2D case i.e. iso-iso-scl, iso-ortho-scl and ortho-ortho-dir (see Table 1 for model abbreviations).The considered isotropic κ iso s ∈ Sym + (3) and orthotropic κ ortho r ∈ Sym + (3) reference conductivity tensors are tabulated in Table 3.Furthermore, by fixing the negative y-axis as the rotational axis, the directional vectors of the model κ ortho r corresponding to the eigenvalues λ (1)  r and (a) iso-iso-scl (b) iso-ortho-scl (c) ortho-ortho-dir Eigenvalues (W/mK) λ (1)  s = 0.54, λ (2)  s = 0.54 λ (3)  s = 0.54   Random scaling with fixed symmetry: A stochastic tensor κ s (ω s ) ∈ Sym + (3) with fixed symmetry (iso-iso-scl), where the random scaling elements λ (1)  s (ω s ) = λ (2)  s (ω s ) = λ (3)  s (ω s ) are modelled by an identical log-normal random variable, is constructed.The schematic representation of the model is shown in Fig. 11, where Fig. 11(a) presents the log-normal PDF of scaling parameters {λ (i) s (ω s )} 3 i=1 , whereas the geometric visualizationin the form of equivalent spheres-of the reference tensor κ iso s and realisations {κ s i } N i=1 is displayed in Fig. 11(b).Here the radius of the sphere is determined by 1/ λ (1) As can be seen that the radius of the spheres in the realisations vary and the spherical shape of the mean is maintained in the realisations, owing to the preservation of isotropic material symmetry in the mean and realisations.
Random scaling with varying symmetry: a random tensor model κ s (ω s ) ∈ Sym + (3) similar to the previous scenario, however with a varying material symmetry (iso-ortho-scl) is showcased in Fig. 12.In Fig. 12(a), the PDF of i.i.d log-normal scaling parameters {λ (i) s (ω s )} 3 i=1 can be seen, whereas Fig. 12(b) shows the geometric interpretation of the mean κ iso s and realisations {κ s i } N i=1 .Looking from the negative y-axis, the eigenvectors of the model κ s (ω s ) conforming to the eigenvalues λ (1)   s and λ (3)     angle of π/4 in an anti-clockwise direction corresponding to x and z-axis, respectively.A change in symmetry from isotropic in the mean (spherical shape) to orthotropic (ellipsoidal shape) in the realisations is evident.Also, the fluctuation in lengths of semi-major axes-calculated by 1/ λ Random orientation with fixed symmetry: to account for directional uncertainty only, a random tensor with fixed symmetry (ortho-ortho-dir) κ r (ω r ) ∈ Sym + (3), which is the exponential of the model described in Eq.( 24), is considered.The eigenvector v(ω r ) of random tensor κ r (ω r ) related to the eigenvalue λ (3)  s is modelled as a von Mises-Fisher distribution (see Eq.( 39)) with the mean direction µ circ (the eigenvector of the reference tensor κ ortho r which corresponds to the eigenvalue λ (3)  s ), and the concentration parameter η = 75.The random rotational angle φ(ω r ) is determined using the dot product operation, φ(ω r ) = arccos(µ circ • v(ω r )), whereas the unit random rotational axis w(ω r )/ w(ω r ) is evaluated by the cross-product operation, w(ω r )/ w(ω r ) = µ circ × v(ω r )/ µ circ × v .The corresponding random rotation matrix R(ω r ) is then modelled as per the Rodrigues rotation formaula, described in Eq.( 29).Fig. 13 summarizes the model κ r (ω r ).Here Fig. 13(a) displays the mean and realisations of von Mises Fisher random variable on unit  sphere; the equivalent ellipsoidal representation of stochastic model κ r (ω r ) is portrayed in Fig. 13(b).It is evident that the shape and size of the ellipsoid in the mean and realisations remain constant; only the orientation of semi-axes of ellipsoids in the realisations fluctuate.
Uncertainty quantification : Figs.14-17 display the MC mean and standard deviation estimates of NT and THFL, whereas Fig. 18 shows the circular MC mean and circular standard deviation estimate of NHFL.Concerning the sensitivity of different stochastic models on the desired quantities, we make similar observations, and hence, draw similar conclusions on the 3D example as in the previous 2D case (from Section 4.1).That is, the findings in Figs.5-9 correspond approximately to the results in Figs.14-18.
The numerical results of both 2D and 3D proximal femur signify the prominence of incorporating different material uncertainties-scaling, orientation and material symmetry-independently into the constitutive model; the results showcase the distinct impact of stochastic models on the desired quantities of interest, such as nodal temperature and heat flux.

Conclusion
Here the task of modelling and representing random symmetric and positive definite (SPD) tensors was addressed, and some desiderata for the modelling and representation were formulated.Namely, the representation has to be SPD even under numerical approximation, and has to be such that one is able to control the material or spatial invariance resp.symmetry, both for each possible realisation and for the the expected value.Regarding the mean or expected value, it turned out that the usual expectation, which is tied to the structure of a linear space, may not be appropriate on the open convex cone of SPD tensors.The more general Fréchet mean is based on distance measurements, possibly different from the usual Euclidean distance, on the set of SPD tensors.Here some desiderata were formulated for the distance measure, which take into account the physical purpose and uncertainty involved in stochastic models of such physically symmetric and positive definite tensors.
Therefore in Appendix A the type of metric that one can consider for symmetric positive tensors and specifically for the set of SPD matrices to compute the correspondin Fréchet mean is discussed.As the Euclidean mean suffers from the so-called swelling effect, we suggest alternatives such as the geometric, log-Euclidean, and scale-rotation separated mean.
It was also proposed to follow a widely used technique to model and represent the logarithm of tensor, such that the symmetric positive definite attribute is automatically satisfied when taking the exponential in the end; the invariance properties automatically translate to the logarithm.To expose the ideas in the simplest possible setting, we consider specifically second-order SPD tensors, where the symmetry classes are quite simple.
More particularly, the focus is on modelling uncertainties in such a way as to have a fine control independently over the strength/scaling and directional attributes of the tensorgiven that the material symmetry is fixed.As we separate the strength and orientations of the random SPD matrices by spectral decomposition, the task boils down to modelling the random scaling parameters of the logarithm of a random diagonal matrix, as well as the random orthogonal matrix of eigenvectors.Such orthogonal matrices are members of the Lie group SO(d), and the well-known exponential map from the corresponding Lie algebra so(d) of skew-symmetric matrices was used to carry random ensembles in the free vector space so(d) onto such ensembles on SO(d).This results in the random orientation being modelled by symmetric circular/spherical distributions.Furthermore, a scenario of fluctuating the material symmetry around the Fréchet mean is also studied, i.e. the mean tensor belongs to a higher order of material or spatial symmetry, whereas each realisation belongs to a lower order of material symmetry.
As an application in this study, the random thermal conductivity tensor of a steadystate heat conduction on a 2D and a 3D model of a proximal femur bone was investigated, a highly anisotropic material.Three different modelling scenarios are considered: random scalings only, random scalings with fluctuating symmetry, and random orientations only.The material uncertainties are then propagated to the output responses such as nodal temperature (NT), total heat flux (THFL) and normalized heat flux (NHFL) using the Monte Carlo method (MC).Subsequently, statistics such as mean and standard deviation estimates of NT and THFL, and circular mean and circular standard deviation estimates of NHFL are determined.
From the numerical results of both 2D and 3D models, it is clear that NT is most stochastically sensitive to the model with random scalings only, and least sensitive to random orientation only.On the other hand, the randomness in the model with random scaling only has almost zero stochastic influence on THFL-determined by the Euclidean norm-but the impact of the other two random models is evident.One may also observe a similar stochastic sensitivity of the stochastic models on NHFL, as well as on THFL.Therefore, the distinctive impact of the three stochastic models on the output response signifies the importance of incorporating different material uncertainties-scaling, orientations and material symmetry-independently into the constitutive model.
In total, we have proposed a modelling scenario for tensors defining positive-definite material properties which takes into account possible material resp.spatial invariances or symmetries.Our approach allows a fine control of the scaling randomness, and independently of the orientation randomness, further separating it into the invariance properties of the mean and those of the random fluctuations.algebra Sym(d) and the Lie group Sym + (d), actually a diffeomorphism and an algebraic isomorphism between the groups (Sym + (d), ) (with the new multiplication) and the additive group (Sym(d), +) of the vector space of symmetric matrices.Thus Sym(d) may be seen as a global chart for Sym + (d).This means in particular that straight lines through the origin-geodesics with the Euclidean or Frobenius distance in Sym(d) -are mapped into geodesics in Sym + (d) with the logarithmic or log-Euclidean metric which also satisfies all desiderata listed in Section 2, and has the boundary ∂ Sym + (d) at an infinite distance from the identity.Among others, this metric is also investigated in [20] (the authors call ϑ L the chaotic geometric distance) and [55].Notice also the equality with ϑ G Eq.( 49) above for the commuting case-if -and also the similarity with the Riemannian metric ϑ R Eq.( 50) on SO(d).This means that the two metrics ϑ G Eq.( 49) and ϑ L Eq.( 52) agree on commuting subsets of Sym + (d).But such a simple Euclidean structure on the Lie algebra mixes the strength and orientational information and disregards the effects of non-commutativity.In our view the non-commutative aspect is part of what is the orientational information in the eigenvectors of a SPD matrix, and for that reason we decide to treat strength and orientation aspects separately.

A.5 Scaling-rotation distance
Thus another possibility arises from the observation that two tensors C 1 , C 2 ∈ Sym + (d) commute only if they have the same invariant subspaces.In that case it can be arranged that the eigenvector matrices according to the spectral decomposition Eq.( 13), In this instance of two commuting tensors, one has As discussed in Section 3.1, the spectral decomposition Eq.( 13) induces the representation of Eq.( 14) on a product of Lie groups, namely on: Diag + (d) × SO(d), which is itself a Lie group; this is the right half of Eq. (15).From this and Eq.( 53) comes the idea [33,22] to measure the distance on Sym + (d) by measuring separately the distance between the eigenspaces (the rotation matrices Q k ), and the distance between the eigenvalues (the scaling matrices Λ k ), and define a distance function by combining ϑ L Eq.( 52) on the commutative factor Diag + (d) with ϑ R Eq.(50) on SO(d) to a product distance for some c > 0 on the product Lie group Diag + (d) × SO(d), i.e. one measures the distance squared between the scales and the orientations separately and adds them both up.
As was pointed in Section 3.1, the representation map in Eq.( 14) is not one-to-one, so this does not carry directly over to a metric on Sym + (d) [22].Thus one sets where (Λ 1 , Q 1 ), (Λ 2 , Q 2 ) ∈ Diag + (d)×SO(d).It is not difficult to see that ϑ S also satisfies all our desiderata listed in Section 2, and the boundary ∂ Sym + (d) is infinitely far away from the identity.As we propose to rather work with H = log(C) ∈ Sym(d), where log(QΛQ T ) = Q log(Λ)Q T = QY Q T , it is important to note that the metrics mentioned before can easily be computed on H = log(C).One has

2 .
Even under numerical approximations like truncation of series resp.numerical quadrature, each realisation C(ω) has to be SPD, i.e. for all z ∈ R d \{0} : z T C(ω)z > 0. 3. ∀ω ∈ Ω : C(ω) ∈ Sym + (d) G has to be invariant under some group of transformations G ⊆ O(R d ), the invariance or symmetry class of each realisation.4. The appropriate mean C ϑ := E ϑ (C(•)) ∈ Sym + (d) has to be invariant under a possibly larger group of transformations G m , with O(d) ⊇ G m ⊇ G, the invariance or symmetry class of the mean C ϑ .

Figure 3 :
Figure 3: Stochastic tensor κ s (ω s ) with varying symmetry (iso-ortho-scl): (a) log-normal PDF of i.i.d.random scaling values {λ (i) s (ω s )} 2 i=1 , (b) geometric visualization of reference tensor κ iso s and realisations of random tensor κ s (ω s ); the geometries are enhanced by scaling the eigenvalues by a factor of 0.5

Figure 4 :
Figure 4: Stochastic tensor κ r (ω r ) with fixed symmetry (ortho-ortho-dir): (a) von Mises PDF of random rotation angle φ(ω r ), (b) realisations of circular random variable φ(ω r ) on unit circle with resultant mean vector (solid straight line), (c) geometric visualization of reference tensor κ ortho r and realisations of random tensor κ r (ω r ); the geometries are enhanced by scaling the eigenvalues by a factor of 0.5, and the orientations with respect to the mean are scaled by a factor of 1.5

Figure 5 :
Figure 5: MC mean estimate of nodal temperature

Figure 6 :
Figure 6: MC standard deviation estimate of nodal temperature

Figure 8 :
Figure 8: MC standard deviation estimate of total heat flux (Euclidean norm)

Figure 9 :Figure 10 :
Figure 9: MC circular mean (vector representation) and circular standard deviation estimates (dimensionless quantity) of normalized heat flux

Figure 12 :
Figure 12: Stochastic tensor κ s (ω s ) with varying symmetry (iso-ortho-scl): (a) log-normal PDF of i.i.d random scaling values {λ (i) s (ω s )} 3 i=1 , (b) geometric visualization of reference tensor κ iso s and realisations of random tensor κ s (ω s ); the geometries are enhanced by scaling the eigenvalues by a factor of 0.5

s
(ω s )-of the ellipsoids represents the aspect of varying scaling parameters.

Figure 13 :
Figure 13: Stochastic tensor κ r (ω r ) with fixed symmetry (ortho-ortho-dir): (a) realisations of von Mises Fisher random variable v(ω r ) on unit sphere with mean µ circ , (b) geometric visualization of reference tensor κ ortho r and realisations of random tensor κ r (ω r ); the geometries are enhanced by scaling the eigenvalues by a factor of 0.5, and the orientations with respect to the mean are scaled by a factor of 1.5

Table 3 :
Reference conductivity tensors with corresponding eigenvalues of 3D femur