¿Cómo generar de manera eficiente matrices de correlación aleatoria positiva-semidefinida?

Me gustaría poder generar eficientemente matrices de correlación semidefinidas positivas (PSD). Mi método se ralentiza dramáticamente a medida que aumento el tamaño de las matrices que se generarán.

1. ¿Podría sugerir alguna solución eficiente? Si conoce algún ejemplo en Matlab, estaría muy agradecido.
2. Al generar una matriz de correlación PSD, ¿cómo elegiría los parámetros para describir las matrices que se generarán? ¿Una correlación promedio, desviación estándar de correlaciones, valores propios?

Puede hacerlo al revés: cada matriz (el conjunto de todas las matrices PSD simétricas ) se puede descomponer como$C \in \mathbb{R}_{++}^p$$p \times p$

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

Para obtener , primero genere una base aleatoria (donde son vectores aleatorios, generalmente en ). A partir de ahí, use el proceso de ortogonalización de Gram-Schmidt para obtener( v 1 , . . . , V p ) v i ( - 1 , 1 ) ( u 1 , . . . . , U p ) = $O$$(v_1,...,v_p)$$v_i$$(-1,1)$$(u_1,....,u_p)=O$

$R$ tiene varios paquetes que pueden hacer la ortogonalización GS de una base aleatoria de manera eficiente, incluso para grandes dimensiones, por ejemplo, el paquete 'lejano'. Aunque encontrará el algoritmo GS en la wiki, probablemente sea mejor no reinventar la rueda y optar por una implementación de matlab (seguramente existe una, simplemente no puedo recomendar ninguna).

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

Un documento que genera matrices de correlación aleatorias basadas en vides y método de cebolla extendida por Lewandowski, Kurowicka y Joe (LKJ), 2009, proporciona un tratamiento unificado y una exposición de los dos métodos eficientes de generar matrices de correlación aleatorias. Ambos métodos permiten generar matrices a partir de una distribución uniforme en un cierto sentido preciso definido a continuación, son fáciles de implementar, rápidos y tienen la ventaja adicional de tener nombres divertidos.

Una matriz simétrica real de tamaño con unos en la diagonal tiene d ( d - 1 ) / 2 elementos únicos fuera de la diagonal y, por lo tanto, se puede parametrizar como un punto en R d ( d - 1 ) / 2 . Cada punto en este espacio corresponde a una matriz simétrica, pero no todos son positivos definidos (como deben ser las matrices de correlación). Las matrices de correlación, por lo tanto, forman un subconjunto de R d ( d - 1 ) /$d \times d$$d(d-1)/2$$\mathbb R^{d(d-1)/2}$$\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 de MATLAB de cada método y los ilustraré con .$d=100$

Método de cebolla

El método de la cebolla proviene de otro artículo (ref. # 3 en LKJ) y su nombre se debe al hecho de que las matrices de correlación se generan comenzando con una matriz y creciendo columna por columna y fila por 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:$1\times 1$

Aquí y debajo, el título de cada subtrama muestra 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 cebolla extendida

LKJ modifica ligeramente este método para poder muestrear matrices de correlación de una distribución proporcional a [ d e $\mathbf C$ . Cuanto mayor sea η , mayor será 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 η = 1 corresponde a una distribución uniforme. En la figura siguiente, las matrices se generan con η = 1 , 10 , 100 , 1000 , $[\mathrm{det}\:\mathbf C]^{\eta-1}$$\eta$$\eta=1$ .$\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 y no η = 1 (como afirma LKJ). No estoy seguro de dónde está el error.$\eta=0$$\eta=1$

%// 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 originalmente sugerido por Joe (J en LKJ) y mejorado por LKJ. Me gusta más, porque es conceptualmente más fácil y también más fácil de modificar. La idea es generar correlaciones parciales (son independientes y pueden tener cualquier valor de [ - 1 , 1 $d(d-1)/2$$[-1, 1]$sin restricciones) y luego conviértalos en correlaciones crudas a través de una fórmula recursiva. Es conveniente organizar el cálculo en un cierto orden, y este gráfico se conoce como "vid". Es importante destacar que si se muestrean correlaciones parciales de distribuciones beta particulares (diferentes para diferentes células en la matriz), la matriz resultante se distribuirá uniformemente. Aquí nuevamente, LKJ introduce un parámetro adicional para muestrear desde una distribución proporcional a [ d e $\eta$ . El resultado es idéntico a la cebolla extendida:$[\mathrm{det}\:\mathbf C]^{\eta-1}$

%// 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 vid con muestreo manual de correlaciones parciales.

$\pm 1$$[0,1]$$[-1, 1]$$\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.

Tenga en cuenta que en este caso no se garantiza que la distribución sea invariante de permutación, por lo que adicionalmente permuto aleatoriamente 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 los histogramas de los elementos fuera de la diagonal buscan las matrices anteriores (la variación de la distribución aumenta monotónicamente):

Actualización: utilizando factores aleatorios

$k<d$$\mathbf W$$k \times d$$\mathbf W \mathbf W^\top$$\mathbf D$$\mathbf B = \mathbf W \mathbf W^\top + \mathbf D$$\mathbf C = \mathbf E^{-1/2}\mathbf B \mathbf E^{-1/2}$$\mathbf E$$\mathbf B$$k={100, 50, 20, 10, 5, 1}$

Y el codigo:

%// 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

Aquí está el código de ajuste utilizado para generar las figuras:

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

$A$$A^TA$$y^T (A^TA) y \ge 0$$y$$y^T (A^TA) y = (Ay)^T Ay = ||Ay||$que no es negativo Entonces, 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 los valores propios que desea; La respuesta de Kwak es mucho mejor en ese sentido. Los valores propios Xproducidos por este fragmento de código deben seguir la distribución Marchenko-Pastur.

Para simular las matrices de correlación de 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)

$\mathbf{D}$$\mathbf{A}=\mathbf{Q}\mathbf{D}\mathbf{Q}^T$$\mathbf{Q}$

No ha especificado una distribución para las matrices. Dos comunes son las distribuciones Wishart y Wishart inversa. La descomposición de Bartlett proporciona una factorización de Cholesky de una matriz aleatoria de Wishart (que también se puede resolver de manera eficiente para obtener una matriz inversa aleatoria de Wishart).

De hecho, el espacio Cholesky es una forma conveniente de generar otros tipos de matrices PSD aleatorias, ya que solo debe asegurarse de que la diagonal no sea negativa.

$U^TSU$

Si desea tener más control sobre su matriz PSD simétrica generada, por ejemplo, generar un conjunto de datos de validación sintética, tiene una serie de parámetros disponibles. Una matriz PSD simétrica corresponde a una hiper-elipse en el espacio N-dimensional, con todos los grados de libertad relacionados:

1. Rotaciones
2. Longitudes de ejes.

Entonces, para una matriz bidimensional (es decir, elipse 2d), tendrá 1 rotación + 2 ejes = 3 parámetros.

Si las rotaciones traen a la mente matrices ortogonales, es un tren correcto, ya que la construcción es nuevamente $\Sigma=ODO^T$, con $\Sigma$ siendo la matriz Sym.PSD producida, $O$ la matriz de rotación (que es ortogonal) y $D$la matriz diagonal, cuyos elementos diagonales controlarán la longitud de los ejes de la elipse.

El siguiente código de Matlab traza 16 conjuntos de datos distribuidos en Gauss bidimensionales basados ​​en $\Sigma$, con un ángulo creciente. El código para la generación aleatoria de parámetros está en los comentarios.

figure;
mu = [0,0];
for i=1:16
    subplot(4,4,i)
    theta = (i/16)*2*pi;   % theta = rand*2*pi;
    U=[cos(theta), -sin(theta); sin(theta) cos(theta)];
    % The diagonal's elements control the lengths of the axes
    D = [10, 0; 0, 1]; % D = diag(rand(2,1));    
    sigma = U*D*U';
    data = mvnrnd(mu,sigma,1000);
    plot(data(:,1),data(:,2),'+'); axis([-6 6 -6 6]); hold on;
end

Para más dimensiones, la matriz Diagonal es directa (como arriba), y el $U$ debe derivar de la multiplicación de las matrices de rotación.

Un enfoque barato y alegre que he usado para probar es generar m N (0,1) n-vectores V [k] y luego usar P = d * I + Suma {V [k] * V [k] '} como una matriz nxn psd. Con m <n esto será singular para d = 0, y para d pequeño tendrá un número de condición alto.