Algoritma PCA terbaik untuk sejumlah besar fitur (> 10K)?

Saya sebelumnya menanyakan ini pada StackOverflow, tetapi sepertinya itu mungkin lebih tepat di sini, mengingat bahwa itu tidak mendapatkan jawaban pada SO. Ini semacam persimpangan antara statistik dan pemrograman.

Saya perlu menulis beberapa kode untuk melakukan PCA (Principal Component Analysis). Saya telah melihat-lihat algoritma yang terkenal dan menerapkan yang ini , yang sejauh yang saya tahu sama dengan algoritma NIPALS. Ini bekerja dengan baik untuk menemukan 2-3 komponen utama pertama, tetapi kemudian tampaknya menjadi sangat lambat untuk bertemu (pada urutan ratusan hingga ribuan iterasi). Berikut ini rincian yang saya butuhkan:

1. Algoritma harus efisien ketika berhadapan dengan sejumlah besar fitur (pesanan 10.000 hingga 20.000) dan ukuran sampel pada urutan beberapa ratus.

2. Itu harus dapat dilaksanakan secara wajar tanpa perpustakaan aljabar / matriks linier yang layak, karena bahasa targetnya adalah D, yang belum memilikinya, dan bahkan jika itu terjadi, saya lebih suka tidak menambahkannya sebagai ketergantungan pada proyek yang dimaksud. .

Sebagai catatan, pada dataset yang sama R tampaknya menemukan semua komponen utama sangat cepat, tetapi menggunakan dekomposisi nilai singular, yang bukan sesuatu yang saya ingin kode sendiri.

Saya telah mengimplementasikan SVD Acak seperti yang diberikan dalam "Halko, N., Martinsson, PG, Shkolnisky, Y., & Tygert, M. (2010). Algoritma untuk analisis komponen utama dari kumpulan data besar. Arxiv pracetak arXiv: 1007.5510, 0526. Diperoleh 1 April 2011, dari http://arxiv.org/abs/1007.5510 . " Jika Anda ingin mendapatkan SVD terpotong, itu benar-benar bekerja jauh lebih cepat daripada variasi svd di MATLAB. Anda bisa mendapatkannya di sini:

function [U,S,V] = fsvd(A, k, i, usePowerMethod)
% FSVD Fast Singular Value Decomposition 
% 
%   [U,S,V] = FSVD(A,k,i,usePowerMethod) computes the truncated singular
%   value decomposition of the input matrix A upto rank k using i levels of
%   Krylov method as given in [1], p. 3.
% 
%   If usePowerMethod is given as true, then only exponent i is used (i.e.
%   as power method). See [2] p.9, Randomized PCA algorithm for details.
% 
%   [1] Halko, N., Martinsson, P. G., Shkolnisky, Y., & Tygert, M. (2010).
%   An algorithm for the principal component analysis of large data sets.
%   Arxiv preprint arXiv:1007.5510, 0526. Retrieved April 1, 2011, from
%   http://arxiv.org/abs/1007.5510. 
%   
%   [2] Halko, N., Martinsson, P. G., & Tropp, J. A. (2009). Finding
%   structure with randomness: Probabilistic algorithms for constructing
%   approximate matrix decompositions. Arxiv preprint arXiv:0909.4061.
%   Retrieved April 1, 2011, from http://arxiv.org/abs/0909.4061.
% 
%   See also SVD.
% 
%   Copyright 2011 Ismail Ari, http://ismailari.com.

    if nargin < 3
        i = 1;
    end

    % Take (conjugate) transpose if necessary. It makes H smaller thus
    % leading the computations to be faster
    if size(A,1) < size(A,2)
        A = A';
        isTransposed = true;
    else
        isTransposed = false;
    end

    n = size(A,2);
    l = k + 2;

    % Form a real n×l matrix G whose entries are iid Gaussian r.v.s of zero
    % mean and unit variance
    G = randn(n,l);


    if nargin >= 4 && usePowerMethod
        % Use only the given exponent
        H = A*G;
        for j = 2:i+1
            H = A * (A'*H);
        end
    else
        % Compute the m×l matrices H^{(0)}, ..., H^{(i)}
        % Note that this is done implicitly in each iteration below.
        H = cell(1,i+1);
        H{1} = A*G;
        for j = 2:i+1
            H{j} = A * (A'*H{j-1});
        end

        % Form the m×((i+1)l) matrix H
        H = cell2mat(H);
    end

    % Using the pivoted QR-decomposiion, form a real m×((i+1)l) matrix Q
    % whose columns are orthonormal, s.t. there exists a real
    % ((i+1)l)×((i+1)l) matrix R for which H = QR.  
    % XXX: Buradaki column pivoting ile yapılmayan hali.
    [Q,~] = qr(H,0);

    % Compute the n×((i+1)l) product matrix T = A^T Q
    T = A'*Q;

    % Form an SVD of T
    [Vt, St, W] = svd(T,'econ');

    % Compute the m×((i+1)l) product matrix
    Ut = Q*W;

    % Retrieve the leftmost m×k block U of Ut, the leftmost n×k block V of
    % Vt, and the leftmost uppermost k×k block S of St. The product U S V^T
    % then approxiamtes A. 

    if isTransposed
        V = Ut(:,1:k);
        U = Vt(:,1:k);     
    else
        U = Ut(:,1:k);
        V = Vt(:,1:k);
    end
    S = St(1:k,1:k);
end

Untuk mengujinya, cukup buat gambar di folder yang sama (seperti halnya matriks besar, Anda bisa membuat matriks sendiri)

% Example code for fast SVD.

clc, clear

%% TRY ME
k = 10; % # dims
i = 2;  % # power
COMPUTE_SVD0 = true; % Comment out if you do not want to spend time with builtin SVD.

% A is the m×n matrix we want to decompose
A = im2double(rgb2gray(imread('test_image.jpg')))';

%% DO NOT MODIFY
if COMPUTE_SVD0
    tic
    % Compute SVD of A directly
    [U0, S0, V0] = svd(A,'econ');
    A0 = U0(:,1:k) * S0(1:k,1:k) * V0(:,1:k)';
    toc
    display(['SVD Error: ' num2str(compute_error(A,A0))])
    clear U0 S0 V0
end

% FSVD without power method
tic
[U1, S1, V1] = fsvd(A, k, i);
toc
A1 = U1 * S1 * V1';
display(['FSVD HYBRID Error: ' num2str(compute_error(A,A1))])
clear U1 S1 V1

% FSVD with power method
tic
[U2, S2, V2] = fsvd(A, k, i, true);
toc
A2 = U2 * S2 * V2';
display(['FSVD POWER Error: ' num2str(compute_error(A,A2))])
clear U2 S2 V2

subplot(2,2,1), imshow(A'), title('A (orig)')
if COMPUTE_SVD0, subplot(2,2,2), imshow(A0'), title('A0 (svd)'), end
subplot(2,2,3), imshow(A1'), title('A1 (fsvd hybrid)')
subplot(2,2,4), imshow(A2'), title('A2 (fsvd power)')

Ketika saya menjalankannya di desktop saya untuk gambar ukuran 635 * 483, saya dapatkan

Elapsed time is 0.110510 seconds.
SVD Error: 0.19132
Elapsed time is 0.017286 seconds.
FSVD HYBRID Error: 0.19142
Elapsed time is 0.006496 seconds.
FSVD POWER Error: 0.19206

Seperti yang Anda lihat, untuk nilai rendah k, itu lebih dari 10 kali lebih cepat daripada menggunakan Matlab SVD. Omong-omong, Anda mungkin memerlukan fungsi sederhana berikut untuk fungsi tes:

function e = compute_error(A, B)
% COMPUTE_ERROR Compute relative error between two arrays

    e = norm(A(:)-B(:)) / norm(A(:));
end

Saya tidak menambahkan metode PCA karena mudah untuk menerapkan menggunakan SVD. Anda dapat memeriksa tautan ini untuk melihat hubungannya.

Anda dapat mencoba menggunakan beberapa opsi.

1- Penalized Matrix penguraian . Anda menerapkan beberapa batasan penalti pada u dan v untuk mendapatkan sparsity. Algoritma cepat yang telah digunakan pada data genomik

Lihat Whitten Tibshirani. Mereka juga memiliki R-pkg. "Dekomposisi matriks yang dihukum, dengan aplikasi untuk jarang komponen utama dan analisis korelasi kanonik."

2- SVD Acak . Karena SVD adalah algoritma utama, cari perkiraan yang sangat cepat mungkin diinginkan, terutama untuk analisis eksplorasi. Menggunakan SVD acak, Anda dapat melakukan PCA pada kumpulan data besar.

Lihat Martinsson, Rokhlin, dan Tygert "Algoritma acak untuk dekomposisi matriks". Tygert memiliki kode untuk implementasi PCA yang sangat cepat.

Di bawah ini adalah implementasi sederhana SVD acak di R.

ransvd = function(A, k=10, p=5) {
  n = nrow(A)
  y = A %*% matrix(rnorm(n * (k+p)), nrow=n)
  q = qr.Q(qr(y))
  b = t(q) %*% A
  svd = svd(b)
  list(u=q %*% svd$u, d=svd$d, v=svd$v)
}

Sepertinya Anda ingin menggunakan Algoritma Lanczos . Jika gagal, Anda mungkin ingin berkonsultasi dengan Golub & Van Loan. Saya pernah mengkodekan algoritma SVD (dalam SML, semua bahasa) dari teks mereka, dan itu bekerja dengan cukup baik.

Saya sarankan mencoba PCA kernel yang memiliki kompleksitas waktu / ruang tergantung pada jumlah contoh (N) daripada jumlah fitur (P), yang saya pikir akan lebih cocok di pengaturan Anda (P >> N)). Kernel PCA pada dasarnya bekerja dengan matriks kernel NxN (matriks kemiripan antara titik-titik data), dan bukan matriks kovarian PxP yang sulit diatasi untuk P. besar. Hal baik lainnya tentang kernel PCA adalah dapat mempelajari proyeksi non-linear. juga jika Anda menggunakannya dengan kernel yang sesuai. Lihat makalah ini pada kernel PCA .

Saya sepertinya ingat bahwa adalah mungkin untuk melakukan PCA dengan menghitung eigen-dekomposisi X ^ TX daripada XX ^ T dan kemudian berubah untuk mendapatkan PC. Namun saya tidak dapat mengingat detailnya secara langsung, tetapi itu ada di buku Jolliffe (luar biasa) dan saya akan mencarinya ketika saya berikutnya bekerja. Saya akan transliterasi rutin aljabar linier dari misalnya Metode Numerik dalam C, daripada menggunakan algoritma lain.

Ada juga metode bootstrap oleh Fisher et al , yang dirancang untuk beberapa ratus sampel berdimensi tinggi.

Gagasan utama dari metode ini dirumuskan sebagai "resampling adalah transformasi dimensi rendah". Jadi, jika Anda memiliki sejumlah kecil (beberapa ratus) sampel berdimensi tinggi, maka Anda tidak bisa mendapatkan lebih banyak komponen utama daripada jumlah sampel Anda. Dengan demikian masuk akal untuk mempertimbangkan sampel sebagai dasar pelit, memproyeksikan data pada subruang linier yang direntang oleh vektor-vektor ini, dan menghitung PCA dalam subruang yang lebih kecil ini. Mereka juga memberikan rincian lebih lanjut bagaimana menangani kasing ketika tidak semua sampel dapat disimpan dalam memori.

Lihat makalah Sam Roweis, Algoritma EM untuk PCA dan SPCA .