Лучший алгоритм PCA для огромного количества функций (> 10K)?

54

Ранее я спрашивал об этом в StackOverflow, но кажется, что это может быть более уместным, учитывая, что он не получил никаких ответов по SO. Это своего рода на пересечении статистики и программирования.

Мне нужно написать код для PCA (Анализ основных компонентов). Я просмотрел известные алгоритмы и реализовал этот , который, насколько я могу судить, эквивалентен алгоритму NIPALS. Он хорошо работает для нахождения первых 2-3 главных компонентов, но затем кажется, что они очень медленно сходятся (порядка сотен и тысяч итераций). Вот детали того, что мне нужно:

  1. Алгоритм должен быть эффективен при работе с огромным количеством функций (порядка 10 000–20 000) и размерами выборки порядка нескольких сотен.

  2. Он должен быть разумно реализуемым без приличной библиотеки линейной алгебры / матрицы, поскольку целевым языком является D, у которого его еще нет, и даже если бы он был, я бы предпочел не добавлять его в качестве зависимости к рассматриваемому проекту. ,

Как примечание: в одном и том же наборе данных R, кажется, находит все основные компоненты очень быстро, но использует разложение по сингулярным значениям, которое я не хочу сам кодировать.

dsimcha
источник
2
Есть множество общедоступных алгоритмов SVD. См. En.wikipedia.org/wiki/… . Вы не можете использовать или адаптировать один из них? Кроме того, R с открытым исходным кодом, и под лицензией GPL, так почему бы не позаимствовать его алгоритм, если он выполняет свою работу?
Роб Хиндман
@Rob: Я хотел бы избежать практически написания библиотеки линейной алгебры, а также я хочу избежать авторского права GPL. Кроме того, я смотрел на кусочки исходного кода R и раньше, и, как правило, он не очень читабелен.
дсимча
4
Я что-то пропустил? У вас есть> 10K функций, но <1K образцов? Это означает, что последние компоненты 9K являются произвольными. Вы хотите все 1К первых компонентов?
Шаббычеф
2
В любом случае, вы не можете избежать использования SVD, хотя благодаря многим исследованиям по численной линейной алгебре, теперь есть много методов на выбор, в зависимости от того, насколько большой / маленький, разреженный / плотный ваша матрица, или если вам нужны только сингулярные значения или полный набор сингулярных значений и левый / правый сингулярные векторы. Алгоритмы не очень сложно понять, ИМХО.
JM не является статистиком
Можете ли вы сказать нам, почему вы хотите сделать PCA?
Робин Жирар

Ответы:

27

Я реализовал рандомизированный SVD, как указано в "Halko, N., Martinsson, PG, Shkolnisky, Y. & Tygert, M. (2010). Алгоритм анализа главных компонентов больших наборов данных. Препринт Arxiv arXiv: 1007.5510, 0526. Получено 1 апреля 2011 г. с сайта http://arxiv.org/abs/1007.5510 . " Если вы хотите получить усеченный SVD, он действительно работает намного быстрее, чем svd-варианты в MATLAB. Вы можете получить его здесь:

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

Чтобы проверить это, просто создайте изображение в той же папке (как большая матрица, вы можете создать матрицу самостоятельно)

% 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)')

Быстрая СВД

Когда я запускаю его на своем рабочем столе для изображения размером 635 * 483, я получаю

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

Как видите, при низких значениях kэто более чем в 10 раз быстрее, чем при использовании Matlab SVD. Кстати, вам может понадобиться следующая простая функция для тестовой функции:

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

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

Я не добавил метод PCA, поскольку его легко реализовать с использованием SVD. Вы можете проверить эту ссылку, чтобы увидеть их отношения.

Petrichor
источник
12

Вы могли бы попытаться использовать несколько вариантов.

1- штрафная матричная декомпозиция . Вы применяете некоторые штрафные ограничения к вам и v, чтобы получить некоторую разреженность. Быстрый алгоритм, который был использован на данных геномики

Смотрите Уиттен Тибширани. У них также есть R-pkg. «Наказанная матричная декомпозиция с приложениями для разреженных главных компонентов и канонического корреляционного анализа».

2- Рандомизированная СВД . Поскольку SVD является основным алгоритмом, может быть желательно найти очень быстрое приближение, особенно для исследовательского анализа. Используя рандомизированный SVD, вы можете делать PCA на огромных наборах данных.

См. Martinsson, Rokhlin, Tygert "Рандомизированный алгоритм разложения матриц". У Тигерта есть код для очень быстрой реализации PCA.

Ниже приведена простая реализация рандомизированного СВД в 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)
}
pslice
источник
+1 за штрафную декомпозицию матрицы. Этот пакет довольно удивительный. Я должен, вероятно, упомянуть, что это написано «Виттен», хотя, в случае, если люди не могут найти ссылку. Наконец, ОП сказал, что они не хотят, чтобы что-то было написано на R, но, по сути, любой большой SVD-пакет будет иметь C, C ++ или Fortran backend для скорости.
Дэвид Дж. Харрис
3

Я бы предложил попробовать ядро PCA, которое имеет временную / пространственную сложность, зависящую от количества примеров (N), а не от числа функций (P), что, я думаю, будет более подходящим в ваших настройках (P >> N)). Ядро PCA в основном работает с матрицей ядра NxN (матрицей сходства между точками данных), а не с ковариационной матрицей PxP, с которой может быть трудно иметь дело для больших P. Еще одна хорошая вещь в ядре PCA заключается в том, что он может изучать нелинейные проекции а также если вы используете его с подходящим ядром. Смотрите эту статью о ядре PCA .

ebony1
источник
2

Кажется, я вспоминаю, что можно выполнить PCA путем вычисления собственного разложения X ^ TX, а не XX ^ T, а затем преобразовать, чтобы получить ПК. Однако я не могу вспомнить подробности, но это в (превосходной) книге Джоллиффа, и я посмотрю ее, когда буду рядом на работе. Я бы транслитерал подпрограммы линейной алгебры, например, из численных методов в C, а не использовал бы любой другой алгоритм.

Дикран Сумчатый
источник
5
Боже мой ... построение ковариационной матрицы никогда не было лучшим способом для SVD. Я показал пример того, почему явное формирование ковариационной матрицы не очень хорошая идея в math.SE: math.stackexchange.com/questions/3869/3871#3871 .
JM не является статистиком
1

Существует также метод начальной загрузки Фишера и др. , Разработанный для нескольких сотен образцов высокой размерности.

Основная идея метода сформулирована как «повторная выборка - преобразование низкой размерности». Таким образом, если у вас есть небольшое (несколько сотен) количество высокоразмерных выборок, то вы не можете получить больше основных компонентов, чем количество ваших выборок. Таким образом, имеет смысл рассматривать выборки как экономную основу, проецировать данные на линейное подпространство, охватываемое этими векторами, и вычислять PCA в этом меньшем подпространстве. Они также предоставляют более подробную информацию о том, как справиться со случаем, когда не все образцы могут быть сохранены в памяти.

Коля Иванков
источник
0

См. Статью Сэма Роуиса « Алгоритмы EM для PCA и SPCA» .

АРС
источник
Алгоритм Википедии цитирует это и эквивалентно этому для случая нахождения одного главного компонента за один раз.
дсимча
Хорошо, теперь я вижу ссылку. Это довольно простой подход, и, как упоминает Википедия, эта базовая идея продвигается вперед. Поразмыслив, вам придется иметь дело с некоторыми компромиссами (в данном случае сходимостью). Интересно, если вы задаете правильный вопрос здесь. Неужели нет хороших привязок к библиотекам linalg для D?
АРС