Bagaimana cara secara efisien menghasilkan matriks korelasi positif-semidefinit acak?

38

Saya ingin dapat menghasilkan matriks korelasi positif-semidefinit (PSD) secara efisien. Metode saya melambat secara dramatis ketika saya meningkatkan ukuran matriks yang akan dihasilkan.

  1. Bisakah Anda menyarankan solusi yang efisien? Jika Anda mengetahui ada contoh di Matlab, saya akan sangat berterima kasih.
  2. Saat membuat matriks korelasi PSD bagaimana Anda memilih parameter untuk menggambarkan matriks yang akan dihasilkan? Korelasi rata-rata, standar deviasi korelasi, nilai eigen?
Eduardas
sumber

Jawaban:

16

Anda dapat melakukannya mundur: setiap matriks CR++p (himpunan semua matriks p×p PSD simetris ) dapat didekomposisi sebagai

C=OTDO mana O adalah matriks ortonormal

Untuk mendapatkan O , pertama-tama buat basis acak (v1,...,vp) (di mana vi adalah vektor acak, biasanya dalam (1,1) ). Dari sana, gunakan proses ortogonisasi Gram-Schmidt untuk mendapatkan (u1,....,up)=O

R memiliki sejumlah paket yang dapat melakukan GS orthogonalization secara acak secara efisien, yang bahkan untuk dimensi besar, misalnya paket 'jauh'. Meskipun Anda akan menemukan algoritma GS di wiki, mungkin lebih baik untuk tidak menemukan kembali roda dan pergi untuk implementasi matlab (satu pasti ada, saya tidak bisa merekomendasikan apapun).

Akhirnya, adalah matriks diagonal yang elemen-elemennya semuanya positif (ini, sekali lagi, mudah untuk menghasilkan: menghasilkan angka acak , kuadratinya, mengurutkannya dan menempatkannya ke diagonal identitas dengan matriks ).p p pDpphal

pengguna603
sumber
3
(1) Perhatikan bahwa dihasilkan tidak akan menjadi matriks korelasi (seperti yang diminta oleh OP), karena tidak akan memiliki matriks korelasi diagonal. Tentu saja hal itu dapat Rescaled untuk memiliki orang-orang di diagonal, dengan menetapkan ke , di mana adalah matriks diagonal dengan diagonal yang sama seperti . (2) Jika saya tidak salah, ini akan menghasilkan matriks korelasi di mana semua elemen off-diagonal terkonsentrasi di sekitar , sehingga tidak ada fleksibilitas yang dicari OP (OP ingin dapat mengatur "korelasi rata-rata , standar deviasi korelasi, nilai eigen " )E - 1 / 2 C E - 1 / 2 E C 0CE-1/2CE-1/2EC0
amoeba berkata Reinstate Monica
@amoeba: Saya akan membahas (2) karena, seperti yang Anda tunjukkan, solusi untuk (1) sepele. Karakterisasi satu angka dari 'bentuk' (hubungan antara masuk dan keluar dari elemen diagonal) dari matriks PSD (dan karenanya kovarians dan juga matriks korelasi) adalah nomor kondisinya. Dan, metode yang dijelaskan di atas memungkinkan kontrol penuh atasnya. 'Konsentrasi elemen diagonal off sekitar 0' bukan fitur metode yang digunakan untuk menghasilkan matriks PSD tetapi, melainkan, konsekuensi dari kendala yang diperlukan untuk memastikan bahwa matriksnya adalah PSD dan fakta bahwa adalah besar. hal
user603
Apakah Anda mengatakan bahwa semua matriks PSD besar memiliki elemen off-diagonal mendekati nol? Saya tidak setuju, tidak begitu. Lihat jawaban saya di sini untuk beberapa contoh: Bagaimana menghasilkan matriks korelasi acak yang memiliki sekitar entri off-diagonal yang terdistribusi normal dengan deviasi standar yang diberikan? Tetapi orang dapat langsung melihat bahwa itu tidak terjadi, karena sebuah matriks persegi memiliki semua yang ada di diagonal dan nilai tetap mana-mana di luar-diagonal adalah PSD dan dapat besar secara sewenang-wenang (tetapi tentu saja di bawah ). ρ 1ρρ1
Amuba mengatakan Reinstate Monica
@amoeba: maka saya salah dalam mengasumsikan bahwa dengan kebutuhan, diagonal dari matriks korelasi besar ketika mereka diizinkan baik positif maupun negatif mendekati 0. Terima kasih atas contoh yang mencerahkan.
user603
1
Saya membaca makalah yang sangat bagus tentang menghasilkan matriks korelasi acak dan memberikan jawaban saya sendiri di sini (serta jawaban lain di utas terkait). Saya pikir Anda mungkin menganggapnya menarik.
Amuba kata Reinstate Monica
27

Sebuah makalah Menghasilkan matriks korelasi acak berdasarkan tanaman merambat dan metode bawang yang diperluas oleh Lewandowski, Kurowicka, dan Joe (LKJ), 2009, memberikan pengobatan terpadu dan paparan dua metode yang efisien untuk menghasilkan matriks korelasi acak. Kedua metode memungkinkan untuk menghasilkan matriks dari distribusi yang seragam dalam arti tertentu yang didefinisikan di bawah ini, mudah diterapkan, cepat, dan memiliki keuntungan tambahan karena memiliki nama yang lucu.

Matriks simetris sungguhan ukuran dengan yang ada di diagonal memiliki elemen off-diagonal yang unik dan dengan demikian dapat ditentukan sebagai titik dalam . Setiap titik dalam ruang ini sesuai dengan matriks simetris, tetapi tidak semua dari mereka adalah positif-pasti (seperti matriks korelasi harus). Matriks korelasi membentuk subset dari (sebenarnya subset cembung yang terhubung), dan kedua metode dapat menghasilkan poin dari distribusi yang seragam pada subset ini.d ( d - 1 ) / 2 R d ( d - 1 ) / 2d×dd(d-1)/2Rd(d-1)/2Rd(d-1)/2

Saya akan memberikan implementasi MATLAB saya sendiri untuk setiap metode dan menggambarkannya dengan .d=100


Metode bawang

Metode bawang berasal dari makalah lain (ref # 3 di LKJ) dan memiliki namanya dengan fakta bahwa matriks korelasi dihasilkan mulai dengan matriks dan menumbuhkannya kolom demi kolom dan baris demi baris. Distribusi yang dihasilkan seragam. Saya tidak begitu mengerti matematika di balik metode ini (dan lebih memilih metode kedua), tetapi inilah hasilnya:1×1

Metode bawang

Di sini dan di bawah judul setiap subplot menunjukkan nilai eigen terkecil dan terbesar, dan penentu (produk dari semua nilai eigen). Ini kodenya:

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

Metode bawang merah yang diperluas

LKJ memodifikasi metode ini sedikit, agar dapat mengambil sampel matriks korelasi dari distribusi yang sebanding dengan . Semakin besar , semakin besar akan menjadi penentu, yang berarti bahwa matriks korelasi yang dihasilkan akan semakin mendekati matriks identitas. Nilai sesuai dengan distribusi seragam. Pada gambar di bawah ini, matriks dihasilkan dengan . [ d e tC η η = 1 η = 1 , 10 , 100 , 1000 , 10[detC]η-1ηη=1η=1,10,100,1000,10000,100000

Metode bawang merah yang diperluas

Untuk beberapa alasan untuk mendapatkan determinan dengan urutan yang sama besarnya seperti pada metode bawang vanili, saya perlu meletakkan dan bukan (seperti yang diklaim oleh LKJ). Tidak yakin di mana kesalahannya.η = 1η=0η=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

Metode anggur

Metode Vine pada awalnya disarankan oleh Joe (J dalam LKJ) dan ditingkatkan oleh LKJ. Saya lebih menyukainya, karena secara konsep lebih mudah dan juga lebih mudah untuk dimodifikasi. Idenya adalah untuk menghasilkan korelasi parsial (mereka independen dan dapat memiliki nilai dari tanpa kendala) dan kemudian mengubahnya menjadi korelasi mentah melalui rumus rekursif. Lebih mudah mengatur komputasi dalam urutan tertentu, dan grafik ini dikenal sebagai "anggur". Yang penting, jika korelasi parsial diambil dari distribusi beta tertentu (berbeda untuk sel yang berbeda dalam matriks), maka matriks yang dihasilkan akan didistribusikan secara seragam. Di sini, LKJ memperkenalkan parameter tambahan untuk sampel dari distribusi yang proporsional[ - 1 , 1 ] η [ d e td(d-1)/2[-1,1]η[detC]η-1 . Hasilnya identik dengan bawang merah yang diperluas:

Metode anggur

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

Metode Vine dengan pengambilan sampel manual korelasi parsial

Seperti yang dapat dilihat di atas, distribusi yang seragam menghasilkan matriks korelasi yang hampir diagonal. Tetapi orang dapat dengan mudah memodifikasi metode anggur untuk memiliki korelasi yang lebih kuat (ini tidak dijelaskan dalam makalah LKJ, tetapi langsung): untuk yang ini harus sampel korelasi parsial dari distribusi yang terkonsentrasi sekitar . Di bawah ini saya sampel dari distribusi beta (dihitung ulang dari menjadi ) dengan . Semakin kecil parameter distribusi beta, semakin terkonsentrasi di dekat tepi.[ 0 , 1 ] [ - 1 , 1 ] α = β = 50 , 20 , 10 , 5 , 2±1[0,1][-1,1]α=β=50,20,10,5,2,1

Metode Vine dengan pengambilan sampel manual

Perhatikan bahwa dalam hal ini distribusi tidak dijamin sebagai permutasi invarian, jadi saya juga secara acak mengubah baris dan kolom setelah generasi.

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

Berikut adalah bagaimana histogram elemen off-diagonal mencari matriks di atas (varians dari distribusi meningkat secara monoton):

Elemen off-diagonal


Pembaruan: menggunakan faktor acak

Salah satu metode yang sangat sederhana untuk menghasilkan matriks korelasi acak dengan beberapa korelasi kuat digunakan dalam jawaban oleh @ shabbychef, dan saya ingin mengilustrasikannya di sini juga. Idenya adalah untuk secara acak menghasilkan beberapa ( ) loadings faktor (matriks acak ukuran ), membentuk matriks kovarians (yang tentu saja tidak akan peringkat penuh ) dan tambahkan padanya matriks diagonal acak dengan elemen positif untuk membuat peringkat penuh. Matriks kovarians yang dihasilkan dapat dinormalisasi menjadi matriks korelasi, dengan membiarkank<dWk×dWWDB=WW+DC=E-1/2BE-1/2, Di mana adalah matriks diagonal dengan diagonal yang sama seperti . Ini sangat sederhana dan bermanfaat. Berikut adalah beberapa contoh matriks korelasi untuk :EBk=100,50,20,10,5,1

matriks korelasi acak dari faktor acak

Dan kodenya:

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

Berikut adalah kode pembungkus yang digunakan untuk menghasilkan angka:

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
amuba kata Reinstate Monica
sumber
2
Ini adalah jadwal yang fantastis, saya senang saya mengatakan sesuatu!
shadowtalker
Ketika saya menerjemahkan kode matlab untuk matriks korelasi berbasis anggur ke R dan mengujinya, kepadatan korelasi dalam kolom 1 selalu berbeda dengan kolom selanjutnya. Mungkin saya salah menerjemahkan sesuatu, tetapi mungkin catatan ini membantu seseorang.
Charlie
3
Untuk pengguna R, fungsi rcorrmatrix di package clusterGeneration (ditulis oleh W Qui dan H. Joe) mengimplementasikan metode vine.
RNM
15

Karakterisasi yang bahkan lebih sederhana adalah bahwa untuk matriks nyata , adalah semidefinit positif. Untuk melihat mengapa hal ini terjadi, kita hanya perlu membuktikan bahwa untuk semua vektor (dengan ukuran yang tepat, tentu saja). Ini sepele:yang tidak negatif. Jadi di Matlab, coba sajaSEBUAHy T ( A T A ) y 0 y y T ( A T A ) y = ( A y ) T A y = | | A y | |SEBUAHTSEBUAHyT(SEBUAHTSEBUAH)y0yyT(SEBUAHTSEBUAH)y=(SEBUAHy)TSEBUAHy=||SEBUAHy||

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

Bergantung pada aplikasinya, ini mungkin tidak memberi Anda distribusi nilai eigen yang Anda inginkan; Jawaban Kwak jauh lebih baik dalam hal itu. Nilai eigen yang Xdihasilkan oleh cuplikan kode ini harus mengikuti distribusi Marchenko-Pastur.

Untuk mensimulasikan matriks korelasi saham, katakanlah, Anda mungkin menginginkan pendekatan yang sedikit berbeda:

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)
shabbychef
sumber
1
apakah Anda bersedia berbagi fungsi ascii_hist Anda secara kebetulan?
btown
@ bandar marginnya terlalu kecil untuk menampungnya!
shabbychef
1
Ada salah ketik di- hilang persegi terakhirnya! yT(ATA)y=(Ay)TAy=||Ay||
Silverfish
8

Sebagai variasi jawaban kwak: hasilkan matriks diagonal dengan nilai eigen nonnegatif acak dari distribusi pilihan Anda, dan kemudian lakukan transformasi kemiripan dengan sebuah matriks orthogonal pseudorandom yang didistribusikan secara Haar .A = Q D Q T QDA=QDQTQ

JM bukan ahli statistik
sumber
M.: Referensi yang bagus: Ini tampaknya merupakan solusi yang paling efisien (tanpa gejala).
whuber
3
@whuber: Heh, saya mengambilnya dari Golub dan Van Loan (tentu saja); Saya menggunakan ini sepanjang waktu untuk membantu dalam menghasilkan matriks pengujian untuk rutinitas pengujian eigenvalue / nilai singular. Seperti yang dapat dilihat dari makalah, ini pada dasarnya setara dengan QR-dekomposisi matriks acak seperti apa yang disarankan kwak, kecuali bahwa itu dilakukan lebih efisien. Ada implementasi MATLAB ini di Higham's Text Matrix Toolbox, BTW.
JM bukan ahli statistik
M.:> Terima kasih atas implementasi matlab. Apakah Anda kebetulan tahu generator matriks pseudo-acak Haar di R?
user603
@ kwak: Tidak tahu, tetapi jika belum ada implementasi, seharusnya tidak terlalu sulit untuk menerjemahkan kode MATLAB ke R (saya bisa mencoba untuk membuat satu jika benar-benar tidak ada); satu-satunya prasyarat adalah generator yang layak untuk varian normal pseudorandom, yang saya yakin R miliki.
JM bukan ahli statistik
M.:> Ya saya mungkin akan menerjemahkannya sendiri. Terima kasih atas tautannya, Terbaik.
user603
4

Anda belum menentukan distribusi untuk matriks. Dua yang umum adalah distribusi Wishart dan terbalik Wishart. The Bartlett dekomposisi memberikan factorisation Cholesky dari matriks acak Wishart (yang juga dapat diselesaikan secara efisien untuk mendapatkan inverse acak Wishart matrix).

Faktanya, ruang Cholesky adalah cara yang mudah untuk menghasilkan jenis matriks PSD acak lainnya, karena Anda hanya perlu memastikan bahwa diagonalnya non-negatif.

Simon Byrne
sumber
> Tidak acak: dua matriks yang dihasilkan dari Whishard yang sama tidak akan independen satu sama lain. Jika Anda berencana untuk mengubah Whishart di setiap generasi, lalu bagaimana Anda berniat untuk menghasilkan Whishart di tempat pertama?
user603
@ kwak: Saya tidak mengerti pertanyaan Anda: dekomposisi Bartlett akan memberikan undian independen dari distribusi Wishart yang sama.
Simon Byrne
> Biarkan saya ulangi ini, dari mana Anda mendapatkan matriks skala dari distribusi whishart Anda?
user603
1
@ kwak: ini adalah parameter distribusi, dan sudah diperbaiki. Anda memilihnya di awal, berdasarkan karakteristik distribusi yang Anda inginkan (seperti rata-rata).
Simon Byrne
3

UTSU

gappy
sumber
Jika entri dihasilkan dari distribusi Normal dan bukan seragam, dekomposisi yang Anda sebutkan harus SO (n) invarian (dan karenanya sama-sama didistribusikan relatif terhadap ukuran Haar).
whuber
Menarik. Bisakah Anda memberikan referensi untuk ini?
gappy
1
> masalah dengan metode ini adalah Anda tidak dapat mengontrol rasio nilai eigen terkecil hingga terbesar (dan saya pikir karena ukuran dataset yang dibuat secara acak tidak terhingga, rasio ini akan menyatu dengan 1).
user603
1

Jika Anda ingin memiliki kontrol lebih besar atas matriks PSD simetris yang Anda hasilkan, mis. Menghasilkan dataset validasi sintetis, Anda memiliki sejumlah parameter yang tersedia. Matriks PSD simetris sesuai dengan hiper-elips di ruang dimensi-N, dengan semua derajat kebebasan yang terkait:

  1. Rotasi.
  2. Panjang sumbu.

Jadi, untuk matriks 2 dimensi (yaitu elips 2d), Anda akan memiliki 1 rotasi + 2 sumbu = 3 parameter.

Σ=HAIDHAITΣHAID

Σ

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

U

PeriRamm
sumber
0

Pendekatan murah dan ceria yang saya gunakan untuk pengujian adalah menghasilkan m N (0,1) n-vektor V [k] dan kemudian menggunakan P = d * I + Jumlah {V [k] * V [k] '} sebagai matriks psn nxn. Dengan m <n ini akan tunggal untuk d = 0, dan untuk d kecil akan memiliki angka kondisi tinggi.


sumber
2
> masalah dengan metode ini adalah Anda tidak dapat mengontrol rasio nilai eigen terkecil hingga terbesar (dan saya pikir karena ukuran dataset yang dibuat secara acak tidak terhingga, rasio ini akan menyatu dengan 1).
user603
> selain itu, metode ini tidak terlalu efisien (dari sudut pandang komputasi)
user603
1
"Matriks acak" Anda adalah matriks terstruktur khusus yang disebut "matriks diagonal plus peringkat-1" (matriks DR1), jadi bukan matriks acak representatif yang baik.
JM bukan ahli statistik