¿Cómo generar eficazmente matrices de correlación aleatorias semidefinidas positivas?

Question

Me gustaría ser capaz de generar eficientemente matrices de correlación positivas-semidefinidas (PSD). Mi método se ralentiza drásticamente a medida que aumenta el tamaño de las matrices a generar.

1. ¿Podría sugerir alguna solución eficaz? Si conoces algún ejemplo en Matlab, te lo agradecería mucho.
2. A la hora de generar una matriz de correlación PSD, ¿cómo elegiría los parámetros para describir las matrices que se van a generar? ¿Una correlación media, la desviación estándar de las correlaciones, los valores propios?

Answer 1

Un papel Generación de matrices de correlación aleatorias basadas en vides y extendidas método de la cebolla de Lewandowski, Kurowicka y Joe (LKJ), 2009, proporciona un tratamiento y exposición unificados de los dos métodos eficientes de generación de matrices de correlación aleatoria. Ambos métodos permiten generar matrices de una distribución uniforme en un cierto sentido preciso que se define a continuación, son sencillas de aplicar, rápidas y tienen la ventaja añadida de tener nombres divertidos.

Una matriz simétrica real de $d \times d$ tamaño con unos en la diagonal tiene $d(d-1)/2$ elementos únicos fuera de diagonal y por lo tanto se puede parametrizar como un punto en $\mathbb R^{d(d-1)/2}$ . Cada punto de este espacio corresponde a una matriz simétrica, pero no todas son definidas positivamente (como tienen que ser las matrices de correlación). Por tanto, las matrices de correlación forman un subconjunto de $\mathbb R^{d(d-1)/2}$ (en realidad un subconjunto convexo conectado), y ambos métodos pueden generar puntos a partir de una distribución uniforme sobre este subconjunto.

Proporcionaré mi propia implementación en MATLAB de cada método y los ilustraré con $d=100$ .

---

Método de la cebolla

El método de la cebolla proviene de otro trabajo (ref #3 en LKJ) y debe su nombre al hecho de que las matrices de correlación se generan a partir de $1\times 1$ matriz y haciéndola crecer columna a columna y fila a fila. La distribución resultante es uniforme. Realmente no entiendo las matemáticas detrás del método (y prefiero el segundo método de todos modos), pero aquí está el resultado:

Aquí y debajo del título de cada subtrama se muestran los valores propios más pequeños y más grandes, y el determinante (producto de todos los valores propios). Aquí está el código:

%// ONION METHOD to generate random correlation matrices distributed randomly
function S = onion(d)
    S = 1;
    for k = 2:d
        y = betarnd((k-1)/2, (d-k)/2); %// sampling from beta distribution
        r = sqrt(y);
        theta = randn(k-1,1);
        theta = theta/norm(theta);
        w = r*theta;
        [U,E] = eig(S);
        R = U*E.^(1/2)*U';             %// R is a square root of S
        q = R*w;
        S = [S q; q' 1];               %// increasing the matrix size
    end
end

Método de la cebolla extendida

LKJ modifica ligeramente este método para poder muestrear las matrices de correlación $\mathbf C$ de una distribución proporcional a $[\mathrm{det}\:\mathbf C]^{\eta-1}$ . Cuanto más grande sea el $\eta$ cuanto mayor sea el determinante, lo que significa que las matrices de correlación generadas se acercarán cada vez más a la matriz de identidad. El valor $\eta=1$ corresponde a una distribución uniforme. En la figura siguiente las matrices se generan con $\eta={1, 10, 100, 1000, 10\:000, 100\:000}$ .

Por alguna razón para obtener el determinante del mismo orden de magnitud que en el método de la cebolla de vainilla, necesito poner $\eta=0$ y no $\eta=1$ (como afirma LKJ). No sé dónde está el error.

%// EXTENDED ONION METHOD to generate random correlation matrices
%// distributed ~ det(S)^eta [or maybe det(S)^(eta-1), not sure]
function S = extendedOnion(d, eta)
    beta = eta + (d-2)/2;
    u = betarnd(beta, beta);
    r12 = 2*u - 1;
    S = [1 r12; r12 1];  

    for k = 3:d
        beta = beta - 1/2;
        y = betarnd((k-1)/2, beta);
        r = sqrt(y);
        theta = randn(k-1,1);
        theta = theta/norm(theta);
        w = r*theta;
        [U,E] = eig(S);
        R = U*E.^(1/2)*U';
        q = R*w;
        S = [S q; q' 1];
    end
end

Método de la vid

El método Vine fue sugerido originalmente por Joe (J en LKJ) y mejorado por LKJ. A mí me gusta más, porque es conceptualmente más fácil y también más fácil de modificar. La idea es generar $d(d-1)/2$ correlaciones parciales (son independientes y pueden tener cualquier valor de $[-1, 1]$ sin ninguna restricción) y luego convertirlas en correlaciones brutas mediante una fórmula recursiva. Es conveniente organizar el cómputo en un determinado orden, y este gráfico se conoce como "vid". Es importante destacar que si las correlaciones parciales se muestrean a partir de distribuciones beta particulares (diferentes para las distintas celdas de la matriz), entonces la matriz resultante se distribuirá uniformemente. También en este caso, LKJ introduce un parámetro adicional $\eta$ para tomar una muestra de una distribución proporcional a $[\mathrm{det}\:\mathbf C]^{\eta-1}$ . El resultado es idéntico al de la cebolla extendida:

%// VINE METHOD to generate random correlation matrices
%// distributed ~ det(S)^eta [or maybe det(S)^(eta-1), not sure]
function S = vine(d, eta)
    beta = eta + (d-1)/2;   
    P = zeros(d);           %// storing partial correlations
    S = eye(d);

    for k = 1:d-1
        beta = beta - 1/2;
        for i = k+1:d
            P(k,i) = betarnd(beta,beta); %// sampling from beta
            P(k,i) = (P(k,i)-0.5)*2;     %// linearly shifting to [-1, 1]
            p = P(k,i);
            for l = (k-1):-1:1 %// converting partial correlation to raw correlation
                p = p * sqrt((1-P(l,i)^2)*(1-P(l,k)^2)) + P(l,i)*P(l,k);
            end
            S(k,i) = p;
            S(i,k) = p;
        end
    end
end

Método de la vid con muestreo manual de correlaciones parciales

Como se puede ver arriba, la distribución uniforme da lugar a matrices de correlación casi diagonales. Pero se puede modificar fácilmente el método de la vid para tener correlaciones más fuertes (esto no se describe en el documento de LKJ, pero es sencillo): para ello hay que muestrear las correlaciones parciales de una distribución concentrada alrededor de $\pm 1$ . A continuación los muestro a partir de la distribución beta (reescalada de $[0,1]$ a $[-1, 1]$ ) con $\alpha=\beta={50, 20, 10, 5, 2, 1}$ . Cuanto más pequeños son los parámetros de la distribución beta, más se concentra cerca de los bordes.

Nótese que en este caso no se garantiza que la distribución sea invariante de la permutación, por lo que permuto adicionalmente las filas y columnas después de la generación.

%// VINE METHOD to generate random correlation matrices
%// with all partial correlations distributed ~ beta(betaparam,betaparam)
%// rescaled to [-1, 1]
function S = vineBeta(d, betaparam)
    P = zeros(d);           %// storing partial correlations
    S = eye(d);

    for k = 1:d-1
        for i = k+1:d
            P(k,i) = betarnd(betaparam,betaparam); %// sampling from beta
            P(k,i) = (P(k,i)-0.5)*2;     %// linearly shifting to [-1, 1]
            p = P(k,i);
            for l = (k-1):-1:1 %// converting partial correlation to raw correlation
                p = p * sqrt((1-P(l,i)^2)*(1-P(l,k)^2)) + P(l,i)*P(l,k);
            end
            S(k,i) = p;
            S(i,k) = p;
        end
    end

    %// permuting the variables to make the distribution permutation-invariant
    permutation = randperm(d);
    S = S(permutation, permutation);
end

Así es como se ven los histogramas de los elementos no diagonales para las matrices anteriores (la varianza de la distribución aumenta monotónicamente):

---

Actualización: uso de factores aleatorios

Un método realmente sencillo para generar matrices de correlación aleatorias con algunas correlaciones fuertes fue utilizado en la respuesta por @shabbychef, y me gustaría ilustrarlo aquí también. La idea es generar aleatoriamente varias ( $k<d$ ) cargas factoriales $\mathbf W$ (matriz aleatoria de $k \times d$ tamaño), forman la matriz de covarianza $\mathbf W \mathbf W^\top$ (que por supuesto no será de rango completo) y añadirle una matriz diagonal aleatoria $\mathbf D$ con elementos positivos para hacer $\mathbf B = \mathbf W \mathbf W^\top + \mathbf D$ rango completo. La matriz de covarianza resultante puede normalizarse para convertirse en una matriz de correlación, dejando que $\mathbf C = \mathbf E^{-1/2}\mathbf B \mathbf E^{-1/2}$ , donde $\mathbf E$ es una matriz diagonal con la misma diagonal que $\mathbf B$ . Esto es muy simple y hace el truco. Aquí hay algunos ejemplos de matrices de correlación para $k={100, 50, 20, 10, 5, 1}$ :

Y el código:

%// FACTOR method
function S = factor(d,k)
    W = randn(d,k);
    S = W*W' + diag(rand(1,d));
    S = diag(1./sqrt(diag(S))) * S * diag(1./sqrt(diag(S)));
end

---

Este es el código de envoltura utilizado para generar las cifras:

d = 100; %// size of the correlation matrix

figure('Position', [100 100 1100 600])
for repetition = 1:6
    S = onion(d);

    %// etas = [1 10 100 1000 1e+4 1e+5];
    %// S = extendedOnion(d, etas(repetition));

    %// S = vine(d, etas(repetition));

    %// betaparams = [50 20 10 5 2 1];
    %// S = vineBeta(d, betaparams(repetition));

    subplot(2,3,repetition)

    %// use this to plot colormaps of S
    imagesc(S, [-1 1])
    axis square
    title(['Eigs: ' num2str(min(eig(S)),2) '...' num2str(max(eig(S)),2) ', det=' num2str(det(S),2)])

    %// use this to plot histograms of the off-diagonal elements
    %// offd = S(logical(ones(size(S))-eye(size(S))));
    %// hist(offd)
    %// xlim([-1 1])
end

Answer 2

Se puede hacer al revés: cada matriz $C \in \mathbb{R}_{++}^p$ (el conjunto de todos los simétricos $p \times p$ Las matrices PSD) pueden descomponerse como

$C=O^{T}DO$ donde $O$ es una matriz ortonormal

Para conseguir $O$ , primero genera una base aleatoria $(v_1,...,v_p)$ (donde $v_i$ son vectores aleatorios, normalmente en $(-1,1)$ ). A partir de ahí, utilice el proceso de ortogonalización de Gram-Schmidt para obtener $(u_1,....,u_p)=O$

$R$ tiene una serie de paquetes que pueden hacer la ortogonalización G-S de una base aleatoria de manera eficiente, es decir, incluso para grandes dimensiones, por ejemplo el paquete 'far'. Aunque encontrarás el algoritmo G-S en la wiki, probablemente sea mejor no reinventar la rueda e ir a por una implementación en matlab (seguro que existe, pero no puedo recomendar ninguna).

Finalmente, $D$ es una matriz diagonal cuyos elementos son todos positivos (esto es, de nuevo, fácil de generar: generar $p$ números al azar, cuadrarlos, ordenarlos y colocarlos en la diagonal de una identidad $p$ por $p$ matriz).

Answer 3

Una caracterización aún más sencilla es que para la matriz real $A$ , $A^TA$ es semidefinido positivo. Para ver por qué esto es así, sólo hay que demostrar que $y^T (A^TA) y \ge 0$ para todos los vectores $y$ (del tamaño adecuado, por supuesto). Esto es trivial: $y^T (A^TA) y = (Ay)^T Ay = ||Ay||$ que no es negativo. Así que en Matlab, simplemente intente

A = randn(m,n);   %here n is the desired size of the final matrix, and m > n
X = A' * A;

Dependiendo de la aplicación, esto puede no darle la distribución de valores propios que desea; la respuesta de Kwak es mucho mejor en ese sentido. Los valores propios de X producido por este fragmento de código debe seguir la distribución Marchenko-Pastur.

Para simular las matrices de correlación de las acciones, por ejemplo, es posible que desee un enfoque ligeramente diferente:

k = 7;      % # of latent dimensions;
n = 100;    % # of stocks;
A = 0.01 * randn(k,n);  % 'hedgeable risk'
D = diag(0.001 * randn(n,1));   % 'idiosyncratic risk'
X = A'*A + D;
ascii_hist(eig(X));    % this is my own function, you do a hist(eig(X));
-Inf <= x <  -0.001 : **************** (17)
-0.001 <= x <   0.001 : ************************************************** (53)
 0.001 <= x <   0.002 : ******************** (21)
 0.002 <= x <   0.004 : ** (2)
 0.004 <= x <   0.005 :  (0)
 0.005 <= x <   0.007 : * (1)
 0.007 <= x <   0.008 : * (1)
 0.008 <= x <   0.009 : *** (3)
 0.009 <= x <   0.011 : * (1)
 0.011 <= x <     Inf : * (1)

Answer 4

Como variación de la respuesta de kwak: generar una matriz diagonal $\mathbf{D}$ con valores propios aleatorios no negativos de una distribución de su elección, y luego realizar una transformación de similitud $\mathbf{A}=\mathbf{Q}\mathbf{D}\mathbf{Q}^T$ avec $\mathbf{Q}$ a Matriz ortogonal pseudoaleatoria distribuida por Haar .

Answer 5

No has especificado una distribución para las matrices. Dos de las más comunes son las distribuciones Wishart y Wishart inversa. La Descomposición de Bartlett da una factorización Cholesky de una matriz Wishart aleatoria (que también puede resolverse eficazmente para obtener una matriz Wishart inversa aleatoria).

De hecho, el espacio Cholesky es una forma conveniente de generar otros tipos de matrices aleatorias PSD, ya que sólo hay que asegurarse de que la diagonal es no negativa.