Чрезвычайно большое средневзвешенное значение

Я использую 64-битный Matlab с 32 г ОЗУ (просто так, чтобы вы знали).

У меня есть файл (вектор) из 1,3 миллиона чисел (целых). Я хочу сделать еще один вектор той же длины, где каждая точка является средневзвешенным значением всего первого вектора, взвешенным по обратному расстоянию от этой позиции (на самом деле это позиция ^-0,1, а не ^-1, но для примера) . Я не могу использовать функцию «фильтр» Matlab, потому что она может только усреднять вещи до текущей точки, верно? Чтобы объяснить более четко, вот пример из 3 элементов

data = [ 2 6 9 ]
weights = [ 1 1/2 1/3; 1/2 1 1/2; 1/3 1/2 1 ]
results=data*weights= [ 8 11.5 12.666 ]
i.e.
8 = 2*1 + 6*1/2 + 9*1/3
11.5 = 2*1/2 + 6*1 + 9*1/2
12.666 = 2*1/3 + 6*1/2 + 9*1

Таким образом, каждая точка в новом векторе является средневзвешенным значением всего первого вектора, взвешенным на 1/(расстояние от этой позиции + 1).

Я мог бы просто переделать вектор весов для каждой точки, а затем поэлементно вычислить вектор результатов, но для этого требуется 1,3 миллиона итераций цикла for, каждая из которых содержит 1,3 миллиона умножений. Я бы предпочел использовать прямое матричное умножение, умножая 1x1,3 мила на 1,3 мила x 1,3 мила, что теоретически работает, но я не могу загрузить матрицу такого размера.

Затем я пытаюсь создать матрицу с помощью сценария оболочки и проиндексировать ее в Matlab, чтобы одновременно вызывался только соответствующий столбец матрицы, но это также занимает очень много времени.

Мне не нужно делать это в Matlab, поэтому любые советы людей по поводу использования таких больших чисел и получения средних значений будут оценены. Поскольку я использую вес ^-0,1, а не ^-1, он не падает так быстро - миллионный балл по-прежнему имеет вес 0,25 по сравнению с исходным весом баллов 1, поэтому я не могу просто сократить его прочь, как он становится большим либо.

Надеюсь, это было достаточно ясно?

Вот код ответа ниже (так что его можно отформатировать?):

data = load('/Users/mmanary/Documents/test/insertion.txt');
data=data.';
total=length(data);
x=1:total;
datapad=[zeros(1,total) data];
weights = ([(total+1):-1:2 1:total]).^(-.4);
weights = weights/sum(weights);
Fdata = fft(datapad);
Fweights = fft(weights);
Fresults = Fdata .* Fweights;
results = ifft(Fresults);
results = results(1:total);
plot(x,results)

person Micah Manary    schedule 25.10.2011    source источник
comment
Всем, кто не знаком с подходом БПФ к решению подобных задач (см. мой ответ ниже), я настоятельно рекомендую бесплатный онлайн-ресурс Scientist и Руководство инженера по цифровой обработке сигналов. Даже для человека, не работающего с DSP, как я, знание основ бесценно. Ускорение при использовании БПФ по сравнению со сверткой грубой силы почти похоже на волшебство. :)   -  person John Colby    schedule 25.10.2011


Ответы (5)


Единственный разумный способ сделать это — использовать свертку БПФ, которая лежит в основе функции filter и подобных. Это очень легко сделать вручную:

% Simulate some data
n = 10^6;
x = randi(10,1,n);
xpad = [zeros(1,n) x];

% Setup smoothing kernel
k = 1 ./ [(n+1):-1:2 1:n];

% FFT convolution
Fx = fft(xpad);
Fk = fft(k);

Fxk = Fx .* Fk;

xk = ifft(Fxk);
xk = xk(1:n);

n=10^6 занимает менее полсекунды!

person John Colby    schedule 25.10.2011
comment
А для слишком длинных сигналов БПФ может разбиваться на кадры. - person Itamar Katz; 25.10.2011
comment
@MicahManary Отлично! Удачи тебе с твоим проектом. - person John Colby; 25.10.2011
comment
@JohnColby Итак, свертка работает, но как мне разделить на общий вес (чтобы получить фактическое среднее значение). Точки по краям нужно разделить на гораздо меньшее число, чем среднюю точку. Мой код сейчас в вопросе (поэтому его можно отформатировать). - person Micah Manary; 25.10.2011
comment
Я считаю, что просто масштабируйте ваше ядро ​​(т.е. weights) так, чтобы его сумма равнялась 1. - person John Colby; 25.10.2011
comment
@JohnColby Конечно, спасибо. Отредактированный код ответа выше, чтобы работать! - person Micah Manary; 25.10.2011

Это, вероятно, не лучший способ сделать это, но с большим количеством памяти вы определенно можете распараллелить процесс.

Вы можете построить разреженные матрицы, состоящие из элементов вашей исходной матрицы, которые имеют значение i^(-1) (где i = 1 .. 1.3 million), умножить их на исходный вектор и суммировать все результаты вместе.

Итак, для вашего примера продукт будет по существу:

a = rand(3,1);
b1 = [1 0 0;
      0 1 0;
      0 0 1];
b2 = [0 1 0;
      1 0 1;
      0 1 0] / 2;
b3 = [0 0 1;
      0 0 0;
      1 0 0] / 3;

c = sparse(b1) * a + sparse(b2) * a + sparse(b3) * a;

Конечно, вы бы не построили разреженные матрицы таким образом. Если вы хотите иметь меньше итераций внутреннего цикла, вы можете иметь более одного из i в каждой матрице.

Посмотрите на цикл parfor в MATLAB: http://www.mathworks.com/help/toolbox/distcomp/parfor.html

person yuzeh    schedule 25.10.2011

Я не могу использовать функцию «фильтр» Matlab, потому что она может только усреднять вещи до текущей точки, верно?

Это неправильно. Вы всегда можете добавить выборки (т. е. добавить или удалить нули) из своих данных или из отфильтрованных данных. Поскольку фильтрация с помощью filter (можно также использовать conv кстати) является линейным действием, она не изменит результат (это как добавление и удаление нулей, которое ничего не делает, а затем фильтрация. Тогда линейность позволяет вам поменять порядок добавить образцы -> фильтровать -> удалить образец).

В любом случае, в вашем примере вы можете использовать ядро ​​​​усреднения:

weights = 1 ./ [3 2 1 2 3]; % this kernel introduces a delay of 2 samples

а потом просто:

result =  filter(w,1,[data, zeros(1,3)]); % or conv (data, w)
% removing the delay introduced by the kernel
result = result (3:end-1);
person Itamar Katz    schedule 25.10.2011

Вы рассматривали только 2 варианта: Умножение матрицы 1.3M*1.3M на вектор один раз или умножение 2 векторов по 1.3M 1.3M раз.

Но вы можете разделить свою весовую матрицу на любое количество подматриц и выполнить умножение матрицы n * 1,3M на вектор 1,3M/n раз.

Я предполагаю, что быстрее всего будет, когда будет наименьшее количество итераций и n будет таким, что создаст самую большую подматрицу, которая помещается в вашей памяти, не заставляя ваш компьютер начинать подкачку страниц на ваш жесткий диск.

с вашим объемом памяти вы должны начать с n = 5000.

вы также можете сделать это быстрее, используя parfor (с делением n на количество процессоров).

person Tal Darom    schedule 25.10.2011

Путь грубой силы, вероятно, сработает для вас с одной небольшой оптимизацией.

Операции ^-0.1 для создания весов займут намного больше времени, чем операции + и * для вычисления взвешенных средних, но вы повторно используете веса во всех миллионах операций взвешенных средних. Алгоритм становится:

  • Создайте вектор весов со всеми весами, которые потребуются для любого вычисления: weights = (-n:n).^-0.1

  • Для каждого элемента вектора:

    • Индексируйте соответствующую часть вектора weights, чтобы считать текущий элемент «центром».

    • Выполните взвешенное среднее с частью весов и целым вектором. Это можно сделать с помощью быстрого векторного умножения с последующим скалярным делением.

Основной цикл выполняет n^2 сложений и вычитаний. Если n равно 1,3 миллиона, это 3,4 триллиона операций. Одно ядро ​​современного процессора с частотой 3 ГГц может выполнять, скажем, 6 миллиардов операций сложения/умножения в секунду, что составляет около 10 минут. Добавьте время для индексации вектора weights и накладных расходов, и, по моим оценкам, вы все еще можете прийти менее чем за полчаса.

person Wesley Hill    schedule 25.10.2011