[Paper Review] An efficient algorithm for sampling from $\sin^k(x)$ for generating random correlation matrices

Enes Makalic, Daniel F. Schmidt|arXiv (Cornell University)|Sep 14, 2018

Cellular Automata and Applications1 references4 citations

TL;DR

This paper presents a highly efficient rejection sampling algorithm for generating random variates from the distribution proportional to $\sin^k(x)$ on $(0, \pi)$, using a beta distribution as an envelope. The method achieves a maximum expected acceptance rate of $\pi/(2\sqrt{2}) \approx 1.11$ iterations per sample across all $k \geq 1$, enabling fast generation of high-dimensional random correlation matrices via hyperspherical parameterization.

ABSTRACT

In this note, we develop a novel algorithm for generating random numbers from a distribution with a probability density function proportional to $\sin^k(x)$, $x \in (0,π)$ and $k \geq 1$. Our algorithm is highly efficient and is based on rejection sampling where the envelope distribution is an appropriately chosen beta distribution. An example application illustrating how the new algorithm can be used to generate random correlation matrices is discussed.

Motivation & Objective

• To develop a computationally efficient method for sampling from the probability density proportional to $\sin^k(x)$ on $(0, \pi)$, which arises in random correlation matrix generation.
• To overcome the inefficiency of inverse transform sampling, which requires $O(k)$ operations and becomes impractical for large $k$.
• To provide a numerically stable alternative to numerical integration, which suffers from precision issues for large $k$.
• To enable scalable generation of high-dimensional correlation matrices by efficiently sampling the required angles $\theta_{ij}$.
• To implement and benchmark the algorithm in widely used platforms like R and MATLAB for practical use.

Proposed method

• The algorithm uses rejection sampling with a symmetric $\text{Beta}(k+1, k+1)$ distribution as the envelope density.
• The acceptance criterion is based on the ratio $f_X(x)/g_Y(x)$, with acceptance if $\frac{1}{k}\log U \leq \log\left(\frac{\pi^2 \sin X}{4X(\pi - X)}\right)$.
• The upper bound $M_k = \frac{f_X(\pi/2)}{g_Y(\pi/2)}$ is derived analytically, ensuring the envelope condition $f_X(x) \leq M_k g_Y(x)$.
• The normalization constants for $f_X(x)$ and $g_Y(x)$ are expressed using the Gamma and Beta functions, enabling exact computation.
• The method is validated by applying it to generate $p \times p$ correlation matrices using the hyperspherical parameterization of Pourahmadi and Wang.
• The algorithm is implemented in both R (via the randcorr package) and MATLAB (File ID: 68810) for public use.

Experimental results

Research questions

• RQ1 Can a rejection sampling approach be designed to efficiently sample from $\sin^k(x)$ on $(0, \pi)$ for all $k \geq 1$?
• RQ2 What envelope distribution yields a tight bound on the acceptance ratio while allowing fast sampling?
• RQ3 What is the theoretical maximum expected number of iterations per sample across all $k \geq 1$?
• RQ4 How does the proposed method compare in performance to inverse transform sampling for high-dimensional correlation matrices?
• RQ5 Can the algorithm be efficiently scaled to generate large $p \times p$ correlation matrices with $p$ up to 1000?

Key findings

• The maximum expected number of iterations per sample is $\pi/(2\sqrt{2}) \approx 1.11$, achieved in the limit as $k \to \infty$.
• For $k=1$, the expected number of iterations is $\pi/3 \approx 1.047$, increasing monotonically with $k$.
• The algorithm generates a $100 \times 100$ random correlation matrix in approximately 0.02 seconds on a commodity laptop.
• A $1000 \times 1000$ correlation matrix is generated in about 0.5 seconds, demonstrating linear scaling with matrix size.
• The proposed method outperforms inverse transform sampling, which required 0.9s for $100 \times 100$ and 42s for $600 \times 600$ matrices on comparable hardware.
• The rejection sampler is numerically stable and avoids the precision issues of numerical integration for large $k$.