Векторизация | Авторская платформа Pandia.ru

Партнерка на США и Канаду по недвижимости, выплаты в крипто

30% recurring commission
Выплаты в USDT
Вывод каждую неделю
Комиссия до 5 лет за каждого referral

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

Индексация. Индекс — это тоже массив. Выбор отдельно стоящих элементов вектора: A([3,5]).

Поменять местами строки:

M = magic(5)

M([1,5],:) = M([5,1],:)

Можно ли применить тот же прием к матрице. Линейный индекс.

sub2ind — переход от индексации вида Aстрока, столбец (нижние индексы — subscripts) к линейной.

idx = sub2ind(size, I,J)

Возвращает линейный индекс idx элементов, расположенных в строках I и столбцах J для матрицы размера size.

Обратный переход осуществляется функцией ind2sub.

Переход к линейному индексу можно реализовать и самому. Для матриц:

idx = i + (j-1)*size(A,1);

i, j — номера строки и столбца заданного элемента.

Логическая индексация. Логический индекс — массив того же размера, что и исходный, состоящий из 1 и 0. Единицы стоят на местах элементов исходного массива, удовлетворяющих логическому условию, нули — на местах элементов, не удовлетворяющих этому условию.

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

Составить массив из положительных элементов заданного массива D:

D = [-0.2 1.0 1.5 3.0 -1.0 4.2 3.14];

Выделим элементы D, отвечающие условию D >= 0:

P = D(D>=0);

Заменим элементы M, большие 4, на обратные. Поэлементный вариант:

НЕ нашли? Не то? Что вы ищете?

for i=1:3

for j=1:3

if M(i, j) > 4

M(i, j) = - M(i, j);

end

Векторизованный с помощью логической индексации код:

M(M > 4) =-M(M > 4);

Вернемся к вычислению «пи» из прошлой лекции. Код с циклами

n = 1e6;

x = rand(1,n);

fun = sqrt(1-x.^2);

y = rand(1,n);

ns = 0; % счетчик попаданий внутрь окружности

for i=1:n

if y(i) <= fun(i)

ns = ns + 1;

end

Выделенный цветом фона код теперь можно заменить одной строкой

ns = sum(y <= fun);

Можно также снабдить эти вычисления красивым графиком

plot(x, fun,'.',x(y <= fun),y(y <= fun),'g.',x(y > fun),y(y > fun),'r.'),axis('square')

Оценка скорости выполнения кода: tic, toc. Не забываем очищать память!

Кусочно-непрерывные функции. Функция sign(x) для скалярного аргумента x:

if x < 0

y = -1;

else

y = 1;

end

Но если x является массивом, то необходимо будет проверить истинность логического условия для каждого элемента массива. Воспользуемся логической индексацией.

Для x = [1 -3 0 4 3 -5 -7] условие x >= 0 дает

ans = [1 0 1 1 1 0 0]

а условие x < 0:

ans = [0 1 0 0 0 1 1]

Объединяя, эти выражения, получим

y = (x >= 0) - (x < 0)

что дает в результате

y = [1 -1 1 1 1 -1 -1]

Функции all и any. all(A) — возвращает логическую 1, если ни один элемент массива A не является нулем. Например, для A = [0.53 0.67 0.01 0.38 0.07 0.42 0.69] вызов all(A < 0.05) даст 0, а вызов all(A>0.001) — 1.

Благодаря этому all удобно использовать при задании различных условий, например

if all(A < 1)

выполнять команды

end

any(A) — возвращает логическую 1, если хотя бы один элемент массива A не является нулем. Например, any(A < 0.05) даст 1.

Задание массивов. Как создать матрицу

1 1 1

2 2 2

3 3 3

4 4 4

5 5 5

чтобы при этом не использовать циклов?

Можно воспользоваться функцией repmat

v = (1:5)';

A = repmat(v,1,3)

Другой вариант — матричное умножение

v = (1:5)';

A = v*ones(1,3)

Трюк Тони (Tony’s Trick).

v = (1:5)';

A = v(:,ones(1,3));

В последней строке можно поставить ones(3,1), так как массив используется как линейный индекс. ones(3,3) дает 9 столбцов.

Сравните самостоятельно скорости работы предложенных вариантов, используя tic, toc.

Вычтем среднее значение из каждого столбца матрицы М

M = rand(10,5);

V = mean(M);

Трудность здесь в том, что из двумерного массива вычитается массив одномерный. Это можно векторизовать так

for i=1:5

M(:,i)=M(:,i)-V(i);

end

Но лучшим решением будет

M=M-V(ones(10,1),:);

Кумулятивные сумма и произведение. Кумулятивной суммой (cumsum) данной последовательности является последовательность частичных сумм элементов данной последовательности. Например, кумулятивные суммы последовательности [a, b, c, d, …] равны [a, a+b, a+b+c, a+b+c+d, …].

cumsum удобно использовать для векторизации циклов, где для накапливания используется суммирование. Так, код

n=1000;

x(1)=100;

for i=2:n

x(i) = x(i-1)+3;

end

дает в результате

100

100 + 3

100 + 3+ 3

…

что, пользуясь cumsum, можно преобразовать как

n=1000;

x=cumsum([100 3*ones(1,n-1)]);

Например, отвечая на вопрос «какая сумма окажется на счете через заданное число периодов?», мы составляли программу с циклом

si = 100; % исходный вклад

p = 0.2; % норма простого процента

s = si; % сумма на счету

n = 6;

for i = 1:n

s = s + si*p;

end

Вместо цикла можно использовать кумулятивную сумму

s2=cumsum([si p*si*ones(1,n)]);

s2(end)

(хотя более простым выходом в этом случае будет воспользоваться формулой для вычисления процента).

Аналогично реализуется кумулятивное произведение (cumprod). Кумулятивные произведения последовательности [a, b, c, d, …] равны [a, ab, abc, abcd, …]. cumprod удобно использовать для векторизации кода цикла с умножением. Например

n=6;

x=100*ones(1,n);

p=rand(1,n-1);

for i = 2:n

x(i) = x(i-1)*(1+p(i-1));

end

можно записать так

n=6;

x=cumprod([100 1+p]);

Вычисление определенных интегралов. Вычислим :