Piotr Graczyk
email: piotr.graczyk@univ-angers.fr

Hideyuki Ishi
email: hideyuki@sci.osaka-cu.ac.jp

Bartosz Kołodziejek
email: b.kolodziejek@mini.pw.edu.pl

† And

Hélène Massam
email: massamh@mathstat.yorku.ca

MODEL SELECTION IN THE SPACE OF GAUSSIAN MODELS INVARIANT BY SYMMETRY

Keywords: colored graph, conjugate prior, covariance selection, invariance, permutation symmetry 1

We consider multivariate centred Gaussian models for the random variable Z = (Z1, . . . , Zp), invariant under the action of a subgroup of the group of permutations on {1, . . . , p}. Using the representation theory of the symmetric group on the field of reals, we derive the distribution of the maximum likelihood estimate of the covariance parameter Σ and also the analytic expression of the normalizing constant of the Diaconis-Ylvisaker conjugate prior for the precision parameter K = Σ -1 . We can thus perform Bayesian model selection in the class of complete Gaussian models invariant by the action of a subgroup of the symmetric group, which we could also call complete RCOP models. We illustrate our results with a toy example of dimension 4 and several examples for selection within cyclic groups, including a high dimensional example with p = 100.

1. Introduction. Let V = {1, . . . , p} be a finite index set and let Z = (Z 1 , . . . , Z p ) be a multivariate random variable following a centred Gaussian model N p (0, Σ). Let S p denote the symmetric group on V , that is, the group of all permutations on {1, . . . , p} and let Γ be a subgroup of S p . A centred Gaussian model is said to be invariant under the action of Γ if for g ∈ Γ, g • Σ • g = Σ (here we identify a permutation g with its permutation matrix).

Given n data points Z (1) , . . . , Z (n) from a Gaussian distribution, our aim in this paper is to do Bayesian model selection within the class of models invariant by symmetry, that is, invariant under the action of some subgroup Γ of S p on V . Given the data, our aim is therefore to identify the subgroup Γ ⊂ S p such that the model invariant under Γ has the highest posterior probability.

Gaussian models invariant by symmetry have been considered by [START_REF] Andersson | Invariant normal models[END_REF] and [START_REF] Andersson | Distribution of eigenvalues in multivariate statistical analysis[END_REF]. Gaussian models invariant under the action of a subgroup Γ ⊂ S p have also been considered in [START_REF] Andersson | Symmetry and lattice conditional independence in a multivariate normal distribution[END_REF], [START_REF] Madsen | Invariant normal models with recursive graphical Markov structure[END_REF] and [START_REF] Højsgaard | Graphical Gaussian models with edge and vertex symmetries[END_REF], but in this paper, we will consider complete Gaussian models without predefined conditional independencies. In [START_REF] Andersson | Symmetry and lattice conditional independence in a multivariate normal distribution[END_REF], the reader can also find references to earlier works dealing with particular symmetry models such as, for example, the circular symmetry model of [START_REF] Olkin | Testing and estimation for a circular and stationary model[END_REF] that we will consider further (Section 5). These works were concentrating on the derivation of the maximum likelihood estimate of Σ and on testing the hypothesis that models were of a particular type. Our work is a first step towards Bayesian model selection in the class of models invariant by symmetry.

Just like the classical papers mentioned above, the fundamental algebraic tool we use in this work is the irreducible decomposition theorem for the matrix representation of the group Γ, which in turn means that, through an adequate change of basis, any matrix X in P Γ , the cone of positive definite matrices invariant under the subgroup Γ of S p , can be written as

(1) X = U Γ •      M K1 (x 1 ) ⊗ I k1/d1
M K2 (x 2 ) ⊗ I k2/d2 . . .

M KL (x L ) ⊗ I kL/dL      • U Γ ,
where U Γ is an orthogonal matrix, k i , d i , r i , i = 1, . . . , L, are integer constants called structure constants we will define later, such that k i /d i are also integers and M Ki (x i ), i = 1, . . . , L is a real matrix representation of an r i × r i Hermitian matrix x i with entries in K i = R, C or H with d i = dim R K i = 1, 2, 4. Moreover the empty entries in the middle matrix of (1) are all equal to 0 so that this middle matrix is block diagonal. Finally, A ⊗ B denotes the Kronecker product of matrices A and B.

Even though some of the results in the theoretical part of the present paper were already developed in [START_REF] Andersson | Invariant normal models[END_REF], we decided to present full proofs for two reasons. First, we think that our arguments are more concrete and should be easier to understand for the reader who is not familiar with representation theory. Second, our results are more explicit and, in particular, we are able to completely solve the case when Γ is a cyclic group (i.e. Γ has one generator), that is, we show how U Γ and all structure constants (k i , d i , r i ) L i=1 can be computed explicitly.

Let us consider the following example.

EXAMPLE 1. For p = 3 and Γ = S 3 , the cone of positive definite matrices X invariant under Γ, that is, such that X ij = X σ(i)σ(j) for all σ ∈ Γ, is

P Γ =     
Γ := v 1 v 2 v 3 ∈ O(3) with v 1 :=   1/ √ 3 1/ √ 3 1/ √ 3   , v 2 :=   2/3 -1/ √ 6 -1/ √ 6   , v 3 :=   0 1/ √ 2 -1 √ 2   , and   a b b b a b a a b   = U Γ •   a + 2b a -b a -b   • U Γ , Here L = 2, k 1 /d 1 = 1, k 2 /d 2 = 2, K 1 = K 2 = R, d 1 = d 2 = 1.
We see immediately in the example above that, following the decomposition (1), the trace Tr [X] = a + 2b + 2(a -b) and the determinant Det (X) = (a + 2b)(a -b) 2 can be readily obtained. Similarly, using (1) allows us to easily obtain Det (X) and Tr [X] in general.

In Section 3, we will see that having the explicit formulas for Det (X) and Tr [X], in turn, allows us to derive the analytic expression of the Gamma function on P Γ , defined as Tr[X] ϕ Γ (X) dX, where ϕ Γ (X) dX is the invariant measure on P Γ (see Definition 13 and Proposition 9) and dX denotes the Euclidean measure on the space Z Γ with the trace inner product.

Γ PΓ (λ) := PΓ Det (X) λ e -
With the Gamma integral on P Γ , we can derive the analytic expression of the normalizing constant I Γ (δ, D) of the Diaconis-Ylvisaker conjugate prior on K = Σ -1 with density, with respect to the Euclidean measure on P Γ , equal to

f (K; δ, D) = 1 I Γ (δ, D) Det (K) (δ-2)/2 e -1 2 Tr[K•D] 1 PΓ (K)
for appropriate values of the scalar hyper-parameter δ and the matrix hyper-parameter D ∈ P Γ . By analogy with the G-Wishart distribution, defined in the context of the graphical Gaussian models, Markov with respect to an undirected graph G on the cone P G of positive definite matrices with zero entry (i, j) whenever there is no edge between the vertices i and j in G, (see [START_REF] Maathuis | Handbook of Graphical Models[END_REF]), we can call the distribution with density f (K; δ, D), the RCOP-Wishart (RCOP is the name coined in [START_REF] Højsgaard | Graphical Gaussian models with edge and vertex symmetries[END_REF] for graphical Gaussian models with restrictions generated by permutation symmetry). It is important to note here that if Σ is in P Γ , so is K = Σ -1 so that K can also be decomposed according to (1). Equipped with all these results, we compute the Bayes factors comparing models pairwise and perform model selection. We will indicate in Section 4 how to travel through the space of cyclic subgroups of the symmetric group.

In Section 3, we also derive the distribution of the maximum likelihood estimate (henceforth abbreviated MLE) of Σ and show that for n ≥ max i=1,...,L ridi ki it has a density equal to

Det (X) n/2 e -1 2 Tr[X•Σ -1 ] Det (2Σ) n/2 Γ PΓ ( n 2 ) ϕ Γ (X)1 PΓ (X).
Clearly, the key to computing the Gamma integral on P Γ , the normalizing constant I Γ (δ, D) or the density of the MLE of Σ is, for each Γ ⊂ S p , to obtain the block diagonal matrix with diagonal block entries M Ki (x i ) ⊗ I ki/di , i = 1, . . . , L, in the decomposition (1). In principle, we have to derive the invariant measure ϕ Γ and find the structure constants k i , d i , r i , i = 1, . . . , L. This goal can be achieved by constructing an orthogonal matrix U Γ and using (1). However, doing so for every Γ visited during the model selection process is computationally heavy.

We will show that for small to moderate dimensions, we can obtain the constants k i , d i , r i , i = 1, . . . , L as well as the expression of Det (X) and ϕ Γ (X) without having to compute U Γ . Indeed, as indicated in Remark 7, for any X ∈ P Γ , Det (X) admits a unique irreducible factorization of the form (2)

Det (X) = L i=1 Det (M Ki (x i )) ki/di = L j=1 f j (X) aj (X ∈ Z Γ ),
where each a j is a positive integer, each f j (X) is an irreducible polynomial of X ∈ Z Γ , and

f i = f j if i = j.
The constants k i , d i , r i are obtained by identification of the two expressions of Det (X) in (2). Factorization of a homogeneous polynomial Det (X) can be performed using standard software such as either MATHEMATICA or PYTHON.

Due to computational complexity, for bigger dimensions, it is difficult to obtain the irreducible factorization of Det (X). For special cases such as the case where the subgroup Γ is a cyclic group, we give (Section 2) a simple construction of the matrix U Γ and thus, for any dimension p, we can do model selection in the space of models invariant under the action of a cyclic group. We argue that restriction to cyclic groups is not as limiting as it may look. The formula for the number of different colorings c p = # { P Γ ; Γ ⊂ S p } for given p is unknown. Obviously, it is bounded from above by the number of all subgroups of S p , because different subgroups may produce the same coloring (e.g. in Example 1 we have P S3 = P (1,2,3) ). On the other hand, it is known (see Lemma 17) that c p is bounded from below by the number of distinct cyclic subgroups, which grows rapidly with p (see OEIS1 sequence A051625). In particular, for p = 182 , we have c p ∈ (7.1 • 10 14 , 7.6 • 10 18 ), see also Table 1. The lower bound for c p indicates that the colorings obtained from cyclic subgroups form a rich subfamily of all possible colorings. Let us consider the more general situation of Gaussian graphical models with conditional independence structure encoded by a non complete graph G. Then one can introduce symmetry restrictions (RCOP) by requiring that the precision matrix K is invariant under some subgroup Γ of S p . However, when G is not complete, not all subgroups are suited to the problem. In such cases, one has to require that Γ belongs to the automorphism group Aut(G) of G. If a graph G is sparse, then Aut(G) may be very small and it is natural to expect that the vast majority of subgroups of Aut(G) are actually cyclic. Moreover, finding the structure constants for a general group is much more expensive and in some situations it may not be worth to consider the problem in its full generality. We consider our work as a first step towards the rigorous analytical treatment of Bayesian model selection in the space of graphical Gaussian models invariant under the action of Γ ⊂ S p when conditional independencies are allowed. Finally, we expect the statistical interpretability of cyclic models to be easier than that of general groups.

The procedure to do model selection will be described in Section 4 and we will illustrate this procedure with Frets' data (see [START_REF] Frets | Heredity of head form in man[END_REF]) and several examples for selection within cyclic groups, including a high dimensional example with p = 100 (Section 5).

Most technical parts of the paper are postponed to the Appendix.

2.

Preliminaries. Representation theory has long been known to be very useful in statistics, cf. [START_REF] Diaconis | Group representations in probability and statistics[END_REF]. However, the representation theory over R that we need in this paper, is less known to the statisticians than the standard one over C. We will thus now recall the basic notions that we need.

2.1. Notation. Let Mat(n, m; R), Sym(n; R) denote the linear spaces of real n × m matrices and symmetric real n × n matrices, respectively. Let Sym + (n; R) be the cone of symmetric positive definite real n × n matrices. A denotes the transpose of a matrix A. Det and Tr denote the usual determinant and trace in Mat(n, n; R). Let p denote the fixed number of vertices of a graph and let S p denote the symmetric group. We write permutations in cycle notation, meaning that (i 1 , i 2 , . . . , i n ) maps i j to i j+1 for j = 1, . . . , r -1 and i n to i 1 . By σ 1 , . . . , σ k we denote the group generated by permutations σ 1 , . . . , σ k . The composition (product) of permutations σ, σ ∈ S p will be denoted by σ • σ . DEFINITION 2. For a subgroup Γ ⊂ S p , we define the space of symmetric matrices invariant under Γ, or the vector space of colored matrices,

Z Γ := x ∈ Sym(p; R) ; x ij = x σ(i)σ(j) for all σ ∈ Γ ,
and the cone of positive definite matrices valued in Z Γ ,

P Γ := Z Γ ∩ Sym + (p; R).
We note that the same colored space and cone can be generated by two different subgroups: in Example 1, the subgroup Γ = (1, 2, 3) generated by the permutation σ

= (1, 2, 3) is such that Γ = Γ but Z Γ = Z Γ . Let us define Γ * = σ * ∈ S p ; x ij = x σ * (i)σ * (j) for all x ∈ Z Γ .
Clearly, Γ is a subgroup of Γ * and Γ * is the unique largest subgroup of S p such that

Z Γ * = Z Γ or, equivalently, such that the Γ * -and Γ-orbits in { {v 1 , v 2 } ; v i ∈ V, i = 1, 2 } are the same. The group Γ * is called the 2 * -closure of Γ. The group Γ is said to be 2 * -closed if Γ = Γ * .
Subgroups which are 2 * -closed are in bijection with the set of colored spaces. These concepts have been investigated in [START_REF] Wielandt | Permutation groups through invariant relations and invariant functions[END_REF]; [START_REF] Serre | On partitions and permutation groups on unordered sets[END_REF] along with a generalization to regular colorings in [START_REF] Siemons | Automorphism groups of graphs[END_REF]. The combinatorics of 2 * -closed subgroups is very complicated and little is known in general, (Graham, Grötschel andLovász, 1995, p. 1502). In particular, the number of such subgroups is not known, but brute-force search for small p indicates that this number is much less than the number of all subgroups of S p (see Table 1). Even though cyclic subgroups of S p are in general not 2 * -closed, each cyclic group corresponds to a different coloring (see Lemma 17).

For a permutation σ ∈ S p , denote its matrix by

(3) R(σ) := p i=1 E σ(i)i ,
where E ab is the p × p matrix with 1 in the (a, b)-entry and 0 in other entries. The condition

x σ(i)σ(j) = x ij is then equivalent to R(σ) • x • R(σ) = x. Consequently, Z Γ = x ∈ Sym(p; R) ; R(σ) • x • R(σ) =
x for all σ ∈ Γ . (4) DEFINITION 3. Let π Γ : Sym(p; R) → Z Γ be the projection such that for any x ∈ Sym(p; R) the element π Γ (x) ∈ Z Γ is uniquely determined by

Tr [x • y] = Tr [π Γ (x) • y] (y ∈ Z Γ ). (5)
In view of (4), it is clear that

π Γ (x) = 1 |Γ| σ∈Γ R(σ) • x • R(σ) (6)
satisfies the above definition. Here |Γ| denotes the order of Γ.

2.2. Z Γ as a Jordan algebra. To prove (1), that is, Theorem 4 below, we need to view P Γ as the cone of squares of a Jordan algebra. We recall here the fundamentals of Jordan algebras, cf. [START_REF] Faraut | Analysis on symmetric cones[END_REF]. A Euclidean Jordan algebra is a Euclidean space A (endowed with the scalar product denoted by •, • ) equipped with a bilinear mapping (product)

A × A (x, y) → x • y ∈ A
such that for all x, y, z in A:

(i) x • y = y • x, (ii) x • ((x • x) • y) = (x • x) • (x • y), (iii) x, y • z = x • y, z .
A Euclidean Jordan algebra is said to be simple if it is not a Cartesian product of two Euclidean Jordan algebras of positive dimensions. We have the following result.

PROPOSITION 1. The Euclidean space Z Γ with inner product x, y = Tr [x • y] and the Jordan product

x • y = 1 2 (x • y + y • x), (7)
is a Euclidean Jordan algebra. This algebra is generally non-simple.

PROOF. Since Z Γ is a subset of the Euclidean Jordan algebra Sym(p; R), if it is endowed with Jordan product (7), conditions (i)-(iii) are automatically satisfied. Moreover, characterization (4) of Z Γ implies that the Jordan product is closed in

Z Γ , that is, R(σ) • (x • y) = (x • y) • R(σ) for all x, y ∈ Z Γ and σ ∈ Γ. The result follows.
Up to linear isomorphism, there are only five kinds of Euclidean simple Jordan algebras. Let K denote the set of either the real numbers R, the complex ones C or the quaternions H. Let us write Herm(r; K) for the space of r × r Hermitian matrices valued in K. Then Sym(r; R), r ≥ 1, Herm(r; C), r ≥ 2, Herm(r; H), r ≥ 2 are the first three kinds of Euclidean simple Jordan algebras and they are the only ones that will concern us. The determinant and trace in Jordan algebras Herm(r; K) will be denoted by det and tr respectively, so that they can be easily distinguished from the determinant and trace in Mat(n, n; R) which we denote by Det and Tr.

To each Euclidean Jordan algebra A, one can attach the set Ω of Jordan squares, that is, Ω = { x • x ; x ∈ A }. The interior Ω of Ω is a symmetric cone, that is, it is self-dual and homogeneous. We say that Ω is irreducible if it is not the Cartesian product of two convex cones. One can prove that an open convex cone is symmetric and irreducible if and only if it is the symmetric cone Ω of some Euclidean simple Jordan algebra. Each simple Jordan algebra corresponds to a symmetric cone. The first three kinds of irreducible symmetric cones are thus, the symmetric positive definite real matrices Sym + (r; R) for r ≥ 1, complex Hermitian positive definite matrices Herm + (r; C), and quaternionic Hermitian positive definite matrices Herm + (r; H), r ≥ 2.

It follows from Definition 2 and Proposition 1 that P Γ is a symmetric cone. In (Faraut and Korányi, 1994, Proposition III.4.5) it is stated that any symmetric cone is a direct sum of irreducible symmetric cones. As it will turn out, only three out of the five kinds of irreducible symmetric cones may appear in this decomposition.

Moreover, we will want to represent the elements of the symmetric cones in their real symmetric matrix representations. So, we recall that both Herm + (r; C) and Herm + (r; H) can be realized as real symmetric matrices, but of bigger dimension. For z = a

+ b i ∈ C define M C (z) = a -b b a
. The function M C is a matrix representation of C. Similarly, any r × r complex matrix can be realized as a (2r) × (2r) real matrix by setting the correspondence Mat(r, r; C) z i,j 1≤i,j≤r M C (z i,j ) 1≤i,j≤r ∈ Mat(2r, 2r; R), that is, an (i, j)-entry of a complex matrix is replaced by its 2 × 2 real matrix representation. Note that M C maps the space Herm(r; C) of Hermitian matrices into the space Sym(2r; R) of symmetric matrices. For example,

M C a c -d i c + d i b =     a 0 c d 0 a -d c c -d b 0 d c 0 b     .
Moreover, by direct calculation one sees that

Det     a 0 c d 0 a -d c c -d b 0 d c 0 b     = det a c -d i c + d i b 2 .
It can be shown that, in general,

Det (M C (Z)) = [det (Z)] 2 and Tr [M C (Z)] = 2 tr[Z] (Z ∈ Herm(r; C)). (8)
Similarly, quaternions can be realized as a 4 × 4 matrix:

a + b i + c j + d k a + b i -c + d i c + d i a -b i     a -b -c -d b a d -c c -d a b d c -b a     .
Then, quaternionic r × r matrices are realized as (4r) × (4r) real matrices. Thus, M H maps Herm(r; H) into Sym(4r; R). Moreover, it is true that

Det (M H (Z)) = [det (Z)] 4 and Tr [M H (Z)] = 4 tr[Z] (Z ∈ Herm(r; H)). (9)
2.3. Basics of representation theory over real fields. We will show, in this section, that the correspondence σ → R(σ) defined in (3) is a representation of Γ and that, as for all representations of a finite group, through an appropriate change of basis, matrices R(σ), σ ∈ Γ can be simultaneously written as block diagonal matrices with the number and dimensions of these block matrices being the same for all σ ∈ Γ. This, in turn, will allow us to write any matrix in Z Γ under the form (1). To do so, we first need to recall some basic notions and results of the representation theory of groups over the reals. For further details, the reader is referred to [START_REF] Serre | On partitions and permutation groups on unordered sets[END_REF].

For a real vector space V , we denote by GL(V ) the group of linear automorphisms on V . Let G be a finite group.

DEFINITION 4. A function ρ : G → GL(V ) is called a representation of G over R if it is a homomorphism, that is ρ(g g ) = ρ(g) ρ(g ) (g, g ∈ G).
The vector space V is called the representation space of ρ.

If dim V = n, taking a basis {v 1 , . . . , v n } of V , we can identify GL(V ) with the group GL(n; R) of all n × n non-singular real matrices. Then a representation ρ : G → GL(V ) corresponds to a group homomorphism B : G → GL(n; R) for which (10) ρ(g)v j = n i=1 B ij (g)v i .
We call B the matrix expression of ρ with respect to the basis {v 1 , . . . , v n }.

By definition, we have R : Γ → GL(p; R) and R(σ

• σ ) = R(σ) • R(σ ) for all σ, σ ∈ S p . Thus, R is a representation of Γ over R. Regarding ρ(σ) = R(σ) ∈ GL(R p
) as an operator on V = R p via the standard basis v i = e i ∈ R p , i = 1, . . . , p, we see that (10) holds trivially with B = R.

Example 6 below gives an illustration of the representation R(σ) and also an illustration of all the notions and results we are about to state now, leading to the expression (1) of a matrix in

Z Γ . DEFINITION 5. A linear subspace W ⊂ V is said to be G-invariant if ρ(g)w ∈ W (w ∈ W, g ∈ G).
A representation ρ is said to be irreducible if the only G-invariant subspaces are non-proper, that is, whole V and {0}. A restriction of ρ to a G-invariant subspace W is a subrepresentation. Two representations, ρ : G → GL(V ) and ρ : G → GL(V ) are equivalent if there exists an isomorphism of vector spaces :

V → V with (ρ(g)v) = ρ (g) (v) (v ∈ V, g ∈ G).
We note that, as we have discussed for the case B = R, a group homomorphism B : G → GL(n; R) defines a representation of G on R n naturally. We see that B is a matrix expression of a representation (ρ, V ) if and only if B and ρ are equivalent via the map :

R n (x i ) n i=1 → n i=1 x i v i ∈ V , that is, (B(g)x) = ρ(g) (x) for x ∈ R n .
Here {v 1 , . . . , v n } denotes a fixed basis of V . Therefore, two representations (ρ, V ) and (ρ , V ) are equivalent if and only if they have the same matrix expressions with respect to appropriately chosen bases. We shall write ρ ∼ B if ρ has a matrix expression B with respect to some basis.

Let (ρ, V ) be a representation of G, and B : G → GL(n; R) be a matrix expression of ρ with respect to a basis {v 1 , . . . , v n } of V . Then it is known that the function χ ρ : G g → Tr B(g) = n i=1 B ii (g) ∈ R is independent of the choice of the basis {v 1 , . . . , v n }. The function χ ρ is called a character of the representation ρ. The function χ ρ characterizes the representation ρ in the following sense.

LEMMA 2. Two representations (ρ, V ) and (ρ , V ) of a group G are equivalent if and only if χ ρ = χ ρ .

We apply this lemma in practice to know whether two given representations are equivalent or not.

It is known that, for a finite group G, the set Λ(G) of equivalence classes of irreducible representations of G is a finite set. We fix the group homomorphisms

B α : G → GL(k α ; R), α ∈ A, indexed by a finite set A so that Λ(G) = { [B α ] ; α ∈ A }, where [B α ] denotes the equivalence class of B α .
Let (ρ, V ) be a representation of G. Then there exists a G-invariant inner product on V . In fact, from any inner product •, • 0 on V , one can define such an invariant inner product

•, • by v, v := g∈G ρ(g)v, ρ(g)v 0 for v, v ∈ V .
In what follows, we fix a G-invariant inner product on V .

If W is a G-invariant subspace, the orthogonal complement W ⊥ is also a G-invariant subspace. Thus, any representation ρ can be decomposed into a finite number of irreducible subrepresentations

ρ = ρ 1 ⊕ . . . ⊕ ρ K (11) along the orthogonal decomposition V = V 1 ⊕ • • • ⊕ V K , where ρ i is the restriction of ρ to the G-invariant subspace V i , i = 1, . . . , K. Let r α be the number of subrepresentations ρ i such that ρ i ∼ B α . Although the irreducible decomposition (11) of V is not unique in general, r α is uniquely determined. We have ρ ∼ rα>0 B ⊕rα α , (12)
where rα>0 r α = K. To see this, let V (B α ) be the direct sum of subspaces V i for which

ρ i ∼ B α . The space V (B α ) is called the B α -component of V . If r α > 0, gathering an appropriate basis of each V i , the matrix expression of the subrepresentation of ρ on V (B α ) becomes (recall that B α (g) ∈ GL(k α ; R)) B α (g) ⊕rα =      B α (g) B α (g) . . . B α (g)      = I rα ⊗ B α (g) ∈ GL(r α k α ; R) (g ∈ G).
Moreover, V is decomposed as V = rα>0 V (B α ). Therefore, taking a basis of V by gathering the bases of V (B α ), we obtain (12).

In particular, for G = Γ ⊂ S p and (ρ,

V ) = (R, R p ), if we let { α ∈ A ; r α > 0 } =: {α 1 , α 2 , .
. . , α L } and if we denote by U Γ an orthogonal matrix whose column vectors form orthonormal bases of V (B α1 ), . . . , V (B αL ) successively, then for σ ∈ Γ, we have ( 13)

U Γ • R(σ) • U Γ =      I r1 ⊗ B α1 (σ) I r2 ⊗ B α2 (σ) . . . I rL ⊗ B αL (σ)      .
Note that, since the left hand side of ( 13) is an orthogonal matrix, matrices B αi (σ), i = 1, . . . , L, are orthogonal. In the general case, B α (g) are orthogonal if we work with a Ginvariant inner product. In what follows, we will consider only the case G = Γ and ρ = R.

Note that the usual inner product on V = R p is clearly Γ-invariant.

The actual formula for B αi (σ) obviously depends on the choice of U Γ and hence, on the choice of orthonormal basis of R p . To ensure simplicity of formulation of our next result (Lemma 3), we will work with special orthonormal bases of V (B α1 ), . . . , V (B αL ), which together constitute a basis of R p . Such bases always exist and will be defined in the next section. Usage of these bases is not indispensable for the proof of our main result (Theorem 4), but simplifies it greatly. EXAMPLE 6. Let p = 4 and let Γ = {id, (1, 2)(3, 4)} be the subgroup of S 4 generated by σ = (1, 2)(3, 4). The matrix representation of σ in the standard basis

(e i ) i of R 4 is R(σ) =     0 1 0 0 1 0 0 0 0 0 0 1 0 0 1 0     ,
which has the two eigenvalues 1 and -1 with multiplicity 2 for each. We choose the following orthonormal eigenvectors of R(σ):

u 1 = 1 √ 2 (e 1 + e 2 ), u 2 = 1 √ 2 (e 3 + e 4 ), u 3 = 1 √ 2 (e 1 -e 2 ), u 4 = 1 √ 2 (e 3 -e 4 )
and let U Γ = (u 1 , u 2 , u 3 , u 4 ). The corresponding eigenspaces V i = Ru i are invariant under R(σ) and R(id) = I 4 . As V i , i = 1, . . . , 4, are 1-dimensional, the subrepresentations defined by

ρ i (γ) = R(γ)| Vi (γ ∈ Γ)
are irreducible. We have the decomposition (11) of R:

R = ρ 1 ⊕ ρ 2 ⊕ ρ 3 ⊕ ρ 4 .
The matrix expressions of ρ 1 and ρ 2 are equal to

B 1 (γ) = (1) for all γ ∈ Γ, since ρ i (γ)v = v for v ∈ V i , i = 1, 2. We have r 1 = 2.
The matrix expressions of ρ 3 and ρ 4 are both equal to B 2 (γ) = sign(γ) for all γ ∈ Γ, since

ρ i (id)v = v and ρ i (σ)v = -v for v ∈ V i for i = 3, 4. We have r 2 = 2.
The representations ρ 1 and ρ 3 are not equivalent, which can be seen by looking at the characters: χ ρ1 = 1, χ ρ3 (γ) = sign(γ), which are not equal.

In the basis u 1 , u 2 , u 3 , u 4 , the matrix of R(γ) is (compare with ( 13))

    1 0 0 0 0 1 0 0 0 0 sign(γ) 0 0 0 0 sign(γ)     = B 1 (γ) ⊕2 ⊕ B 2 (γ) ⊕2 = U Γ • R(γ) • U Γ .
This is the decomposition (12) of R in the basis (u 1 , u 2 , u 3 , u 4 ).

2.4. Block diagonal decomposition of Z Γ . So far, we have shown that through an appropriate change of basis, the representation (R, R p ) of Γ can be expressed as the direct sum (12) of irreducible subrepresentations. We now want to turn our attention to the elements of Z Γ .

A linear operator T : V → V is said to be an intertwining operator of the representation (ρ, V ) if T • ρ(g) = ρ(g) • T holds for all g ∈ G. In our context, since (4) can be rewritten as

Z Γ = { x ∈ Sym(p; R) ; x • R(σ) = R(σ) • x for all σ ∈ Γ } , Z Γ is the set of symmetric intertwining operators of the representation (R, R p ).
Let End Γ (R p ) denote the set of all intertwining operators of the representation (R, R p ) of Γ. Recall that the set A enumerates the elements of Λ(Γ), the finite set of all equivalence classes of irreducible representations of Γ. From ( 12) and ( 13), it is clear that to study End Γ (R p ), it is sufficient to study the sets,

End Γ (V α ) = { T ∈ Mat(k α , k α ; R) ; T • B α (σ) = B α (σ) • T for all σ ∈ Γ } , α ∈ A,
of all intertwining operators of the irreducible representation B α , where V α := R kα is the representation space of B α equipped with a Γ-invariant inner product. Indeed, we have

V (B α ) = I rα ⊗ V α .
To reach our main result in this section, Theorem 4, we use the result from [START_REF] Serre | On partitions and permutation groups on unordered sets[END_REF], Page 108) that, since the representation B α is irreducible, the space End Γ (V α ) is isomorphic either to R, C, or the quaternion algebra H. Let

f α : K α → End Γ (V α ), denote this isomorphism, where K α is R, C, or H. Let d α := dim R End Γ (V α ) = dim R K α ∈ {1, 2, 4}.
The representation space V α becomes a vector space over

K α of dimension k α /d α via q • v := f α (q)v (q ∈ K α , v ∈ V α ).
Clearly the space RI kα of scalar matrices is contained in End

Γ (V α ). If d α = 1 = dim R RI kα , we have End Γ (V α ) = RI kα . Further, if d α = 2, take a C-basis {v 1 , . . . , v kα/2 } of V α in such a way that {v 1 , . . . , v kα/2 , i • v 1 , . . . , i • v kα/2 } is an orthonormal R-basis of V α .
We identify R kα and V α via this R-basis. Then, the action of q = a + bi ∈ C on w ∈ R kα V α is expressed as

q • w = aI kα/2 -bI kα/2 bI kα/2 aI kα/2 w = M C (a + bi) ⊗ I kα/2 w. Thus, if d α = 2, then End Γ (V α ) = M C (q) ⊗ I kα/2 ; q ∈ C = M C (C) ⊗ I kα/2 .
Similarly, when

K α = H, take an H-basis {v 1 , . . . , v kα/4 } of V α so that {v 1 , . . . , v kα/4 , i • v 1 , . . . , i • v kα/4 , j • v 1 , . . . , j • v kα/4 , k • v 1 , . . . , k • v kα/4 } is an orthonormal R-basis of V α . The action of Q ∈ H on V α is expressed as M H (Q) ⊗ I kα/4
with respect to this basis.

In this way we have proved the following result.

LEMMA 3. For each α ∈ A, one has

End Γ (V α ) = M Kα (K α ) ⊗ I kα/dα . ( 14
)
Recall the orthogonal matrix U Γ in (13). It will be shown in the Appendix that, using Lemma 3, we can represent Z Γ as follows.

THEOREM 4. The Jordan algebra Z Γ is isomorphic to L i=1 Herm(r i ; K i ) through the map ι :

L i=1 Herm(r i ; K i ) (x i ) L i=1 → X ∈ Z Γ given by (15) X := U Γ      M K1 (x 1 ) ⊗ I k1/d1 M K2 (x 2 ) ⊗ I k2/d2 . . . M KL (x L ) ⊗ I kL/dL      U Γ .
2.5. Determining the structure constants and invariant measure on P Γ . As mentioned in the introduction, in order to derive the analytic expression of the Gamma-like functions on P Γ , we need the structure constants r i , d i , k i as well as the invariant measure ϕ Γ . However, due to Proposition 9 below, ϕ Γ (X) is expressed in terms of the polynomials det(x i ), where x i ∈ Herm(r i ; K i ), i = 1, . . . , L, coming from decomposition (15). These can be derived from the decomposition of Z Γ as indicated in Sections 2.1-2.4 above. Let us note that the constants d i and k i depend only on the group Γ, while r i depend on a particular representation of Γ, which is R.

In view of decomposition ( 15), for

X ∈ Z Γ , define φ i (X) = x i ∈ Herm(r i ; K i ) for i = 1, . . . , L. COROLLARY 5. For X ∈ Z Γ , one has (16) Det (X) = L i=1 det(φ i (X)) ki .
PROOF. By (15), we have

Det (X) = L i=1 Det M Ki (x i ) ⊗ I ki/di = L i=1 Det (M Ki (x i )) ki/di = L i=1 det(x i ) di ki/di = L i=1 det(x i ) ki ,
whence follows the formula. We have used ( 8) and ( 9) for the third equality above.

REMARK 7. Let us note that ( 16) gives us a simple way to find the structure constants {k i , r i , d i } L i=1 . Indeed, assume that we have an irreducible factorization ( 17)

Det (X) = L j=1 f j (X) aj (X ∈ Z Γ ),
where each a j is a positive integer, each f j (X) is an irreducible polynomial of X ∈ Z Γ , and

f i = f j if i = j.
Since the determinant polynomial of a simple Jordan algebra is always irreducible (Upmeier, 1986, Lemma 2.3 (1)), comparing ( 16) and ( 17), we obtain that L = L , and that, for each j, there exists i such that f j (X) aj = det(φ i (X)) ki . Then k i = a j and r i is the degree of f j (X) = det(φ i (X)). From the structure theory of Jordan algebras, we see that, if E j ∈ Z Γ is the gradient of f j (X) at X = I p , then the linear operator

P j : Z Γ x → E j •x ∈ Z Γ coincides with the projection L i=1 Herm(r i ; K i ) → Herm(r i ; K i ). In particular, dim R Herm(r i ; K i ) = r i + d i r i (r i -1)/2 is
equal to the rank of P j , from which we can determine d i when r i > 1.

If r i = 1, the determination of d i is not needed for writing the block decomposition of Z Γ , since in this case R = Herm(1; R) = Herm(1; C) = Herm(1; H) and, if k i is divisible by 2 or by 4, we have

M Ki (x i ) ⊗ I ki/di = x i I ki .
The practical significance of the method proposed in this Remark is that neither representation theory nor group theory is used. It is a strong advantage when we consider colorings corresponding to a large number of different groups, for which finding structure constants is very complicated. REMARK 8. The factorization of multivariate polynomials over an algebraic number field can be done for example in PYTHON (see sympy.polys.polytools.factor). Indeed, in our setting, the irreducible factorization over the real number field coincides with the one over the real cyclotomic field

Q ζ + 1 ζ =    ϕE(M )/2-1 k=0 q k ζ + 1 ζ k ; q k ∈ Q, k = 0, 1, . . . , ϕ E (M )/2 -1    ,
where ζ is the primitive M -th root e 2πi/M of unity with M being the least common multiple of the orders of elements σ ∈ Γ, and ϕ E (M ) is the number of positive integers up to M that are relatively prime to M (Serre, 1977, Section 12.3).

An example showing the utility of Remark 7 can be found in the Appendix (see Example 18). 2.6. Construction of the orthogonal matrix U Γ when Γ is cyclic. We now show that, when the group Γ is generated by one permutation σ ∈ S p , the orthogonal matrix U Γ can be constructed explicitly, and we obtain the structure constants r i , k i and d i easily.

Let us consider the Γ-orbits in {1, 2, . . . , p}. Let {i 1 , . . . , i C } be a complete system of representatives of the Γ-orbits, and for each c = 1, . . . , C, let p c be the cardinality of the Γorbit through i c . The order N of Γ equals the least common multiple of p 1 , p 2 , . . . , p C and one has Γ = {id, σ, σ 2 , . . . , σ N -1 }. THEOREM 6. Let Γ = σ be a cyclic group of order N . For α = 0, 1, . . . , N 2 set

r * α = # { c ∈ {1, . . . , C} ; α p c is a multiple of N } d * α = 1 (α = 0 or N/2) 2 (otherwise).
Then we have

L = # { α ; r * α > 0 }, r = (r * α ; r * α > 0) and k = d = (d * α ; r * α > 0). For c = 1, . . . , C define v (c) 1 , . . . , v (c) pc ∈ R p by v (c) 1 := 1 p c pc-1 k=0 e σ k (ic) , v (c) 2β :
= 2 p c pc-1 k=0 cos 2πβk p c e σ k (ic) (1 ≤ β < p c /2), v (c) 2β+1 := 2 p c pc-1 k=0 sin 2πβk p c e σ k (ic) (1 ≤ β < p c /2), v (c) pc := 1 p c pc-1 k=0 cos(πk)e σ k (ic) (if p c is even).
THEOREM 7. The orthogonal matrix U Γ from Theorem 4 can be obtained by arranging column vectors {v

(c) k }, 1 ≤ c ≤ C, 1 ≤ k ≤ p c in the following way: we put v (c) k earlier than v (c ) k if (i) [k/2] pc < [k /2] p c , or (ii) [k/2] pc = [k /2] p c
and c < c , or

(iii) [k/2] pc = [k /2] p c
and c = c and k is even and k is odd.

Proofs of the above results are postponed to the Appendix. We shall see there that R(σ) acts on the 2-dimensional space spanned by v EXAMPLE 9. Let us consider σ = (1, 2, 3)(4, 5)(6) ∈ S 6 . The three Γ-orbits are {1, 2, 3}, {4, 5} and {6}. Set i 1 = 1, i 2 = 4, i 3 = 6. Then p 1 = 3, p 2 = 2, p 3 = 1. We have

N = 6. We count r * 0 = 3, r * 1 = 0, r * 2 = 1, r * 3 = 1, so that Z Γ Sym(3; R) ⊕ Herm(1; C) ⊕ Sym(1; R). According to Theorem 7, U Γ = v (1) 1 , v (2) 1 , v (3) 1 , v (1) 2 , v (1) 3 , v (2) 2 =         1/ √ 3 0 0 2/3 0 0 1/ √ 3 0 0 -1/6 1/ √ 2 0 1/ √ 3 0 0 -1/6 -1/ √ 2 0 0 1/ √ 2 0 0 0 1/ √ 2 0 1/ √ 2 0 0 0 -1/ √ 2 0 0 1 0 0 0        
.

Then we have (cf. ( 13))

U Γ • R(σ k ) • U Γ =   I 3 ⊗ B 0 (σ k ) B 2 (σ k ) B 3 (σ k )   , where B 0 (σ k ) = 1, B 2 (σ k ) = Rot( 2πk 3 ) ∈ GL(2; R) and B 3 (σ k ) = (-1) k .
Here Rot(θ) denotes the rotation matrix cos θ -sin θ sin θ cos θ for θ ∈ R.

The block diagonal decomposition of Z Γ is

U Γ • Z Γ • U Γ =      x 1 x 2 I 2 x 3   ; x 1 ∈ Sym(3; R), x 2 , x 3 ∈ R    .
2.7. The case of a non-cyclic group. In Section 2.5, we gave the general algorithm for determining the structure constants as well as the invariant measure ϕ Γ . In principle, the factorization of a determinant can be done in PYTHON, however there are some limitations regarding the dimension of a matrix. If the p × p matrix is not sparse, then the number of terms in the usual Laplace expansion of a determinant produces a polynomial with p! terms. The RAM memory requirements for calculating such a polynomial would be in excess of p!, which cannot be handled on a standard PC even for moderate p. Depending on the subgroup and the method of calculating the determinant, we were able to obtain the determinant for models of dimensions up to 10-20. In order to factorize the determinant for moderate to high dimensions, we want to find an orthogonal matrix U such that U • X • U is sparse enough for a computer to calculate its determinant Det U • X • U = Det (X). The matrix U Γ from (1) is in general very hard to obtain, but we propose an easy surrogate.

Let Γ be a subgroup of S p , which is not necessarily a cyclic group. Let {i 1 , i 2 , . . . , i C } be a complete family of representatives of the Γ-orbits in V = {1, . . . , p}. Take σ 0 ∈ S p for which the cyclic group Γ 0 := σ 0 generated by σ 0 has the same orbits in V as Γ does. In general, there are no inclusion relations between Γ and Γ 0 . However, we observe that a vector

v ∈ R p is Γ-invariant, (i.e. R(σ)v = v for all σ ∈ Γ) if and only if R(σ 0 )v = v if
and only if v is constant on Γ-orbits (i.e. v i = v j if i and j belong to the same orbit of Γ); see Example 10 below.

Let U Γ0 be the orthogonal matrix constructed as in Theorem 7 from the cyclic group Γ 0 . Note that the first C column vectors of U Γ0 are v

(1) 1 , v (2) 1 , . . . , v (C) 1 , which are Γ-invariant. The space V 1 := span v (c) 1 ; c = 1, . . . , C ⊂ R p is the trivial-representation-component of Γ as explained after (12). Therefore, if α 1 is the trivial representation of Γ, then r 1 = C and d 1 = k 1 = 1. The orthogonal complement V ⊥ 1 of V 1 is spanned by the rest of v (c) β , 1 ≤ c ≤ C, 1 < β ≤ 2[p c /2]. For X ∈ Z Γ , we see that X • v ∈ V 1 for v ∈ V 1 and that X • w ∈ V ⊥ 1 for w ∈ V ⊥ 1 . In other words, U Γ0 • X • U Γ0 = x 1 0 0 y (18)
with x 1 ∈ Sym(C; R) and y ∈ Sym(p -C; R). Actually, the correspondence φ 1 : Z Γ X → x 1 ∈ Sym(C; R) is exactly the Jordan algebra homomorphism defined before Corollary 5. Therefore Det (X) = Det (x 1 ) Det (y) , (19) while the factor Det (x 1 ) = det(φ 1 (X)) is an irreducible polynomial of degree r 1 = C. In this way, for any subgroup Γ, we are able to factor out the polynomial of degree equal to the number of Γ-orbits in V easily. On the other hand, the factorization of Det (y) requires study of the subrepresentation R of Γ on V ⊥ 1 , where the group Γ 0 is useless in general.

EXAMPLE 10. Let Γ = (1, 2, 3), (4, 5, 6) ⊂ S 6 , which is not a cyclic group. The space Z Γ consists of symmetric matrices of the form and moreover, Z Γ does not coincide with Z σ for any σ ∈ S 6 . Noting that the group Γ has two orbits: {1, 2, 3} and {4, 5, 6}, we define σ 0 := (1, 2, 3)(4, 5, 6). Taking i 1 = 1 and i 2 = 4, we have

X =        
U Γ0 =         1/ √ 3 0 2/3 0 0 0 1/ √ 3 0 -1/ √ 6 1/ √ 2 0 0 1/ √ 3 0 -1/ √ 6 -1/ √ 2 0 0 0 1/ √ 3 0 0 2/3 0 0 1/ √ 3 0 0 -1/ √ 6 1/ √ 2 0 1/ √ 3 0 0 -1/ √ 6 -1/ √ 2        
.

Note that the first two column vectors 1/ √ 3 1/ √ 3 1/ √ 3 0 0 0 and We cannot expect that the matrix y in ( 18) is always of a nice form, as in the example above. However, we note that in many examples we considered, the matrix y was sparse, which also makes the problem of calculating Det (X) much more feasible on a standard PC.

0 0 0 1/ √ 3 1/ √ 3 1/ √ 3 of U Γ0 are Γ-invariant. By direct calculation we verify that U Γ0 • X • U Γ0 is of the form         A B 0 0 0 0 B C 0 0 0 0 0 0 D 0 0 0 0 0 0 D 0 0 0 0 0 0 E 0 0 0 0 0 0 E         ,
In general Γ 0 defined above is not a subgroup of Γ. As we argue below, valuable insight about the factorization of Z Γ can be obtained by studying cyclic subgroups of Γ. In general, if Γ 1 is a subgroup of Γ, then Z Γ is a subspace of Z Γ1 .

Let Γ 1 be a cyclic subgroup of Γ and let U Γ1 be the orthogonal matrix constructed in Theorem 7. By Theorem 4, for any

X ∈ Z Γ ⊂ Z Γ1 the matrix U Γ1 • X • U Γ1 belongs to the space                        M K1 (x 1 ) ⊗ I k1 d1 M K2 (x 2 ) ⊗ I k2 d2 . . . M KL (x L ) ⊗ I kL dL         ; x i ∈ Herm(r i ; K i ) i = 1, . . . , L               
, where the structure constants (k i , d i , r i ) L i=1 can be easily calculated using Theorem 6. In particular, we have k 1 = d 1 = 1 and r 1 is the number of Γ 1 -orbits in {1, . . . , p}. Thus, we have M K1 (x 1 ) ⊗ I k1/d1 = x 1 ∈ Sym(r 1 ; R). In contrast to (18), x 1 in general can be further factorized and we know that Det (x 1 ) from ( 19) is an irreducible factor of Det (x 1 ). In conclusion, each cyclic subgroup of the general group Γ brings various information about the factorization. EXAMPLE 11. We continue Example 10. Let Γ 1 = (1, 2, 3) , which is a subgroup of Γ. There are four Γ 1 -orbits in V , that is, {1, 2, 3}, {4}, {5}, and {6}. We have

U Γ1 =         1/ √ 3 0 0 0 2/3 0 1/ √ 3 0 0 0 -1/ √ 6 1/ √ 2 1/ √ 3 0 0 0 -1/ √ 6 -1/ √ 2 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0         . For X ∈ Z Γ , we see that U Γ1 • X • U Γ1 is of the form         A 11 A 21 A 31 A 41 0 0 A 21 A 22 A 32 A 42 0 0 A 31 A 32 A 33 A 43 0 0 A 41 A 42 A 43 A 44 0 0 0 0 0 0 D 0 0 0 0 0 0 D        
, where A ij are linear functions of a, b, . . . , e, but they are not linearly independent. Indeed, we have

Det     A 11 A 21 A 31 A 41 A 21 A 22 A 32 A 42 A 31 A 32 A 33 A 43 A 41 A 42 A 43 A 44     = E 2 det A B B C ,
which exemplifies the fact that Det (x 1 ) is an irreducible factor of Det (x 1 ).

3. Gamma integrals.

Gamma integrals on irreducible symmetric cones.

Let Ω be one of the first three kinds of irreducible symmetric cones, that is, Ω = Herm + (r; K), where K ∈ {R, C, H}. Determinant and trace on corresponding Euclidean Jordan algebras will be denoted by det and tr. Then, we have the relation

dim Ω = r + r(r -1) 2 d, where d = 1 if K = R, d = 2 if K = C and d = 4 if K = H.
Let m(dx) denote the Euclidean measure associated with the Euclidean structure defined on A = Herm(r;

K) by x, y = tr[x • y] = tr[x • y]. The Gamma integral Γ Ω (λ) := Ω det(x) λ e -tr[x] det(x) -dim Ω/r m(dx) is finite if and only if λ > 1 2 (r -1)d = dim Ω/r -1 and in such case Γ Ω (λ) = (2π) (dim Ω-r)/2 Γ(λ)Γ(λ -d/2) . . . Γ(λ -(r -1)d/2). (20)
Moreover, one has

Ω det(x) λ e -tr[x•y] det(x) -dim Ω/r m(dx) = Γ Ω (λ) det(y) -λ (21)
for any y ∈ Ω.

Let us explain the role of the measure µ Ω (dx) = det(x) -dim Ω/r m(dx). Let G(Ω) be the linear automorphism group of Ω, that is, the set { g ∈ GL(A) ; g Ω = Ω }, where A is the associated Euclidean Jordan algebra. Then, the measure µ Ω is a G(Ω)-invariant measure in the sense that for any Borel measurable set B one has

µ Ω (g -1 B) = µ Ω (B) (g ∈ G(Ω)).
3.2. Gamma integrals on the cone P Γ . We endow the space Z Γ with the scalar product

x, y = Tr [x • y] (x, y ∈ Z Γ ).
Let dX denote the Euclidean measure on the Euclidean space (Z Γ , •, • ). Let us note that this normalization is not important in the model selection procedure as there we always consider quotients of integrals.

EXAMPLE 12. Consider p = 3 and Γ = S 3 . The space Z Γ is 2-dimensional and it consists of matrices of the form (see Example 1)

X =   a b b b a b b b a   for a, b ∈ R. Since X 2 = Tr X 2 = 3a 2 + 6b 2 = v v with v = ( √ 3a, √ 6b), we have dX = √ 3 √ 6 da db = 3 √ 2 da db.
Generally, if m i denotes the Euclidean measure on A i := Herm(r i ; K i ) with the inner product defined from the Jordan algebra trace (recall ( 8) and ( 9)), then (1) implies that for X ∈ Z Γ we have

X 2 = X, X = L i=1 k i d i Tr M Ki (x i ) 2 = L i=1 k i tr[(x i ) 2 ],
which implies that

(22) dX = L i=1 k i dim Ωi m i (dx i ) = L i=1 k dim Ωi/2 i m i (dx i ).
The key ingredient to compute the Gamma integral on P Γ is Theorem 4. Let Ω i denote the symmetric cone of the simple Jordan algebra A i = Herm(r i ; K i ) and

d i = dim R K i , i = 1, . . . , L. Recall that dim Ω i = r i + r i (r i -1)d i /2.
We start with the following more general integral. Recall that, for X ∈ Z Γ represented as in (1), we write φ i (X) = x i ∈ A i for i = 1, . . . , L.

LEMMA 8. For any Y ∈ P Γ and λ i > -1, i = 1, . . . , L, we have

PΓ L i=1 det(φ i (X)) λi e -Tr[Y •X] dX = e -BΓ L i=1 k -riλi i L i=1 Γ Ωi (λ i + dim Ω i /r i ) det(φ i (Y )) λi+dim Ωi/ri , ( 23
)
where

B Γ := 1 2 L i=1 (dim Ω i )(log k i ). ( 24
)
The proof is postponed to the Appendix. DEFINITION 13. Let G(P Γ ) = { g ∈ GL(p; R) ; gP Γ = P Γ } be the linear automorphism group of P Γ . We define the G(P Γ )-invariant measure ϕ Γ (X) dX by

ϕ Γ (X) = e BΓ P * Γ e -Tr[X•Z] dZ. PROPOSITION 9. We have ϕ Γ (X) = L i=1 det(φ i (X)) -dim Ωi/ri . (25)
PROOF. Recall that P Γ is a symmetric cone, so that it coincides with its dual cone, P * Γ . Thus,

ϕ Γ (Y ) = e BΓ PΓ e -Tr[Y •X] dX (Y ∈ P Γ ).
Setting λ 1 = . . . = λ L = 0 in ( 23) we obtain the expression of ϕ Γ (Y ) above.

DEFINITION 14. The Gamma function of P Γ is defined by the following integral Tr[X] ϕ Γ (X) dX, whenever it converges.

(26) Γ PΓ (λ) := PΓ Det (X) λ e -
THEOREM 10. The integral (26) converges if and only if

λ > max i=1,...,L (r i -1)d i 2k i (27)
and, for these values of λ, we have

Γ PΓ (λ) = e -AΓλ+BΓ L i=1 Γ Ωi (k i λ) (28)
where Γ Ωi is given in (20), B Γ in (24) and

A Γ := L i=1 r i k i log k i . (29) Moreover, if Y ∈ P Γ and (27) holds true, then PΓ Det (X) λ e -Tr[Y •X] ϕ Γ (X) dX = Γ PΓ (λ)Det (Y ) -λ . ( 30
)
We also have the following result

THEOREM 11. If Y ∈ P Γ and λ > max i=1,...,L - 1 k i , (31) then PΓ Det (X) λ e -Tr[Y •X] dX = e -AΓλ-BΓ L i=1 Γ Ωi k i λ + dim Ω i r i Det (Y ) -λ ϕ Γ (Y ). (32)
PROOF OF THEOREM 10 AND THEOREM 11. Recall that for X ∈ P Γ we have Det (X) = L i=1 det(φ i (X)) ki , where the map φ i : Z Γ → Herm(r i ; K i ) is a Jordan algebra homomorphism, i = 1, . . . , L.

If λ i = k i λ -dim Ω i /r i and ( 27) holds, then (23) implies

L i=1 k -riλi i = e -AΓλ+2BΓ and L i=1 det(φ i (Y )) -λi-dim Ωi/ri = L i=1 det(φ i (Y )) ki -λ
If λ i = k i λ with (31), then by ( 23) we obtain

L i=1 k -riλi i = e -AΓλ and L i=1 det(φ i (Y )) -λi-dim Ωi/ri = L i=1 det(φ i (Y )) ki -λ ϕ Γ (Y ).
3.3. RCOP-Wishart laws on P Γ . Let Σ ∈ P Γ ⊂ Sym + (p; R) and consider i.i.d. random vectors Z (1) , . . . , Z (n) following the N p (0, Σ) distribution. Define U i = Z (i) • Z (i) , i = 1, . . . , n, and U = n i=1 U i . Our aim is to analyse the probability distribution of the random matrix

W n = π Γ (U ) = π Γ (U 1 + • • • + U n ) = π Γ (U 1 ) + . . . + π Γ (U n ).
In the rest of this section, we find n 0 such that for n ≥ n 0 the random matrix W n follows an absolutely continuous law, and we compute its density. Further, we extend the shape parameter to a continuous range and define the RCOP-Wishart law on P Γ .

We start with the following easy result.

LEMMA 12. For any θ ∈ Sym + (p; R) we have

Ee -Tr[θ•πΓ(U1)] = Det (I p + 2 Σ • π Γ (θ)) -1/2 .
PROOF. Using (5) repeatedly we have

Tr [θ • π Γ (U 1 )] = Tr [π Γ (θ) • π Γ (U 1 )] = Tr [π Γ (θ) • U 1 ] .
The assertion follows from the usual multivariate Gauss integral.

PROPOSITION 13. The law of W n is absolutely continuous on P Γ if and only if

n ≥ n 0 := max i=1,...,L r i d i k i . ( 33
)
If n ≥ n 0 , then its density function with respect to dX is given by

Det (X) n/2 e -1 2 Tr[X•Σ -1 ] Det (2Σ) n/2 Γ PΓ ( n 2 ) ϕ Γ (X)1 PΓ (X). (34) PROOF. With λ = n/2, condition (27) becomes n > max i=1,...,L (r i -1)d i k i .
Since the quotient k i /d i is an integer, the last condition is equivalent to (33).

In view of Lemma 12, it is enough to show that W n has density (34) if and only if for any θ ∈ P Γ ,

Ee -Tr[θ•Wn] = Det (I p + 2 Σ • θ) -n/2 . ( 35
)
This follows directly from (30).

It is known that the MLE exists and is unique if and only if the sufficient statistic lies in the interior of its convex support, see [START_REF] Barndorff-Nielsen | Information and exponential families in statistical theory[END_REF]. It is clear that if (33) is not satisfied, then the support of W n is contained in the boundary of P Γ . COROLLARY 14. If the number of samples n satisfies (33), then the MLE of Σ exists and is given by

Σ = 1 n π Γ (U 1 + • • • + U n ).
The above result has been already proven in (Andersson, 1975, Theorem 5.9) (see also (Andersson and Madsen, 1998, Sec. A.3, A.4)). Moreover, in (Andersson and Madsen, 1998, Sec. A.5) a formula for E Det Σ α has been given. The right hand side of their formula (A.4) coincides with

2 αp ΓP Γ (α+ n 2 ) n αp ΓP Γ ( n
2 ) in our notation (substitute (d i , nk i /d i , r i ) i for their (d µ , n µ , p µ ) µ ). However, their formula does not imply an explicit expression for I Γ (δ, D).

Let us recall that the MLE of Σ in the standard normal model exists if and only if n ≥ p. We recover this case for Γ = {id}, since then we have L = 1, r 1 = p and k 1 = d 1 = 1.

When n < n 0 , the law of W n is singular, and it can be described as a direct product of the singular Wishart laws on the irreducible symmetric cones Ω i , see e.g. [START_REF] Hassairi | Riesz exponential families on symmetric cones[END_REF]. DEFINITION 15. Let η > max (r i -1) di ki ; i = 1, . . . , L and Σ ∈ P Γ . The RCOP-Wishart law W Γ η,Σ is defined by its density

W Γ η,Σ (dX) = Det (X) η/2 e -1 2 Tr[X•Σ -1 ] Det (2Σ) η/2 Γ PΓ ( η 2 ) ϕ Γ (X)1 PΓ (X) dX. (36)
With this new notation, we see that if (33) is satisfied, then W n ∼ W Γ n,Σ .

LEMMA 15. The Jacobian of the transformation

P Γ X → X -1 ∈ P Γ (37) equals ϕ Γ (X -1 ) 2 .
Proof of the lemma is postponed to the Appendix. The lemma gives the following result.

PROPOSITION 16. Let W ∼ W Γ η,Σ with η > max { (r i -1)d i /k i ; i = 1, . . . , L } and Σ ∈ P Γ . Then its inverse Y = W -1 has density Det (y) -η/2 e -1 2 Tr[y -1 •Σ -1 ] Det (2Σ) η/2 Γ PΓ ( η 2 )
ϕ Γ (y)1 PΓ (y).

3.4. The Diaconis-Ylvisaker conjugate prior for K. The Diaconis-Ylvisaker conjugate prior [START_REF] Diaconis | Conjugate priors for exponential families[END_REF]) for the canonical parameter K = Σ -1 is given by

f (K; δ, D) = 1 I Γ (δ, D) Det (K) (δ-2)/2 e -1 2 Tr[K•D] 1 PΓ (K) dK,
for hyper-parameters δ > 2 max { 1 -1/k i ; i = 1, . . . , L } and D ∈ P Γ . By (32), the normalizing constant is equal to

(38) I Γ (δ, D) = e -AΓ(δ-2)/2-BΓ L i=1 Γ Ωi k i δ -2 2 + dim Ω i r i Det (D) -(δ-2)/2 ϕ Γ (D),
where A Γ , B Γ and ϕ Γ are given in ( 29), ( 24) and (25).

4. Model selection. Bayesian model selection on all colored spaces seems at the moment intractable. This is due in great part to a poor combinatorial description of the colored spaces Z Γ . In particular, the number of such spaces, that is, # { Z Γ ; Γ ∈ S p } is generally unknown for large p. It was shown in [START_REF] Gehrmann | Lattices of graphical Gaussian models with symmetries[END_REF] that these colorings constitute a lattice with respect to the usual inclusion of subspaces. However the structure of this lattice is rather complicated and is unobtainable for big p. This, in turn, does not allow to define a Markov chain with known transition probabilities on such colorings. Finally, the fundamental problem which prevents us from doing Bayesian model selection on all colored spaces for arbitrary p is the following. In order to compute Bayes factors, one has to be able to find the structure constants (k i , d i , r i ) L i=1 for arbitrary subgroups of S p . This is equivalent to finding irreducible representations over reals for an arbitrary finite group, which is very hard in general, although general algorithms have been developed for this issue (see [START_REF] Plesken | Constructing rational representations of finite groups[END_REF]).

In this section, we are making a step forward in the problem of model selection for colored models in two ways. In Section 4.1, we use the results of Section 2.6, to obtain the structure constants when we restrict our search to the space of colored models generated by a cyclic group, that is, when Γ = σ for σ ∈ S p and we propose a model selection procedure restricted to the cyclic colorings. In Section 4.2, we use Remarks 7 and 8 to obtain the irreducible representations of Z Γ and the structure constants by factorization of the determinant. We apply this technique to do model selection for the four-dimensional example given by Frets' data since, in that case, there are only 22 models and we can compute all the Bayes factors.

4.1. Model selection within cyclic groups. The smaller space of cyclic colorings has a much better combinatorial description. In particular, the following result can be proved.

LEMMA 17. If Z σ = Z σ for some σ, σ ∈ S p , then σ = σ .
This result allows us to calculate the number of different colorings corresponding to cyclic groups, that is, the number of labeled cyclic subgroups of the symmetric group S p , which can be found in OEIS, sequence A051625 (see the last column of Table 1).

We will present two applications of the Metropolis-Hastings algorithm. In the first one, the Markov chain will move on the space of cyclic groups. The drawback of this first approach is that we need to compute the proposal distribution g, whose computational complexity grows faster than quadratically as p increases (see ( 40)). In the second algorithm, we consider a larger state space S p , which allows us to consider an easy proposal distribution. However, this comes at the cost of slower convergence of the posterior probabilities (see Theorem 18). 4.1.1. First approach. Each cyclic subgroup Γ can be uniquely represented by a permutation, which is minimal in the lexicographic order within permutations generating Γ. Let ν(Γ) ∈ S p be such a permutation, that is,

ν(Γ) = min { σ ∈ S p ; σ = Γ } . Define c t := ν(c t-1 ) • x t , (39)
where c 0 is a fixed cyclic subgroup and (x t ) t∈N is a sequence of i.i.d. random transpositions distributed uniformly, that is, P(x t = α) = 1/ p 2 for any α ∈ T := {(i, j) ∈ S p }. Clearly, the sequence (c t ) t is a Markov chain. Its state space is the set of all cyclic subgroups of S p . Moreover, the trivial subgroup {id} can be reached from any subgroup c t (and vice versa) in a finite number of steps with positive probability. Thus the chain (c t ) t is irreducible. The proposal distribution in the Metropolis-Hastings algorithm is the conditional distribution of c t |c t-1 . It is proportional to the number of possible transitions from c to c , that is,

g c |c := # { (i, j) ∈ S p ; c = ν(c) • (i, j) } p 2 , ( 40
)
where c and c are cyclic subgroups.

We follow the principles of Bayesian model selection for graphical models, presented, for example, in (Maathuis et al., 2018, Chapter 10, p.247). Let Γ be uniformly distributed on the set C := { σ ; σ ∈ S p } of cyclic subgroups of S p . We assume that K|{Γ = c}, c ∈ C, follows the Diaconis-Ylvisaker conjugate prior distribution on P c with hyper-parameters δ and D, that is,

f K|Γ=c (k) = 1 I c (δ, D) Det (k) (δ-2)/2 e -1 2 Tr[D•k] 1 Pc (k),
where the normalizing constant is given in (38). Suppose that Z 1 , . . . , Z n given {K = k, Γ = c} are i.i.d. N p (0, k -1 ) random vectors with k ∈ P c . Then, it is easily seen that we have

P (Γ = c|Z 1 , . . . , Z n ) ∝ I c (δ + n, D + U ) I c (δ, D) (c ∈ C) (41) with U = n i=1 Z i • Z i .
Φ(c) = #C c (c ∈ C).
For c ∈ C, we denote

π c = P(Γ = c|Z 1 , . . . , Z n ),
which we want to approximate. In our model we have (see ( 41))

π c ∝ I c (δ + n, D + U ) I c (δ, D) (c ∈ C). (43)
In order to find π = (π c ; c ∈ C) let us consider π = (π σ ; σ ∈ S p ), a probability distribution on S p such that

πσ ∝ I σ (δ + n, D + U ) I σ (δ, D) (σ ∈ S p ). (44)
Since ( 43) and ( 44) imply that πσ ∝ π σ , we have

πσ = π σ c∈C Φ(c)π c (σ ∈ S). ( 45
)
As before, let (x t ) t∈N be a sequence of i.i.d random transpositions distributed uniformly on T = {(i, j) ∈ S p }. We define a random walk on S p by

s t+1 = s t • x t+1 , (t = 0, 1, . . .).
Then, (s t ) t is an irreducible Markov chain with symmetric transition probability

g(σ |σ) = 1 ( p 2 ) , if σ -1 • σ ∈ T , 0, if σ -1 • σ / ∈ T .
We note that ( s t ) t is not a Markov chain on the space of cyclic subgroups. Indeed, it can be shown that the necessary conditions for (f (s t )) t to be a Markov chain (see (Burke and Rosenblatt, 1958, Eq. (3))) are not satisfied for f (σ) := σ if p > 4. A remedy for this fact was introduced in (39). Indeed, the sequence ( s t ) t is very similar to the sequence (c t ) t defined previously. Both move along cyclic subgroups and their definitions are very similar. However, ( s t ) t is not a Markov chain, whereas (c t ) t is a Markov chain. We took care of this problem by using the minimal generator ν(•) as in definition (39) of c t . We use the Metropolis-Hastings algorithm with the above proposal distribution to approximate π.

ALGORITHM 17. Starting from a permutation σ 0 ∈ S p , repeat the following two steps for t = 1, 2, . . .:

1. Sample x t uniformly from the set T of all transpositions and set σ = σ t-1 • x t ; 2. Accept the move σ t = σ with probability min 1,

I σ (δ + n, D + U ) I σt-1 (δ, D) I σ (δ, D) I σt-1 (δ + n, D + U ) .
If the move is rejected, set σ t = σ t-1 .

By the ergodicity of the Markov chain (σ t ) t constructed above, as the number of steps T → ∞, we have

T t=1 1 σ=σt T a.s. -→ πσ (σ ∈ S p ). ( 46
)
This fact allows us to develop a scheme for approximating the posterior probability π.

THEOREM 18. We have as T → ∞, -→ πσ . Moreover,

1 Φ(c) T t=1 1 c= σt T t=1 1 Φ( σt ) = 1 Φ(c) σ∈Cc n (T ) σ T t=1 γ∈C 1 Φ(γ) 1 γ= σt = 1 Φ(c) σ∈Cc n (T ) σ T γ∈C 1 Φ(γ) σ∈Cγ n (T ) γ T a.s.
-→

1 Φ(c) σ∈Cc πσ γ∈C 1 Φ(γ) σ∈Cγ πγ .
Finally, by (45) we have

1 Φ(c) σ∈Cc πσ = π c γ∈C Φ(γ)π γ ∝ π c ,
which completes the proof.

In order to approximate the posterior probability π, we allowed the Markov chain to travel on the larger space S p . In particular, each state c ∈ C was multiplied Φ(c) ≥ 1 times, where Φ(c) is the number of permutations generating c. This procedure should result in slower convergence to the stationary distribution in (47). By comparing with (46), we see that (47) can be interpreted as follows: let us assign to each cyclic subgroup c a weight 1/Φ(c) ≤ 1. Then, the denominator N T := T t=1 1/Φ( σ t ) can be thought of as an "effective" number of steps and the numerator is the number of "effective" steps spent in state c. In general, for large T we expect N T T (see an example in Section 5.2).

4.2. Model selection for p = 4. Our numbering of colored models on four vertices is in accordance with (Gehrmann, 2011, Figures 15 and 16, p 674-675). However, we identify models by the largest group with the same coloring Γ * rather than the smallest as in [START_REF] Gehrmann | Lattices of graphical Gaussian models with symmetries[END_REF]. There are 30 different subgroups of S 4 , which generate 22 different colored spaces. Up to conjugacy (renumbering of vertices), there are 8 different conjugacy classes. Within a conjugacy class, the sets of constants {k i , r i , d i } L i=1 remain the same. Groups Γ * k for k = 1, . . . , 17 correspond to cyclic colorings. We apply our results and methods in order to do Bayesian model selection for the celebrated example of Frets' heads, Frets (1921); [START_REF] Whittaker | Graphical models in applied multivariate statistics[END_REF]. The head dimensions (length L i and breadth B i , i = 1, 2) of 25 pairs of first and second sons were measured. Thus we have n = 25 and p = 4. The following sample covariance matrix is obtained (we have 2287.04 1268.84 1671.88 1106.68 1268.84 1304.64 1231.48 841.28 1671.88 1231.48 2419.36 1356.96 1106.68 841.28 1356.96

Z = (L 1 , B 1 , L 2 , B 2 ) ), U = n i=1 Z (i) • Z (i) =    
1080.56     .
We perform Bayesian model selection within all RCOP models, not just the ones corresponding to cyclic subgroups. In Table 2 we list all RCOP models on full graph with four vertices, along with corresponding structure constants. Structure constants remain the same within a conjugacy class, however the invariant measure ϕ Γ is always different. Since there are only 22 such models, we calculate all exact posterior probabilities. The Table 2 and the invariant measures ϕ Γ were obtained by using Remark 7.

(ki) (ri) (di) Γ * 1 = {id} (1) (4) (1) Γ * 2 = (1, 2) (1,1) (3,1) (1,1) Γ * 3 = (1, 3) Γ * 4 = (1, 4) Γ * 5 = (2, 3) Γ * 6 = (2, 4) Γ * 7 = (3, 4) Γ * 8 = (1, 2, 3), (1, 2) (1,2) (2,1) (1,1) Γ * 9 = (1, 2, 4), (1, 2) Γ * 10 = (1, 3, 4), (1, 3) Γ * 11 = (2, 3, 4), (2, 3) Γ * 12 = (1, 2)(3, 4) (1,1) (2,2) (1,1) Γ * 13 = (1, 3)(2, 4) Γ * 14 = (1, 4)(2, 3) Γ * 15 = (1, 2, 3, 4), (1, 3) (1,1,2) (1,1,1) (1,1,1) Γ * 16 = (1, 2, 4, 3), (1, 4) Γ * 17 = (1, 3, 2, 4), (1, 2) Γ * 18 = (1, 2), (3, 4) (1,1,1) (2,1,1) (1,1,1) Γ * 19 = (1, 3), (2, 4) Γ * 20 = (1, 4), (2, 3) Γ * 21 = (1, 2)(3, 4), (1, 4)(2, 3) (1,1,1,1) (1,1,1,1) (1,1,1,1) Γ * 22 = S4 (1,3) (1,1) (1,1)
In Table 3 we summarize the results when δ = 3, giving the three best coloring models with the highest posterior probability, for each given D. Results are very similar for δ = 10 and the given values of D.

For different values of D = dI 4 , the only models that have highest posterior probability are the 4 models:

Γ * 22 = S 4 , Γ * 19 = (1, 3), (2, 4) , Γ * 13 = (1, 3)(2, 4) , Γ * 1 = {id}. These four subgroups form a path in the Hasse diagram of subgroups of S * 4 , i.e. Γ * 22 ⊃ Γ * 19 ⊃ Γ * 13 ⊃ Γ *
1 . Thus the four selected colorings, corresponding to the permutation groups are in some way consistent. Moreover, each of them has a good statistical interpretation. Let us interpret models Γ * 13 and Γ * 19 . Recall the enumeration of vertices (1, 2, 3, 4) = (L 1 , B 1 , L 2 , B 2 ). The invariance with respect to the transposition (1, 3) means that L 1 is exchangeable with L 2 and, similarly, the invariance with respect to the transposition (2, 4) implies exchangeability of B 1 and B 2 . Both together correspond to the fact that sons should be exchangeable in some way.

We observe that only the Γ * 22 model appeared in former attempts of model selection for Frets' heads data. It was selected in [START_REF] Højsgaard | Graphical Gaussian models with edge and vertex symmetries[END_REF], by a likelihood ratio test. Note that the only complete RCOP model selected in Gehrmann (2011) (who used the Edwards-Havranek model selection procedure) among the 9 minimally accepted models on p. 676 of her article is Γ * 10 , which is not selected by our exact Bayesian procedure for any choice 5. Simulations. Let the covariance matrix Σ be the symmetric circulant matrix

Σ =        c 0 c 1 . . . c 2 c 1 c 1 c 0 c 1 . . . c 2 . . . c 1 c 0 . . . . . . c 2 . . . . . . c 1 c 1 c 2 . . . c 1 c 0       
, with c 0 = 1 + 1/p and c k = 1 -k/p for k = 1, . . . , p/2 . It is easily seen that this matrix belongs to P σ * with σ * = (1, 2, . . . , p -1, p).

5.1. First approach. For p = 10 and n = 20, we sampled Z (1) , . . . , Z (n) from the N p (0, Σ) distribution and obtained U = n i=1 Z (i) • Z (i) depicted in Fig. 1 (b). We run the Metropolis-Hastings algorithm starting from the group σ 0 = {id} with hyperparameters δ = 3 and D = I 10 . After 1 000 000 steps, the five most visited states are given in the Tab. 4.

The Metropolis-Hastings (M-H) algorithm recovered the true pattern of the covariance matrix. The acceptance rate was 2.5% and the Markov chain visited 746 different cyclic groups. The acceptance rate can be increased by a suitable choice of the hyper-parameters (e.g. for D = 10I 10 the acceptance rate is around 10%).

In order to grasp how randomness may influence results, we performed 100 simulations, where each time we sample Z (1) , . . . , Z (n) from N p (0, Σ) and we run M-H for 100 000 steps with the same parameters as before. In Table 5 we present how many times a given cyclic subgroup was most visited during these 100 simulations (second column). There were 53 distinct cyclic subgroups, which were most visited at least in one of the 100 simulations; below we present 10 such subgroups. The average acceptance rate is 1.4% (see the histogram (1,2,3,4,5,6,7,8,9,10) 25 1.00 (1, 3, 5, 7, 9)(2, 4, 6, 8,10) 13 0.60 (1,2,4,3,5,6,7,9,8,10) 3 0.43 (1,2,4,3,5,6,7,8,9,10) 2 0.46 (1,3,2,4,5,6,8,7,9,10) 2 0.43 (1,3,5,9,2,6,8,10,4,7) 2 0.43 (1,4,3,5,2,6,9,8,10,7) 2 0.35 (1,4,5,7,8)(2,3,6,9,10) 2 0.24 (1,8,10,9)(2,7)(3,5,4,6 in Fig. 2). When we regard colorings as partitions of the set V ∪ E according to group orbit decomposition, the two colorings may be compared using the so-called adjusted Rand index (ARI, see [START_REF] Hubert | Comparing partitions[END_REF]), a similarity measure comparing partitions which takes values between -1 and 1, where 1 stands for perfect match and independent random labelings have score close to 0. In the third column of Table 5, we give the adjusted Rand index between the colorings generated by given cyclic subgroup and the true coloring. We see that groups which were most visited by the Markov chain have positive ARI and the true pattern was recovered in a quarter of cases. We stress that even though the colorings generated by (1, 2, 3, 4, 5, 6, 7, 8, 9, 10) and (1, 3, 5, 7, 9)(2, 4, 6, 8, 10) are very similar, the distance between these subgroups is 9, that is, the Markov chain (C t ) t needs at least 9 steps to get from one subgroup to the other.

This indicates that the Markov chain may encounter many local maxima and one should always tune the hyper parameters in order to have higher acceptance rate or to allow the Markov chain (C t ) t to make bigger steps. 5.2. Second approach. We performed T = 100 000 steps of Algorithm 17 with σ 0 = id, p = 100, n = 200, δ = 3 and D = I 100 . Let us note that for p = 100, there are about 4 • 10 155 cyclic subgroups and this is the number of models we consider in our model search.

We have used Theorem 18 to approximate the posterior probability distribution (π c ; c ∈ C) (see ( 41)). The highest estimated posterior probability was obtained for c * := σ * , where σ * =(1, 2, 3, 4) (6, 8, 15)(7, 10, 9)(11, 16, 12)(13, 17, 14)(18, 19, 20, 22, 21)(23, 26) (24, 42, 28, 44)(25, 31, 30, 32)(27, 34)(29, 37)(33, 45)(35, 39, 36, 40) (38, 47, 41, 48)(43, 51, 46, 49)(50, 52, 53, 54)(56, 58, 57)(59, 66, 67) (60, 65, 63)(61, 62, 64)(68, 71, 72, 70, 69)(73, 93)(74, 77)(75, 98, 81, 100) is depicted, which can be thought of as an "effective" number of steps of the algorithm (for an explanation, see the paragraph at the end of Subsection 4.1.2). In blue, we present a sequence , which represents the number of weighted accepted steps, where the weight of the kth step equals 1 Φ( σk ) . We restricted the plot to steps k = 1, . . . , 10 000, because after 10 000 steps, the Markov chain (σ t ) 10 000≤t≤100 000 changed its state only 9 times. For k = 100 000, the value of the blue curve is 25.75, while the value of red one is 6 381.5.

The model suffers from poor acceptance rate, which could be improved by an appropriate choice of the hyper-parameter D or by allowing the Markov chain to do bigger steps.

SUPPLEMENTARY MATERIAL

Appendix: Proofs of Theorems 4, 6 and 7. Example for Remark 7.

For the proof of Theorem 4, we will need the following Lemma.

LEMMA 19. Let i, j = 1, 2, . . . , L, and assume that Y ∈ Mat(r i k i , r j k j ; R) satisfies the condition

Y L1 Y L2 • • • Y LL     
with Y ij ∈ Mat(r i k i , r j k j ; R). Then x := U Γ • y • U Γ belongs to Z Γ , so that (4) implies

R(σ) • U Γ • y • U Γ • R(σ) = U Γ • y • U Γ
for σ ∈ Γ, and this equality is rewritten as

[U Γ • R(σ) • U Γ ] • y = y • [U Γ • R(σ) • U Γ ].
By (13), we have

[I ri ⊗ B i (σ)] • Y ij = Y ij • [I rj ⊗ B j (σ)].
Lemma 19 tells us that Y ij = 0 if i = j, and that Y ii = M Ki (x i ) ⊗ I ki/di with some x i ∈ Mat(r i , r i ; K i ). Since y is a symmetric matrix, the block Y ii is also symmetric, which implies that x i ∈ Herm(r i ; K i ). Hence, the map ι : L i=1 Herm(r i ; K i ) (x i ) L i=1 → X ∈ Z Γ given by (15) gives a Jordan algebra isomorphism. EXAMPLE 18. In this example we present a colored space Z Γ ⊂ Sym(16; R), which has a component Herm(2; H). Let Γ = σ 1 , σ 2 be the subgroup of S 16 generated by the two permutations

  a > 0 and b ∈ (-a/2, a)

For A ∈

 ∈ Mat(m, n; R) and B ∈ Mat(m , n ; R), we denote by A⊕B the matrix A O O B ∈ Mat(m + m , n + n ; R), and by A ⊗ B the Kronecker product of A and B. For a positive integer r, we write B ⊕r for I r ⊗ B ∈ Mat(rm , rn ; R)

  rotation with the angle 2πβ/p c . The condition (i) means that the angle for v (c) k is smaller than the one for v (c ) k .

  where A, B, • • • , E are linear functions of a, b, • • • , e. The matrices x 1 and y are A . The matrix y is of such simple form, because Γ 0 is a subgroup of Γ in this case.

  These derivations allow us to run the Metropolis-Hastings algorithm restricted to cyclic groups, as follows.ALGORITHM 16. Starting from a cyclic group C 0 ∈ C, repeat the following two steps for t = 1, 2, . . .:1. Sample x t uniformly from the set T of all transpositions and set c = ν(C t-1 ) • x t ; 2. Accept the move C t = c with probability min 1,I c (δ + n, D + U ) I Ct-1 (δ, D) I c (δ, D) I Ct-1 (δ + n, D + U ) g (C t-1 |c ) g (c |C t-1 )If the move is rejected, set C t = C t-1 .4.1.2. Second approach. It is known that σ = σ if and only if σ = σ k for some k ∈ β(|σ|), whereβ(n) = { k ∈ {1, . . . , n} ; k and n are relatively prime } (42) and |σ| denotes the order of σ. Let C = { σ ; σ ∈ S p } denote the set of cyclic subgroups of S p . For c ∈ C we define Φ(c) := #β(|c|) and C c := { σ ∈ S p ; σ = c }, the set of permutations, which generate the cyclic subgroup c. We have

  , σ ∈ S p . We have T = σ∈Sp n

FIG 1 .

 1 FIG 1. Heat map of matrix Σ (a) and matrix U/n (b).

FIG 2 .

 2 FIG 2. Histogram of acceptance rates in 100 simulations of Metropolis-Hastings algorithm.

FIG 3 .

 3 FIG 3. Heat map of matrix Σ (a) and matrix U/n (b) and projection of U/n onto Zc * .

FIG 4 .

 4 FIG 4. Number of "effective" steps (red) and number of "effective" accepted steps (blue).

(

  76, 84, 78, 83)(79, 85)(80, 94, 82, 91)(86, 92, 87, 90)(88, 96, 89, 97)(95, 99). The order of c * is |c * | = 60 and Φ(c * ) = 16. The estimate of the posterior probability π c * is equal to (recall (47))The true covariance matrix Σ, the data matrix U/n and the projection Π c * (U/n) are illustrated in Fig.3.We visualize the performance of the algorithm on Fig 4. In red color, a sequence k t=1 1 Φ( σt ) k

  ) 1 σt = σt-1 k

Y,Y

  ⊗ B i (σ)] • Y = Y • [I rj ⊗ B j (σ)] (σ ∈ Γ).If i = j, then there exists C ∈ Mat(r i , r i ;K i ) such that Y = M Ki (C) ⊗ I ki/di. On the other hand, if i = j, then Y = 0. PROOF. Let us consider a block decomposition of Y as 11 Y 12 . . . Y 1,rj Y 21 Y 22 . . . Y 2,rj . . . . . . . . . . . . Y ri,1 Y ri,2 . . . Y ri,rjwhere each Y ab is a k i × k j matrix. Then (48) implies that(49) B i (σ) • Y ab = Y ab • B j (σ) (σ ∈ Γ) for all a, b. If i = j, then Y ab ∈ End Γ (R ki ), so that there exists C ab ∈ K i for which Y ab = M Ki (C ab ) ⊗ I ki/di thanksto Lemma 3. Let us consider the case i = j. Eq. (49) tells us that Ker Y ab ⊂ R kj is a Γ-invariant subspace, which then equals {0} or R kj because of the irreducibility of B j . Similarly, since Image Y ab ⊂ R ki is a Γ-invariant subspace by (49), Image Y ab equals {0} or R ki . Now suppose that Y ab = 0. Then Ker Y ab = {0} and Image Y ab = R ki by the argument above, and it means that Y ab induces an isomorphism from (B j , R kj ) onto (B i , R ki ). But this contradicts the fact that the representations B i and B j are not equivalent for i = j. Hence we get Y ab = 0. PROOF OF THEOREM 4. Take y ∈ U Γ • Z Γ • U Γ and consider the block decomposition of y as 11 Y 12 . . . Y 1L Y 21 Y 22 . . . Y 2L .

TABLE 1

 1 Number of all subgroups of a symmetric group, number of their conjugacy classes, number of different colorings and a number of cyclic groups p #subgroups of Sp #conjugacy classes of Sp #different ZΓ #cyclic groups

	1	1	1	1	1
	2	2	2	2	2
	3	6	4	5	5
	4	30	11	22	17
	5	156	19	93	67
	6	1 455	56	739	362
	7	11 300	96	4 508	2039
	8	151 221	296	?	14 170
	9	1 694 723	554	?	109 694
	10	29 594 446	1 593	?	976 412
	18	≈ 7.6 • 10 18	7.3 • 10 6	?	≈ 7.1 • 10 14

TABLE 2

 2 Structure constants for all colorings with four vertices

	Group

TABLE 3

 3 Posterior probabilities in Frets' heads for three best models, δ = 3 and given D.

	D	Best model	2nd best	3rd best
	I4 Γ * 22	(95.2%) Γ * 16	(2.5%)	Γ * 17	(1.3%)
	50I4 Γ * 19	(33.8%) Γ * 13	(29.6%) Γ * 8	(13.3%)
	100I4 Γ * 13	(39.6%) Γ * 19	(29.8%) Γ * 8	(7.2%)
	1000I4 Γ * 1	(38.9%) Γ * 13	(10.5%) Γ * 3	(10.3%)

TABLE 5

 5 Cyclic subgroups which were chosen by M-H algorithm most often generator of a cyclic group #most visited ARI

The On-Line Encyclopedia of Integer Sequences, https://oeis.org/.

The number of subgroups of Sp is unknown for p > 18, see[START_REF] Holt | Enumerating subgroups of the symmetric group[END_REF] and OEIS sequence A005432.

Acknowledgements. The authors would like to thank Steffen Lauritzen for his interest and encouragements.

* Supported by JSPS KAKENHI Grant Number 16K05174 and JST PRESTO. † Supported by grant 2016/21/B/ST1/00005 of the National Science Center, Poland. ‡ Supported by an NSERC Discovery Grant. MSC 2010 subject classifications: Primary 62H99, 62F15 ; secondary 20C35

Thus, Remark 7 gives us that L = 5 and

This in turn implies

As a matter of fact, the group Γ has four 1-dimensional representations and one 4dimensional real irreducible representation, and each representation appears twice in R 16 .

PROOF OF THEOREM 6 AND THEOREM 7. Let M := N 2 and denote the irreducible representations of Γ by

where Rot(θ) denotes the rotation matrix cos θ -sin θ sin θ cos θ for θ ∈ R. Then all the equivalence classes of the irreducible representations of Γ are

2 (otherwise). Recall that {i 1 , . . . , i C } is a complete system of representatives of the Γ-orbits, and, for each c = 1, . . . , C, p c is the cardinality of the Γ-orbit through i c . Let

where we have used σ pc (i c ) = i c and ζ pc c = 1 at the second equality. It follows that

with β pc = α N . Similarly, we have

pc

(if p c and N are even).

Therefore, for α = 0, . . . , [N/2], the multiplicity r α of the representation B α of Γ in (R, R p ) is equal to the number of c such that β pc = α N with some β ∈ N. In other words, (51)

Then we have

PROOF OF LEMMA 8. Denote the left hand side of (23) by I. Let us change variables x i = φ i (X) for i = 1, . . . , L. By ( 16) and ( 22) we obtain

Each integral can be calculated using (21) for λ i > -1 and φ i (Y ) ∈ Ω i , i = 1, . . . , L. Hence,

and so we obtain (23).

PROOF OF LEMMA 15. First observe that

so that, the Jacobian of (37) equals Det End (P X -1 ), where Det End is the determinant in the space of endomorphisms of Z Γ and for any X ∈ Z Γ by P X we denote the linear map on Z Γ to itself defined by P X Y = X • Y • X. It is easy to see that for any X ∈ P Γ we have P X ∈ G(P Γ ). Indeed, since P X Y is positive definite for Y ∈ P Γ , it is enough to verify that

This follows quickly by the fact that X, Y ∈ P Γ . Further, by the G(P Γ ) invariance of ϕ Γ , we have ϕ Γ (gX) = |Det End (g)| -1 ϕ Γ (X) (g ∈ G(P Γ )).

Taking g = P X -1 , we eventually obtain

where the latter inequality can be easily verified by ( 25).

().