Форум программистов, компьютерный форум, киберфорум
Наши страницы
Matlab
Войти
Регистрация
Восстановить пароль
 
AndreyL
1 / 1 / 1
Регистрация: 11.09.2013
Сообщений: 50
#1

Можно ли векторизовать операцию, избавившись от цикла for?

14.02.2017, 07:29. Просмотров 384. Ответов 15
Метки нет (Все метки)

Дамы и Господа!

Возникла такая задача. Есть массив Х размерностью (n,m), два массива Y и V размерностью (k,m) и массив Р размерностью (k,1). Необходимо из каждой строки массива Х вычесть каждую строку массива Y, возвести это все в квадрат и поделить на строки массива V. После чего сложить элементы строк получившегося массива и сумму логарифмов строк массива V и прибавить к этому элемент массива Р. На выходе должен получится массив L размерностью (n,k). Довольно сумбурное объяснение, проще показать, как это сейчас делается в цикле:
Matlab M
1
2
3
4
5
for i=1:k
  xProm=bsxfun(@minus,X,Y(i,:)).^2;
  xProm=bsxfun(@rdivide,xProm,V(i,:));    
  L(:,i)=-0.5*(sum(xProm,2)+sum(log(V(i,:)),2))-sumconst+P(i);
end;
умножение на -0.5 и скалярное слагаемое sumconst в этой операции - это фактически вычисление логарифма произведения плотностей нормальных распределений с центрами Y, дисперсиями V и логарифмами априорных вероятностей Р. Разбивка на три строки сделана исключительно для наглядности, можно все записать в одну строку, смысл от этого не изменится.

Вопрос такой - можно ли векторизовать эту операцию, избавившись от цикла for? В этом цикле алгоритм живет львиную долю времени - хотелось бы его разогнать.
0
Надоела реклама? Зарегистрируйтесь и она исчезнет полностью.
Similar
Эксперт
41792 / 34177 / 6122
Регистрация: 12.04.2006
Сообщений: 57,940
14.02.2017, 07:29
Ответы с готовыми решениями:

Как векторизовать этот цикл?
Всем привет. Есть проблема с векторизацией одного цикла в программке,...

Можно ли оптимизировать код, избавившись от его повторяющихся (похожих) участков
Здравствуйте, помогите пожалуйста расписать код что сначыт каждая строчка, те...

Показать, что любой оператор цикла while можно записать с помощью цикла repeat
Показать, что любой оператор цикла с предусловием можно записать с помощью...

Векторизовать размытую линию
Доброго времени! Собственно сабж. Проблема в том что линия на картинке не...

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

15
Nick07
415 / 329 / 35
Регистрация: 17.07.2013
Сообщений: 1,729
14.02.2017, 08:54 #2
Цитата Сообщение от AndreyL Посмотреть сообщение
цикле алгоритм живет львиную долю времени - хотелось бы его разогнать.
Это хорошо выглядит на простеньких примера из документации и учебников.
В реалии вопрос не однозначный. Где-то векторизация ускоряет, но у меня были ситуации, когда векторизация тормозила.
Однозначное позитивное решение - аппаратное: увеличить ОЗУ, применить SSD. Иногда хороший эффект дает GPU.
Все определяется объемом и структурой данных.
См:
Altman Y. Accelerating MATLAB 2015

Добавлено через 1 минуту
ЕХЕ-шники тоже могут дать эффект ускорения.

Добавлено через 2 минуты
Для начала выполните общепринятые правила:
Стр. 243:
_Дьяконов А Г АНАЛИЗ ДАННЫХ MATLAB ВМК МГУ 2010 г 278 стр
0
AndreyL
1 / 1 / 1
Регистрация: 11.09.2013
Сообщений: 50
14.02.2017, 10:41  [ТС] #3
Цитата Сообщение от Nick07 Посмотреть сообщение
Однозначное позитивное решение - аппаратное: увеличить ОЗУ, применить SSD. Иногда хороший эффект дает GPU.
8-ми ядерный Атлон, при счете задействуются только 3 ядра, и то на 30%, общая загрузка процессора 12%, когда не считается не более 1%. Оперативки 16 ГБ, общее задействование памяти менее 10%. Про GPU, конечно, здорово, но покупать ради этого Теслу мне никто не даст.
Цитата Сообщение от Nick07 Посмотреть сообщение
Все определяется объемом и структурой данных.
Примерно так n=70-100, m=8-12, k=5-9.
Цитата Сообщение от Nick07 Посмотреть сообщение
ЕХЕ-шники тоже могут дать эффект ускорения.
Не хотелось бы связываться с ЕХЕ-шниками - задача меняется постоянно, постоянно вводятся коррективы в код.
Цитата Сообщение от Nick07 Посмотреть сообщение
Для начала выполните общепринятые правила:
Место под L выделено заранее, можно, конечно, сумму логарифмов заменить логарифмом произведения, но это не сильно спасет. К тому же даже Дьяконов пишет "Попробуйте несколько вариантов решения одной задачи и выберите оптимальный", вот и хочу заменить цикл векторной операцией. Если, конечно, это вообще возможно.
0
Nick07
415 / 329 / 35
Регистрация: 17.07.2013
Сообщений: 1,729
14.02.2017, 12:04 #4
Цитата Сообщение от AndreyL Посмотреть сообщение
Примерно так n=70-100, m=8-12, k=5-9.
И сколько времени считает?
0
AndreyL
1 / 1 / 1
Регистрация: 11.09.2013
Сообщений: 50
14.02.2017, 13:38  [ТС] #5
Цитата Сообщение от Nick07 Посмотреть сообщение
И сколько времени считает?
Единичный вызов этого цикла 0.32 секунды, таких вызовов на один расчет примерно 60 тысяч (моделирование Монте-Карлой с ЕМ окончанием). Итого 5 с копейками часов. На все остальное остается 30% времени
0
Nick07
415 / 329 / 35
Регистрация: 17.07.2013
Сообщений: 1,729
14.02.2017, 17:18 #6
А профилирование по времени смотрели, на что идет основное время?
0
Excalibur921
750 / 425 / 68
Регистрация: 12.10.2013
Сообщений: 2,837
14.02.2017, 18:40 #7
Цитата Сообщение от Nick07 Посмотреть сообщение
ЕХЕ-шники тоже могут дать эффект ускорения.
Попробуйте скомпилировать. Это многократно ускорит.
Цитата Сообщение от AndreyL Посмотреть сообщение
Не хотелось бы связываться с ЕХЕ-шниками - задача меняется постоянно
Каждый раз и компилить а как иначе? Помню в mathematica рисовал поверхность думал компиляция это ерунда…ан нет =).
И это был еще не exe а некий внутренний компилятор. Exe было лень разбирать итак пахало в разы быстрей голого кода.

Как вариант вашей сортировки\ даже всей проги на языке программирования с нуля с вкуриванием как включить все ядра.
0
Nick07
415 / 329 / 35
Регистрация: 17.07.2013
Сообщений: 1,729
14.02.2017, 20:44 #8
Цитата Сообщение от Excalibur921 Посмотреть сообщение
с вкуриванием как включить все ядра.
Очень здравая мысль, хотя и не матлабовская.
0
Centurio
Модератор
724 / 653 / 187
Регистрация: 13.09.2015
Сообщений: 2,355
14.02.2017, 21:20 #9
AndreyL, во-первых, поправлю вас, вы путаете размерность и размер массива. Почти все ваши массивы имеют размерность 2, и только массив Р - размерность 1. А то, что вы называете размерностью - это на самом деле размер.
Во-вторых, векторизацию провести, конечно можно. Вот вам код, который считает гораздо быстрее цикла:
Matlab M
1
2
3
4
5
6
7
8
9
[n,m]=size(X);
k=size(Y,1);
Xe=repmat(X,k,1);
Ye=repelem(Y,n,1);
Ve=repelem(V,n,1);
Pe=repelem(P,n,1);
xProm=(Xe-Ye).^2./Ve;
Le=-0.5*(sum(xProm,2)+log(prod(Ve,2)))-sumconst+Pe;
L=reshape(Le,n,k);
0
AndreyL
1 / 1 / 1
Регистрация: 11.09.2013
Сообщений: 50
19.02.2017, 10:16  [ТС] #10
Цитата Сообщение от Centurio Посмотреть сообщение
поправлю вас, вы путаете размерность и размер массива
Да, конечно, Спасибо!
Цитата Сообщение от Centurio Посмотреть сообщение
Вот вам код, который считает гораздо быстрее цикла
Спасибо, ускорение раза в три. Пришлось, правда, обновить МАТЛАБ (я работал на 2013, который меня до сих пор устраивал), поэтому долго молчал. Идея развернуть трехмерный массив в двумерный мне не приходила.
0
Excalibur921
750 / 425 / 68
Регистрация: 12.10.2013
Сообщений: 2,837
19.02.2017, 14:37 #11
Цитата Сообщение от AndreyL Посмотреть сообщение
ускорение раза в три
Exe в 100 раз?
0
AndreyL
1 / 1 / 1
Регистрация: 11.09.2013
Сообщений: 50
20.02.2017, 06:16  [ТС] #12
Цитата Сообщение от Centurio Посмотреть сообщение
Вот вам код, который считает гораздо быстрее цикла
Если позволите, задам еще вопрос. Попытался аналогичным способом векторизовать цикл расчета средних и дисперсий - не получилось. Не подскажете, можно ли векторизовать такой цикл? Сначала экспонируем матрицу L и приводим к единице каждую строку.
Matlab M
1
2
3
4
5
6
7
maxll = max(L,[],2);
%minus maxll to avoid underflow
post = exp(bsxfun(@minus, L, maxll));
%density(i) is \sum_j \alpha_j P(x_i| \theta_j)/ exp(maxll(i))
density = sum(post,2);
%normalize posteriors
L = bsxfun(@rdivide, post, density);
Это что бы был понятен смысл самого цикла - L после этой операции содержит статистические веса, т.е. вероятность принадлежности строки матрицы X распределению с центром и дисперсией в строках матриц Y и V соответственно (и таких распределений k штук). Цикл расчета взвешенных средних и дисперсий такой
Matlab M
1
2
3
4
5
6
7
pSum1=1./sum(L,1);
for i=1:k;
  Y(i,:)=L(:,i)'*X*pSum1(i);
  xProm=bsxfun(@minus,X,Y(i,:));
  xProm=bsxfun(@times, xProm.*xProm, L(:,i));
  V(i,:)=sum(xProm,1)*pSum1(i);
end;
Сложность в том, что если развернуть матрицы и они будут содержать по n*k строк, то суммы придется считать не по всему столбцу, а по кусочкам длиной n.
0
Centurio
Модератор
724 / 653 / 187
Регистрация: 13.09.2015
Сообщений: 2,355
20.02.2017, 07:32 #13
AndreyL, у вас матрица L имеет размер (n,k), а матрица Х - размер (n,m). k и m разные, поэтому строка из матрицы L не может быть поэлементно умножена на строку из матрицы Х, т.к. размеры не совпадают. Быть может, каждую строку матрицы Х нужно умножить на один элемент из соответствующей строки матрицы L и таких умножений сделать k раз, по количеству элементов в строках матрицы L?
0
AndreyL
1 / 1 / 1
Регистрация: 11.09.2013
Сообщений: 50
20.02.2017, 07:40  [ТС] #14
Совершенно верно, если Вы имеете ввиду строку 3 цикла, т.е. расчет взвешенных средних. Но матрица L(:,i)' имеет размер 1*n, и умножение на Х дает вектор 1*m.
0
Centurio
Модератор
724 / 653 / 187
Регистрация: 13.09.2015
Сообщений: 2,355
20.02.2017, 08:54 #15
Цитата Сообщение от AndreyL Посмотреть сообщение
Но матрица L(:,i)' имеет размер 1*n, и умножение на Х дает вектор 1*m.
Почему же вам, в таком случае, не умножить сразу транспонированную матрицу L на матрицу Х?
0
AndreyL
1 / 1 / 1
Регистрация: 11.09.2013
Сообщений: 50
21.02.2017, 14:34  [ТС] #16
Цитата Сообщение от Centurio Посмотреть сообщение
Почему же вам, в таком случае, не умножить сразу транспонированную матрицу L на матрицу Х?
Для среднего, наверное, так
Matlab M
1
Y=L'*X.*repmat(pSum1,m,1)';
А как быть со взвешенными дисперсиями?

Добавлено через 25 минут
Да, я же обновился до 2016b, со средним можно и так
Matlab M
1
Y=L'*X.*pSum1';
Добавлено через 18 часов 5 минут
Посмотрел оригинал - функция gmcluster.m в стандартном MATLAB 2016b. Там тоже делается циклом
Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
S.PComponents = sum(post,1);
S.Sigma = zeros(1,d,'like',X);
for j = 1:k
    if S.PComponents(j) == 0 
        %When the small posterior probablities are set to zero,
        % it's possilble to get a cluster with zero prior.
        continue;
    end
    S.mu(j,:) = post(:,j)' * X / S.PComponents(j);
    Xcentered =  X - S.mu(j,:);
    S.Sigma = S.Sigma + post(:,j)' *(Xcentered.^2);
end
S.Sigma = S.Sigma/sum(S.PComponents)+regVal;
S.mu у меня называется Y, S.Sigma у меня называется V, post - это L после экспонирования и приведения к 1. Я только не понял 11 строку - неужели действительно проще сложить два вектора целиком, чем одному элементу вектора присвоить значение? regVal по дефолту равен 0.
0
21.02.2017, 14:34
MoreAnswers
Эксперт
37091 / 29110 / 5898
Регистрация: 17.06.2006
Сообщений: 43,301
21.02.2017, 14:34

Как ускорить работу кода, избавившись от Memo?
Всем привет! Есть у меня парсер данных со спорт сайта, так вот в нем я паршу...

какой программой можно произвести побитово операцию обратную XOR?
Есть программа WinHex которая позволяет работать с информацией побитово. ...

Можно ли сделать так, чтобы во время цикла можно было вводить какую нибудь символ?
Можно ли сделать так, чтобы во время цикла можно было вводить какую нибудь...


Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
16
Ответ Создать тему
Опции темы

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2018, vBulletin Solutions, Inc.
Рейтинг@Mail.ru