Korelasi aneh dalam hasil SVD dari data acak; apakah mereka memiliki penjelasan matematis atau apakah itu bug LAPACK?

Saya mengamati perilaku yang sangat aneh dalam hasil SVD dari data acak, yang dapat saya tiru di Matlab dan R. Sepertinya beberapa masalah numerik di perpustakaan LAPACK; Apakah itu?

Saya menarik $n=1000$ sampel dari $k=2$ dimensi Gaussian dengan nol mean dan kovarian identitas: $X\sim \mathcal N (0, \mathbf I)$ . Aku mengumpulkan mereka di $1000 \times 2$ matriks data $\mathbf X$ . (Saya dapat opsional pusat $\mathbf X$ atau tidak, itu tidak mempengaruhi berikut.) Lalu aku melakukan dekomposisi nilai singular (SVD) untuk mendapatkan $\mathbf X=\mathbf{USV}^\top$ . Mari kita ambil dua elemen tertentu dari $\mathbf U$ , misalnya $U_{11}$ dan $U_{22}$ , dan tanyakan apa korelasi di antara mereka di berbagai imbang berbeda$\mathbf X$. Saya akan berharap bahwa jika jumlah$N_\mathrm{rep}$ dari menarik cukup besar, maka semua korelasi tersebut akan menjadi sekitar nol (yaitu korelasi populasi harus nol, dan korelasi sampel akan menjadi kecil).

Namun, saya mengamati beberapa korelasi aneh yang kuat (sekitar $\pm0.2$ ) antara $U_{11}$ , $U_{12}$ , $U_{21}$ , dan $U_{22}$ , dan hanya di antara elemen-elemen ini. Semua pasangan elemen lainnya memiliki korelasi sekitar nol, seperti yang diharapkan. Berikut adalah bagaimana matriks korelasi untuk $20$ elemen "atas" dari $\mathbf U$ terlihat ( $10$ elemen pertama dari kolom pertama, kemudian $10$ elemen pertama dari kolom kedua):

Perhatikan nilai aneh tinggi di sudut kiri atas setiap kuadran.

Itu komentar ini @ whuber ini yang membawa efek ini ke perhatian saya. @whuber berpendapat bahwa PC1 dan PC2 tidak independen dan disajikan korelasi kuat ini sebagai bukti untuk itu. Namun, kesan saya adalah dia tidak sengaja menemukan bug numerik di perpustakaan LAPACK. Apa yang terjadi disini?

Berikut adalah kode R @ whuber:

stat <- function(x) {u <- svd(x)$u; c(u[1,1], u[2, 2])};
Sigma <- matrix(c(1,0,0,1), 2);
sim <- t(replicate(1e3, stat(MASS::mvrnorm(10, c(0,0), Sigma))));
cor.test(sim[,1], sim[,2]);

Ini kode Matlab saya:

clear all
rng(7)

n = 1000;     %// Number of variables
k = 2;        %// Number of observations
Nrep = 1000;  %// Number of iterations (draws)

for rep = 1:Nrep
    X = randn(n,k);
    %// X = bsxfun(@minus, X, mean(X));
    [U,S,V] = svd(X,0);

    t(rep,:) = [U(1:10,1)' U(1:10,2)'];
end

figure
imagesc(corr(t), [-.5 .5])
axis square
hold on
plot(xlim, [10.5 10.5], 'k')
plot([10.5 10.5], ylim, 'k')

Ini bukan bug.

Seperti yang telah kami jelajahi (ekstensif) dalam komentar, ada dua hal yang terjadi. Yang pertama adalah bahwa kolom $U$ dibatasi untuk memenuhi persyaratan SVD: masing-masing harus memiliki panjang unit dan ortogonal untuk yang lainnya. Melihat $U$ sebagai variabel acak yang dibuat dari matriks $X$ acak melalui algoritma SVD tertentu, dengan demikian kami mencatat bahwa k ini ($k(k+1)/2$ kendala independen fungsional ini menciptakan dependensi statistik di antara kolom$U$ .

Ketergantungan ini mungkin terungkap sampai tingkat yang lebih besar atau lebih kecil dengan mempelajari korelasi antara komponen $U$ , tetapi fenomena kedua muncul : solusi SVD tidak unik. Minimal, setiap kolom $U$ dapat dinegasikan secara independen, memberikan setidaknya $2^k$ solusi berbeda dengan $k$ kolom. Korelasi yang kuat (melebihi $1/2$ ) dapat diinduksi dengan mengubah tanda-tanda kolom tepat. (Salah satu cara untuk melakukan ini diberikan dalam komentar pertama saya untuk jawaban Amoeba di utas ini: Saya memaksa semua $u_{ii},i=1,\ldots, k$ untuk memiliki tanda yang sama, membuat semuanya negatif atau semuanya positif dengan probabilitas yang sama.) Di lain pihak, semua korelasi dapat dibuat menghilang dengan memilih tanda-tanda secara acak, mandiri, dengan probabilitas yang sama. (Saya memberikan contoh di bawah ini di bagian "Edit".)

Dengan hati-hati, kita dapat sebagian membedakan kedua fenomena ini ketika membaca matriks sebar komponen $U$ . Karakteristik tertentu - seperti penampakan titik-titik yang terdistribusi secara seragam di wilayah melingkar yang terdefinisi dengan baik - meyakini kurangnya kemandirian. Lainnya, seperti scatterplots yang menunjukkan korelasi bukan nol yang jelas, jelas tergantung pada pilihan yang dibuat dalam algoritma - tetapi pilihan tersebut hanya mungkin karena kurangnya independensi di tempat pertama.

Tes akhir dari algoritma penguraian seperti SVD (atau Cholesky, LR, LU, dll) adalah apakah ia melakukan apa yang diklaimnya. Dalam keadaan ini cukup untuk memeriksa bahwa ketika SVD mengembalikan tiga matriks $(U, D, V)$ , bahwa $X$ dipulihkan, hingga kesalahan floating point yang diantisipasi, oleh produk $UDV^\prime$ ; bahwa kolom $U$ dan $V$ adalah ortonormal; dan bahwa $D$ adalah diagonal, elemen-elemen diagonalnya adalah non-negatif, dan disusun dalam urutan menurun. Saya telah menerapkan tes tersebut ke svdalgoritma di Windows 7Rdan tidak pernah menemukan kesalahan. Meskipun itu bukan jaminan, itu sepenuhnya benar, pengalaman seperti itu - yang saya percaya dibagikan oleh banyak orang - menunjukkan bahwa bug apa pun akan membutuhkan beberapa jenis input yang luar biasa agar dapat terwujud.

Berikut ini adalah analisis yang lebih rinci tentang poin-poin spesifik yang diangkat dalam pertanyaan.

Menggunakan R's svdprosedur, pertama Anda dapat memeriksa bahwa sebagai $k$ meningkat, korelasi antara koefisien $U$ tumbuh lebih lemah, namun mereka masih nol. Jika Anda hanya melakukan simulasi yang lebih besar, Anda akan menemukan mereka signifikan. (Ketika $k=3$ , 50000 iterasi harus mencukupi.) Bertentangan dengan pernyataan dalam pertanyaan, korelasi tidak "menghilang sepenuhnya."

Kedua, cara yang lebih baik untuk mempelajari fenomena ini adalah dengan kembali ke pertanyaan mendasar tentang independensi koefisien. Meskipun korelasi cenderung mendekati nol dalam banyak kasus, kurangnya independensi jelas terlihat. Ini dibuat paling jelas dengan mempelajari distribusi multivariat penuh koefisien $U$ . Sifat distribusi muncul bahkan dalam simulasi kecil di mana korelasi nol tidak dapat (belum) terdeteksi. Misalnya, periksa matriks sebaran koefisien. Untuk membuat ini bisa dilakukan, saya mengatur ukuran setiap dataset yang disimulasikan menjadi $4$ dan tetap $k=2$ , sehingga menggambar $1000$realisasi dari matriks $4\times 2$$U$ , menciptakan matriks $1000\times 8$ . Berikut ini adalah matriks sebar penuhnya, dengan variabel-variabel yang dicantumkan berdasarkan posisi mereka di dalam $U$ :

Memindai bawah kolom pertama mengungkapkan kurangnya menarik kemerdekaan antara $u_{11}$ dan lainnya $u_{ij}$ : lihat bagaimana kuadran atas dari sebar dengan $u_{21}$ hampir kosong, misalnya; atau periksa awan miring ke atas yang menggambarkan hubungan $(u_{11}, u_{22})$ dan awan miring ke bawah untuk $(u_{21}, u_{12})$ . Pandangan yang dekat menunjukkan kurangnya kemandirian di antara hampir semua koefisien ini: sangat sedikit dari mereka yang terlihat independen, walaupun kebanyakan dari mereka menunjukkan korelasi yang hampir nol.

(NB: Sebagian besar awan melingkar adalah proyeksi dari hypersphere yang diciptakan oleh kondisi normalisasi yang memaksa jumlah kuadrat dari semua komponen dari setiap kolom menjadi satu.)

Matriks scatterplot dengan $k=3$ dan $k=4$ menunjukkan pola yang sama: fenomena ini tidak terbatas pada$k=2$ , juga tidak bergantung pada ukuran setiap set data yang disimulasikan: mereka semakin sulit untuk dihasilkan dan diperiksa.

Penjelasan untuk pola-pola ini mengacu pada algoritma yang digunakan untuk mendapatkan $U$ dalam dekomposisi nilai singular, tetapi kita tahu pola non-independensi tersebut harus ada oleh sifat-sifat yang sangat menentukan dari $U$ : karena setiap kolom berturut-turut adalah (secara geometris) ortogonal dengan pendahulunya yang, kondisi ortogonalitas ini memaksakan dependensi fungsional di antara koefisien, yang dengan demikian menerjemahkan dependensi statistik di antara variabel acak yang sesuai.

Edit

Menanggapi komentar, mungkin perlu berkomentar sejauh mana fenomena ketergantungan ini mencerminkan algoritma yang mendasari (untuk menghitung SVD) dan seberapa banyak mereka melekat dalam sifat proses.

The spesifik pola korelasi antara koefisien tergantung banyak pada pilihan sewenang-wenang yang dilakukan oleh algoritma SVD, karena solusinya adalah tidak unik: kolom $U$ dapat selalu independen dikalikan dengan $-1$ atau $1$ . Tidak ada cara intrinsik untuk memilih tanda. Jadi, ketika dua algoritma SVD membuat pilihan tanda yang berbeda (sewenang-wenang atau bahkan mungkin acak), mereka dapat menghasilkan pola pola sebar yang berbeda dari nilai $(u_{ij}, u_{i^\prime j^\prime})$ . Jika Anda ingin melihat ini, ganti statfungsi dalam kode di bawah ini dengan

stat <- function(x) {
  i <- sample.int(dim(x)[1]) # Make a random permutation of the rows of x
  u <- svd(x[i, ])$u         # Perform SVD
  as.vector(u[order(i), ])   # Unpermute the rows of u
}

Ini pertama secara acak memesan kembali pengamatan x, melakukan SVD, kemudian menerapkan urutan terbalik uagar sesuai dengan urutan pengamatan asli. Karena efeknya adalah untuk membentuk campuran versi yang dipantulkan dan diputar dari scatterplots asli, scatterplots dalam matriks akan terlihat jauh lebih seragam. Semua korelasi sampel akan sangat mendekati nol (dengan konstruksi: korelasi yang mendasarinya persis nol). Namun demikian, kurangnya independensi masih akan jelas (dalam bentuk lingkaran seragam yang muncul, khususnya antara $u_{i,j}$ dan $u_{i,j^\prime}$ ).

Kurangnya data di beberapa kuadran dari beberapa plot sebar asli (ditunjukkan pada gambar di atas) muncul dari bagaimana Ralgoritma SVD memilih tanda untuk kolom.

Tidak ada yang berubah tentang kesimpulan. Karena kolom kedua dari $U$ adalah ortogonal dengan yang pertama, maka (dianggap sebagai variabel acak multivariat) bergantung pada yang pertama (juga dianggap sebagai variabel acak multivariat). Anda tidak dapat membuat semua komponen dari satu kolom tidak tergantung pada semua komponen yang lain; yang dapat Anda lakukan adalah melihat data dengan cara yang mengaburkan dependensi - tetapi ketergantungan akan tetap ada.

Berikut adalah Rkode yang diperbarui untuk menangani kasing $k\gt 2$ dan menggambar sebagian dari matriks sebar.

k <- 2    # Number of variables
p <- 4    # Number of observations
n <- 1e3  # Number of iterations
stat <- function(x) as.vector(svd(x)$u)
Sigma <- diag(1, k, k); Mu <- rep(0, k)
set.seed(17)
sim <- t(replicate(n, stat(MASS::mvrnorm(p, Mu, Sigma))))
colnames(sim) <- as.vector(outer(1:p, 1:k, function(i,j) paste0(i,",",j)))
pairs(sim[, 1:min(11, p*k)], pch=".")

Jawaban ini menyajikan replikasi hasil @ whuber di Matlab, dan juga demonstrasi langsung bahwa korelasi adalah "artefak" bagaimana implementasi SVD memilih tanda untuk komponen.

Mengingat rantai panjang komentar yang berpotensi membingungkan, saya ingin menekankan untuk pembaca masa depan yang saya setujui sepenuhnya dengan yang berikut:

1. $\mathbf U$
2. $\mathbf U$$1$$U_{i1}$$U_{j1}$$i\ne j$$N_\mathrm{rep}$
3. $\mathbf U$$U_{i1}$$U_{j2}$

$\sim 0.2$$N_\mathrm{rep}=1000$

$n=4$$k=2$$N_\mathrm{rep}=1000$

Di sebelah kiri adalah matriks korelasi, di sebelah kanan - sebar plot mirip dengan @ whuber's. Kesepakatan antara simulasi kami tampaknya sempurna.

$U$

[U,S,V] = svd(X,0);

Saya menambahkan dua baris berikut:

U(:,1) = U(:,1) * sign(randn(1));
U(:,2) = U(:,2) * sign(randn(1));

Inilah hasilnya:

Semua korelasi menghilang, persis seperti yang saya harapkan dari awal !

$1$$1$

$\mathbf U$$\mathbf U$

PS. Selamat kepada @whuber karena telah melampaui reputasi 100rb hari ini!

Check the norm of your singular vectors U and V, it's 1 by definition. You don't need to go through SVD to get the same exact matrix you plot by simply generating two random variables $x$ and $y$ with the constraint that the sum of their squares is 1:

Assume that the means are zero, then

This will not be equal to zero. You can plug the standard normal into $x$ and $y$ to see what's the value of covariance to be expected here.