Generating random correlation matrices based on partial correlation vines and the onion method
Harry Joe
Department of Statistics
University of British Columbia

Abstract
Partial correlation vines and the onion method are presented for generating random correlation matrices. Asa special case, a uniform distribution over the set of d × d positive definite correlation matrices obtains.
Byproducts are: (a) For a uniform distribution over the space of d × d correlation matrices, the marginal
distribution of each correlation is Beta(d/2, d/2) on (−1, 1). (b) An identity is obtained for the determinant ofa correlation matrix R via partial correlations¡ in¢ a vine. (c) A formula is obtained for the volume of the set ofd × d positive definite correlation matrices in d2 -dimensional space.
Outline1. Statement of key results on generating random correlation matrices
2. Parametrization with partial correlations
3. Regular vines, D-vines, C-vines
4. Random correlation matrices based on partial correlations
5. Random correlation matrices based on onion method
Key results: R is a d × d correlation matrix
1. Several simple ways to generate a random R that is uniform over set of positive definite d × d correlation
matrices; more generally with density ∝ [det(R)]α−1 .
2. Identity: (1 − ρ212 )(1 − ρ223 )(1 − ρ213;2 ) for d = 3det(R) =d−1
Yi=1

(1 − ρ2i,i+1 ) ×

d−1Y d−kY(1 − ρ2j,j+k;j+1...j+k−1 ).

k=2 j=1[partial correlations in the double product]
3. Volume of d × d correlation matrices in 2d(d−1)/2 dimensional space is:Q(d−1)/2Γ(2m) π (d2 −1)/4 (d−1)m=1
if d is odd;
2 /4 d−1 d+1 ,2Γ( 2 )2 −4d)/4 d d Q(d−2)/2(3d
Γ(2m)
Γ (2)
 π d(d−2)/4 2
m=1, if d is even.
Γd−1 (d)Partial correlations
correlations ρi,i+1 for i = 1, . . . , d − 1
the partial correlations ρij;i+1,...j−1 for j − i ≥ 2
For example, ρ12 , ρ23 , ρ34 , ρ13;2 , ρ24;3 , ρ14;23 for d = 4.
For multivariate normal random variables, partial correlations are conditional correlations. In general, they
are functions of a correlation
matrix.
¡ ¢Advantage: the d2 parameters can independently take values in the interval (−1, 1); different vines can beused.
Disadvantage: the reparametrization is non-unique, and depends on the order of indexing, e.g.: ρ12 , ρ13 , ρ23;1or ρ12 , ρ23 , ρ13;21

d = 5: 3 different vines + matrix for onion method are shown.
d-vine which is like Markovian dependence c-vine
regular partial correlation vine
onion method#distinct vines increases quickly as d increases,partial correlations independently in (−1, 1).
1
ρ
 12

 ρ31

ρ41
ρ51

ρ121
ρ32
ρ42
ρ52



1
 ρ12

 ρ31

ρ41
ρ51

ρ13;2
ρ23
1
ρ43
ρ53
ρ12
1
ρ32
ρ42
ρ52

ρ14;23
ρ24;3
ρ34
1
ρ54

ρ13;2
ρ23
1
ρ43
ρ53


ρ15;234
ρ25;34 ρ35;4 

ρ45
1



1
ρ
 12

 ρ31

ρ41
ρ51


ρ15;24
ρ25;4 
ρ35;124 

ρ45
1

ρ14;2
ρ24
ρ34;121
ρ54



ρ121
ρ32
ρ42
ρ52

1
 ρ12

 ρ31

ρ41
ρ51

ρ13
ρ23;11
ρ43
ρ53
ρ121
ρ32
ρ42
ρ52


ρ15
ρ25;1 
ρ35;12 ρ45;123
1

ρ14
ρ24;1
ρ34;121
ρ54

ρ13
ρ23
1
ρ43
ρ53


ρ15
ρ25 

ρ35 

ρ45
1

ρ14
ρ24
ρ34
1
ρ54

D-vine for d = 5: specify densities separately for each partial correlation; then density of the correlation
matrix R isfR (r) =×××fρ12 (r12 )fρ23 (r23 )fρ34 (r34 )fρ45 (r45 )fρ13;2 (r13;2 )fρ24;3 (r24;3 )fρ35;4 (r35;4 )fρ14;23 (r14;23 )fρ25;34 (r25;34 )fρ15;234 (r15;234 )Jacobian

Under appropriate choices, this leads to a uniform density.
References for vines
D. Kurowicka, R. Cooke, Uncertainty Analysis with High Dimensional Dependence Modelling, Wiley, 2006.
Vines are a graphical method for modeling multivariate dependencies
3{23|1}

2{34|12}

{12}

{12}

1

2{23}

3{13}

1

{14}

{13|2}

4{24|3}{14|23}{24|1}

2{34}

4

3
{13|2}

1

{12}

{23}
{24}

2

4{45}

5{34|12}{25|4}{14|2}{15|24}{35|124}

Random correlation matrices [Joe, 2006, JMVA]
Partial correlation approach: distributions separately for ρi,i+1 for i = 1, . . . , d − 1, and ρij;i+1,...j−1 forj − i ≥ 2.
This leads to a joint distribution of {ρij : 1 ≤ i < j ≤ d}.
More generally, choose distributions separately for partial correlation in a regular vine.
Previous papers : combination of random eigenvalues and random (orthogonal) matrices, or a random T
matrix leading to T 0 T , but cannot do distribution theory?
Motivations:• generation of ellipsoidal clusters to study clustering methods• efficiency of estimating eq’n methods for multiv probit model• assess probability that a correlation matrix is a rank correlation matrix.
• prior for correlation matrix in Bayesian statistics?
Main result
With appropriate choices of beta distributions for ρi,i+1 for i = 1, . . . , d − 1, and ρij;i+1,...j−1 for j − i ≥ 2,can get a joint density for {ρij : 1 ≤ i < j ≤ d} that is ∝ |R|α−1 , where α > 0.
More generally, same for partial correlations on a regular vine.
In this case, the joint density is invariant to the order of indexing of variables for the partial correlations(and to the vine), and each ρij marginally has a Beta(α − 1 + d/2, α − 1 + d/2) distribution on (−1, 1).
Special case with α = 1: uniform density over the set of positive definite correlation matrices.
Anderson (1958, p. 80): partial corr ρj,j+k;j+1,...,j+k−1 (2 ≤ k ≤ d − 1) ishwhereρj,j+k − r01 (j, k)(R2 (j, k))−1 r3 (j, k)
i1/2 hi1/2
1 − r01 (j, k)(R2 (j, k))−1 r1 (j, k)
1 − r03 (j, k)(R2 (j, k))−1 r3 (j, k)1r01 (j, k) ρj,j+k
R[j : j + k] =  r1 (j, k) R2 (j, k) r3 (j, k)  ,ρj+k,j r03 (j, k)
1with r01 (j, k) = (ρj,j+1 , . . . , ρj,j+k−1 ), r03 (j, k) = (ρj+k,j+1 , . . . , ρj+k,j+k−1 ), and R2 (j, k) consisting of the middle
k − 1 rows / columns. Henceρj,j+k = r01 (j, k)(R2 (j, k))−1 r3 (j, k) + ρj,j+k;j+1,...,j+k−1 Djk ,3wherehihi2Djk= 1 − r01 (j, k)(R2 (j, k))−1 r1 (j, k) 1 − r03 (j, k)(R2 (j, k))−1 r3 (j, k) .

d = 3 [ρ for rv, r for dummy variable of function]
Consider (ρ12 , ρ23 , ρ13·2 ) → (ρ12 , ρ23 , ρ13 ). Sinceρ13 − ρ12 ρ23

ρ13·2 = p

(1 − ρ212 )(1 − ρ223 ),the Jacobian of the transformation is:22)(1 − r23)]−1/2 .
∂r13·2 /∂r13 = [(1 − r12Note that2
1 − r13·2=222det(R)1 − r12− r23− r13+ 2r12 r23 r13=222 )(1 − r 2 )(1 − r12 )(1 − r23 )(1 − r1223222det(R(r12 , r23 , r23 )) = 1 − r12− r23− r13+ 2r12 r23 r13222))(1 − r13·2)(1 − r23= (1 − r12(D1)

1−1(1 − u2 )α−12 [B(α, α)]Beta(α, α) density on (−1, 1):
or U = 2V − 1, V ∼ Beta(α, α) on (0, 1).
ρ12 , ρ23 ∼ Beta(α1 , α1 ),ρ13·2 ∼ Beta(α2 , α2 ) independently.

222)]α1 −3/2)(1 − r23)α2 −1 [(1 − r12fρ12 ,ρ23 ,ρ13 (r12 , r23 , r13 ) ∝ (1 − r13·2If α1 = α2 + 12 , then the density is∝ [det(R)]α2 −1(D2)This is symmetric in r12 , r13 , r23 , which means the same joint density obtains from ρ12 , ρ13 ∼ Beta(α1 , α1 ), ρ23·1 ∼Beta(α2 , α2 ) etc.
α2 = 1, α1 = 3/2 ⇒ Uniform over set of correlation matrices.
d > 3 : computer proof before math proof
Based on the results for the d = 3 case, conjectured an extension of (D1) as an identity for det(R), andconjectured the Beta distributions on ρij;i+1,...,j−1 needed to get joint density of (ρij ) to be ∝ [det(R)]α−1 .
The identity was verified numerically and the univariate margins of the random correlation matrix were also
checked numerically before proving the results for general d.
Results for d > 3.
Thm 1:
det(R) =d−1Y

(1 − ρ2i,i+1 ) ×i=1d−1
Y d−kY(1 − ρ2j,j+k;j+1...j+k−1 ).

k=2 j=1

Thm 4: The determinant |Jd | of the Jacobian for the transform of (ρ12 , ρ23 , ρ13 , ρ34 , ρ24 , ρ14 , ρ45 , . . . , ρ1d ) to(ρ12 , ρ23 , ρ13;2 , ρ34 , ρ24;3 , ρ14;23 , ρ45 , . . . , ρ1d;2...d−1 ) is:
hd−1
Y
i=1

(1 − ρ2i,i+1 )d−2 ×d−2Y d−kY

(1 − ρ2i,i+k;i+1...i+k−1 )d−1−kk=2 i=14

i−1/2.

Thm 5: If αk = αd−1 + 21 (d − 1 − k), k = 1, . . . , d − 1, and ρi,i+k;i+1...i+k−1 is Beta(αk , αk ) on (−1, 1) for1 ≤ i < i + k ≤ d, then the joint density becomes£ ©
ª¤αd−1 −1
c−1det (rij )1≤i,j≤d,dwhere the normalizing constant cd is
Pd−1

2

k=1

(2αd−1 −2+d−k)(d−k)×d−1Y£ ¡¢¤d−k
B αd−1 + 12 (d − 1 − k), αd−1 + 12 (d − 1 − k).

k=1

If αd−1 = 1 and αk = 12 (d + 1 − k) for k = 1, . . . , d − 2, leading to uniform joint density for {ρij , i < j}, then thenormalizing constant is
cd = 2

Pd−1
k=1(d−k)2

×d−1
Y£ ¡1
¢¤d−k
B 2 (d − k + 1), 21 (d − k + 1)

k=1

=2

Pd−1 2
k=1

k×d−1
Y£ ¡1
¢¤kB 2 (k + 1), 12 (k + 1) ,k=1

and the recursion is£¤d−12
cd = cd−1 × 2(d−1) × B( 12 d, 12 d).

¡d¢ Byproduct: normalizing constant is the volume of the set of d-dim positive definite correlation matrices in
2 -dim spaced cd
cd /2d(d−1)/22 2
3 4.934802 = π 2 · (1/2)
0.617
4 11.69731 = π 2 · (32/27)
0.183
5 22.53256 = π 6 · (3/128)
0.022
6 31.11388 = π 6 · (8192/253125) 0.00095
7 27.85823 = π 12 · (n7 /d7 )0.000013
8 14.87740 = π 12 · (n8 /d8 )5.5 × 10−820
9 4.411544 = π · (n9 /d9 )6.4 × 10−1120
10 0.682269 = π · (n10 /d10 )1.9 × 10−14
cd /2d(d−1)/2 decreases quickly as d increases.
Random correlation matrix based on vines [Lewandowski, Kurowicka and Joe, 2007]
Extensions of Theorems 1,4,5 in Joe (2006, JMVA).
V = (e) is a partial correlation vine;
for example, V = (12, 13, 14, 23; 1, 24; 1, 34; 12) withdepth=(order−1)=(#conditioned variables): (0, 0, 0, 1, 1, 2).
For a vine based on d variables, there are d − 1 edges with depth ke = 0, d − 2 edges with depth ke = 1, . . ., and1 edge with depth ke = d − 2.
Let R be a (random) correlation matrix based on distributions for {ρe : e ∈ V }.
T1: det(R) =Q2e∈V (1 − ρe ) [Kurowicka and Cooke, 2006, LAA] (product over edges of the vine)5

T4: Jacobian of transform from correlations {ρij } to partial correlations {ρe } is(Y)−1/2(1 − ρ2e )d−2−kee∈V

T5. Let ρe ∼ Beta(αke +1 , αke +1 ) on (−1, 1) independently, where αke +1 = αd−1 + 21 (d − 2 − ke ). Get
uniform over space of correlation matrices if αd−1 = 1, and then marginal distribution for each correlation isBeta(d/2, d/2).
The C-vine algorithm for generating correlation matrices with density ∝ [det(r)]η−1 , η > 1.
1. Initialization β = η + (d − 1)/2.
2. Loop for k = 1, . . . , d − 1.
a) β ← β − 12 ;b) Loop for i = k + 1, . . . , d;i) generate partial corr pk,i;1,...,k−1 ∼ Beta(β, β) on (−1, 1);ii) use recursive formula on pk,i;1,...,k−1 to get rk,i = ri,k [no matrix inversions];
3. Return r, a d × d correlation matrix.

Onion method for random correlation matrices[Ghosh & Henderson, 2003; extended LKJ 2007]
Result 1. Consider the spherical density c(1−w0 w)β−1 for w ∈ <m , w0 w ≤ 1, where c = Γ(β+m/2)π −m/2 /Γ(m/2).
If W has this density, then W = V U where V ∼ Beta(m/2, β) and U is uniform on the surface of the mdimensional hypersphere. If Q = AW, where A is an m × m non-singular matrix, then the density of Q isc[det(AA0 )]−1/2 (1 − q0 [AA0 ]−1 q)β−1 ,q 3 q0 [A0 A]−1 q ≤ 1.

µ¶rm qResult 2. Partition rm+1 =where rm is an m × m correlation matrix and q is a m-vector ofq0 1
correlations such that rm+1 is an (m + 1) × (m + 1) correlation matrix. Then¡¢det(rm+1 ) = det(rm ) · 1 − q0 r−1
m q .
We use upper case letter of rm , q, rm+1 to denote random vectors and matrices. Let β, βm > 0. If Rm hasdensity ∝ [det(rm )]βm −1 andβ−1
Q given Rm = rm has density ∝ [det(rm )]−1/2 (1 − q0 r−1,m q)βm −3/2
0 −1 β−1then the density of Rm+1 is ∝ [det(rm )](1 − q rm q).
If βm = β + 21 , then the density of Rm+1 is ∝ [det(rm+1 )]β−1 .
Algorithm for the extended onion method to get random correlation matrices in dimension d with density∝ [det(r)]η−1 , η > 0.
µ¶1 r12
1. Initialization. β = η + (d − 2)/2, r12 ← 2u − 1, where u ∼ Beta(β, β), r ←r12 1
2. Loop for m = 2, . . . , d − 1.
6

(a) β ← β − 12 ;(b) generate y ∼ Beta(m/2, β), generate u = (u1 , . . . , um )0 uniform on the surface of m-dimensional
hypersphere;(c) w ← y 1/2 u, obtain Cholesky decomp AA0 = r, set q ← Aw;µ¶
r q(d) r ←.
q0 1
3. Return r, a d × d correlation matrix.

Why is this called the onion method ???
By the symmetry, each ρjk = Rjk in the correlation matrix has a marginal Beta(η + [d − 2]/2, η + [d − 2]/2)
density on (−1, 1). For the special case of η = 1 leading to uniform over the space of correlation matrices, themarginal density of each Rjk is Beta(d/2, d/2) on (−1, 1).
Computational time:
In C, onion method with incremental Cholesky is fastest, and C-vine is faster than onion method without incremental Cholesky.
In Matlab, C-vine (which avoids matrix inversions) is fastest.
D-vine is slower because of matrix inversions.
Other comments:
If partial correlations in a particular vine are needed, then use the appropriate vine.
Choosing the vine and the densities for the partial corr in this vine, one could get random correlation matricesthat have larger correlations at a few particular pairs.
References
H. Joe, Generating random correlation matrices based on partial correlations, Journal of Multivariate Analysis
97 (2006) 2177–2189.
S. Ghosh, S.G. Henderson, Behavior of the NORTA method for correlated random vector generation as thedimension increases. ACM Transactions on Modeling and Computer Simulation (TOMACS) 13 (2003), 276–294.
T. Bedford, R. Cooke, Vines - a new graphical model for dependent random variables, Annals of Statistics 30(2002) 1031–1068.
D. Kurowicka, R. Cooke, Uncertainty Analysis with High Dimensional Dependence Modelling, Wiley, 2006.
D. Kurowicka, R. M. Cooke, Completion problem with partial correlation vines, Linear Algebra and Its
Applications 418 (2006) 188–200.
D. Lewandowski, D. Kurowicka, H. Joe, Generating random correlation matrices based on vines and extendedOnion method, (2007), submitted.

7