A direct method to generate correlation matrices with specified eigenvalues

In a previous article, I implemented an algorithm due to Niels Waller (TAS, 2020) that uses the method of alternating projections (MAP) to generate random correlation matrices that have a specified set of eigenvalues. The algorithm is iterative, and the MAP method is not guaranteed to converge, although Waller claims that the method worked well in a simulation study. MAP and other indirect methods tend to be computationally expensive, so I was happy to discover a paper by Marsaglia and Olkin (1984) that uses a direct method to generate random correlation matrices. The Marsaglia and Olkin algorithm is straightforward to implement in a high-level language like SAS IML. This article presents an IML function that implements the Marsaglia and Olkin direct method for generating a random correlation matrix.

Correlation matrices that have a common spectrum

A matrix is positive semi-definite if all eigenvalues are greater than or equal to zero. All correlation matrices are symmetric and positive semi-definite.

The spectrum of the matrix is the set of eigenvalues: Λ = {λ₁, λ₂, ..., λ_d}. For a positive definite matrix, λ_i > 0 for all i. For a semi-definite matrix, λ_i ≥ 0. For any square matrix, the sum of the eigenvalues equals the trace of the matrix. Consequently, for an d x d correlation matrix, Σ_i λ_i = d.

As shown in the previous article, there are many different correlation matrices that have the same spectrum.

The Marsaglia and Olkin method for random correlation matrices

The Marsaglia-Olkin (1984) paper is behind a paywall, so I have not read it. Fortunately, Gentle (2003, p. 200) presents pseudo-code for their algorithm. Gentle's pseudo-code contains two typos (Steps 5 and 6), which I have corrected in the following description.

The input to the Marsaglia-Olkin routine is Lambda, which is the spectrum of a d x d SPD correlation matrix. Lambda is a d x 1 vector {λ₁, λ₁, ..., λ_d}, where λ_i > 0 and Σ_i λ_i = d. The algorithm outputs R, which is a random correlation matrix such that the eigenvalues of R are the specified spectrum. The algorithm is as follows:

Set E = I_d and k=1.
Generate a d x 1 vector, w, of i.i.d. standard normal random variates.
Form x = E*w and compute the scalar a = (1-Lambda)` * x##2
Generate a d x 1 vector, z, of i.i.d. standard normal random variates.
Form y = E*z and compute the scalars:
- b = (1-Lambda)` * (x#y)
- c = (1-Lambda)` * y##2
- e² = b##2 - a*c
If e² < 0, goto Step 2.
Choose a random sign, s ∈ {-1, 1}. Set r = (b+s*e)/a * x - y.
Choose another random sign, s ∈ {-1, 1}. Set p_k = s*r/norm(r). (Fixes typo in Gentle.) This vector will be a row in the P matrix. It is a unit vector in the hypercone defined by the equation p_k`(I - *Lambda;)p_k = 0.
Perform a rank-one update: E → E - p_k`*p_k (fixes typo), and increment k=k+1. At each step, E is the projection matrix onto the subspace orthogonal to the rows of P that have been constructed so far.
If k < d, goto Step 1.
Generate a d x 1 vector, w, of i.i.d. standard normal random variates.
Form x = E*w and set p_d = x / norm(x).
Construct the matrix P by using the vectors p_k as its rows.
Return P*diag(Lambda)*P` as the random correlation matrix.

I always say that the main benefit of a matrix language such as MATLAB, R, and SAS IML is the ease of converting a high-level pseudo-code description of an algorithm into a ready-to-run program. The Appendix to this article contains my implementation of the algorithm in a SAS IML moduled called CorrWithEigen_MO. The next section shows how to use the function to generate a random correlation matrix.

Generate random correlation matrices in SAS

First, run the program in the Appendix, which stores the CorrWithEigen_MO function. You can then load the function and use it in other programs. The following call to PROC IML generates two random 6 x 6 correlation matrices that have the eigenvalues Lambda = {1.8, 1.5, 1, 1, 0.5, 0.2}.

proc iml;
load module=(CorrWithEigen_MO);
call randseed(123);
reset fuzz;   /* print tiny numbers as 0 */
 
/* Specify the spectrum for a 6x6 corr matrix */
Lambda = {1.8, 1.5, 1, 1, 0.5, 0.2};
/* get random 6x6 correlation matrix with this spectrum */
Corr1 = CorrWithEigen_MO( Lambda );
/* call again, make sure we get a different answer */
Corr2 = CorrWithEigen_MO( Lambda );
print Corr1[F=best6.], Corr2[F=best6.];
 
/* verify the spectrum of the random matrices */
v1 = eigval(Corr1);
v2 = eigval(Corr2);
print Lambda v1 v2;

The output shows two different random 6x6 correlation matrices. Each matrix has the same set of eigenvalues, as demonstrated by the calls to the EIGVAL function.

Summary

This article provides an IML implementation of the Marsaglia and Olkin (1984) algorithm for the direct generation of random correlation matrices that have a specified spectrum. You can use the algorithm to generate random correlation or covariance matrices in SAS. You can use this function in simulation studies where you need a positive definite matrix.

Appendix: The Marsaglia and Olkin Algorithm in IML

This section defines an IML function (and some helper functions) that implement the Marsaglia and Olkin (1984) algorithm for the direct generation of random correlation matrices that have a specified spectrum.

/* From Gentle (2003, p. 200) Random Number Generation and Monte Carlo Methods
   Algorithm 5.9 
   Marsaglia-Olkin (1984) Method for Random Correlation Matrices with Given Eigenvalues
   Ref: Marsaglia, G. and Olkin, I. (1984) "Generating Correlation Matrices", 
        SIAM J. on Sci. and Stat. Comp., 5(2), pp 470-475. doi 10.1137/0905034.
 
   INPUT:
   Lambda = the spectrum of a dxd SPD correlation matrix. Lambda is a vector
            (lambda_1, ..., lambda_d), where lambda_i > 0 and \sum lambda_i = d.
   OUTPUT:
   R = dxd correlation matrix, R, such that the eigenvalues of R are the specified spectrum:
        Lambda(R) = (lambda_1, ..., lambda_d)
   ALGORITHM:
   0. Set E = I_d and k=1.
   1. Generate d-vector, w, of i.i.d. std normal deviate.
      Form x = E*w and compute the scalar a = (1-Lambda)` * x##2
   2. Generate a d-vector, z, of i.i.d. std normal variates. 
      Form y = E*z and compute the scalars 
      b = (1-Lambda)` * (x#y)
      c = (1-Lambda)` * y##2
      e^2 = b##2 - a*c
   3. If e^2 < 0, goto Step 2.
   4. Choose a random sign, s \in {-1, 1}. Set r = (b+s*e)/a * x - y.
   5. Choose another random sign, s \in {-1, 1}. Set p_k = s*r/norm(r). (fixes typo in Gentle)
   6. Perform a rank-one update: E -> E - p_k`*p_k (fixes typo), and increment k=k+1.
   7. If k < d, goto Step 1.
   8. Generate a d-vector, w, of i.i.d. std normal deviate.
      Form x = E*w and set p_d = x / norm(x).
   9. Construct the matrix P by using the vectors p_k as its rows.
      Return P*diag(Lambda)*P` as the random correlation matrix.
*/
 
proc iml;
/* return k random signs with values +1 or -1.
   Let B be a binary random variable chosen uniformly at random from {0,1}.
   Then C = 2*(B - 1/2) is random uniform in {-1,1}.
*/
start RandSign(k);
   return( 2*(randfun(k, "Bernoulli", 0.5) - 0.5) );
finish RandSign;
 
/* for column vector v, standardize so sum(v)=nrow(v) */
start StdizeEigenval(v);
    L = colvec(v);
    return ( L * nrow(L)/sum(L) );  /* make sum(v)=dimension */
finish StdizeEigenval;
 
/* The  Marsaglia-Olkin (1984) method for creating a random correlation 
   matrix with a specified set of eigenvalues */
start CorrWithEigen_MO(targetLambda, maxIters=100);
   lambda = StdizeEigenval(targetLambda);
   d = nrow(Lambda);
   if any(Lambda < 0) then do; 
      print "ERROR: The spectrum must contain only nonnegative values"; 
      return( J(d,d,.) ); 
   end;
   /* If the spectrum is {1,1,...,1}, the only solution is I(d).
      The do-while logic in the M-O method will enter an infinite loop
      in this case, so detect this edge case and return the identity matrix */
   if norm(1 - Lambda, "LInf") < 1e-8 then return (I(d));
 
   /* use column vectors for Lambda, w, and z */
   w = j(d, 1, .);
   z = j(d, 1, .);
   P = j(d, d, .);
   /* 0. Set E = I_d and k=1. */
   E = I(d);
 
   /* compute the k_th row of P */
   do k = 1 to d-1;
      /* 1. Generate dx1 vector, w, of i.i.d. std normal variates.
            Form x = E*w and compute the scalar a = (1-Lambda)` * x##2. */
      call randgen(w, "Normal");
      x = E*w;
      a = sum((1-Lambda) # x##2);
      /* 2. Generate a dx1 vector, z, of i.i.d. std normal variates. 
            Form y = E*z and compute the scalars 
            b = (1-Lambda)` * (x#y)
            c = (1-Lambda)` * y##2
            e2 = b##2 - a*c
         3. If e2 < 0, goto Step 2.     */
      e2 = -1;
      do cnt = 1 to maxIters until(e2 > 0);
         call randgen(z, "Normal");
         y = E*z;
         b = sum((1-Lambda) # x # y);
         c = sum((1-Lambda) # y##2);
         e2 = b##2 - a*c;
      end;
      if e2 <= 0 then do;
         print "ERROR:  Marsaglia-Olkin method failed after 100 iterations.";
         return (J(d, d, .)); 
      end;
      /* 4. Choose two random signs, s \in {-1, 1} */
      s = RandSign(2);
      r = (b+s[1]*sqrt(e2))/a * x - y;   /* Set r = (b+s*e)/a * x - y (Note Pr(a=0)=0) */
      p_k = s[2] * r/norm(r);            /* 5. Set p_k = s/norm(r) * w */
 
      /* 6. Perform a rank-one update: E -> E - r*r` */
      E = E - p_k*p_k`;
      P[,k] = p_k;
   end;   /* 7. If k < d, goto Step 1. */
 
   /* 8. Generate a dx1 vector, w, of i.i.d. std normal variates.
         Form x = E*w and set p_d = x / norm(x). */
   call randgen(w, "Normal");
   x = E*w;
   P[,d] = x / norm(x);
   /* 9. Construct the matrix P by using the vectors p_k as its rows.
         Return P*diag(Lambda)*P` as the random correlation matrix.
         (Note: Instead, I will construct P by using p_k as columns, then return P`*Lambda*P) */
   G = P` # sqrt(Lambda`);
   Corr = G*G`;
   return( Corr );
finish CorrWithEigen_MO;
store module=(RandSign StdizeEigenval CorrWithEigen_MO);
QUIT;

Blogs