Векторизация: как ускорить алгоритм машинного обучения с помощью x78

Учитывая уравнение, мы увидим, как шаг за шагом можно добиться более эффективного кода не только в 78 раз с точки зрения скорости, но и используя всего 3 строки кода! Давайте погрузимся в это ...

Вступление

Как интерпретируемый язык Python for циклы по своей сути медленнее, чем их аналог на C. Это большое узкое место для печально известного языка программирования, поскольку алгоритмы глубокого и машинного обучения сильно зависят от матричных операций, которые выполняются через for циклы.

Вот почему разработчики разработали пакеты, такие как numpy, которые предлагают векторизованные действия с numpy массивами. Это означает, что он сдвигает цикл for, который обычно выполняется в Python, до уровня C, что намного быстрее.

Python + скорость уровня C = рай

Эта проблема

(вы можете пропустить части объяснения, если вы согласны с алгоритмом EM)

Мы хотим использовать алгоритм Expectation-Maximization (EM) для задачи обучения без учителя (скажем, например, распознавание рукописных цифр в наборе данных MNIST), и наши данные являются двоичными (например, двоичные изображения). Естественный способ - смоделировать наши данные как модель смеси Бернулли; взвешенная сумма распределений Бернулли, где каждое распределение имеет свой собственный скалярный вес π и собственный средний вектор μ и представляет собой кластер данных (например, если наши данные представляют собой изображения цифр 2, 3 и 4 и мы используем 3 Бернулли для их моделирования, один Бернулли будет цифрами 2, другой - цифрами 4 и т. д.). В целом это делает первое вектором, а второе - матрицей.

Пусть N = количество наблюдений, D = размерность одного наблюдения и K = количество кластеров. Поскольку это важно для нашей проблемы, мы имеем следующие типы случайных величин:
• X; наши данные представляют собой матрицу NxD
(N количество изображений, D их размерность → 5 изображений размером 28 * 28 образуют матрицу X 5x784)
• π; вектор K, скаляр для каждого распределения, представляющего его вес.
(например, три вектора Бернулли могут иметь π = [0,2, 0,75, 0,05] взвешенный вектор)
• μ; матрица среднего KxD для каждого кластера.
(одно изображение имеет размерность D = 28 * 28 = 784, где каждое из них представляет значение пикселя. Принимая среднее значение для каждого пикселя изображений, принадлежащих одному кластеру скажем цифра 2, мы получаем средний вектор 784. Таким образом, μ будет матрицей KxD)

На шаге E нас особенно интересует ожидаемое значение апостериорной части скрытых переменных, или так называемые обязанности.

γ фактически возвращает ожидаемое значение наблюдения (изображения) n, принадлежащего кластеру k.
γ - это матрица NxK; для каждого наблюдения мы назначаем вероятность принадлежности к каждому кластеру. Тот, у которого есть максимум, - это тот, который мы назначаем.

Зачем я все это рассказываю?

«Самое важное в векторизации - понять размерность переменных».

Расчет ответственности - это тот, который мы собираемся векторизовать.

Подводя итог:
• X: матрица NxD
• π: вектор 1xK
• μ: матрица KxD
• γ: матрица NxK

Трубопровод

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

observations = [5, 10, 20, 50, 100, 200, 500, 1000]
for n in observations:
    X_test = bin_train_data[:n]
    D_test, K_test = X_test.shape[1], 10
    mu_test = np.random.uniform(low=.25, high=.75, 
                                size=(K_test,D_test))
    pi_test = np.ones(K_test) / K_test
    t0 = time.time()
    gamma_test = E_step_1(X_test, mu_test, pi_test)
    runtime = time.time() - t0
    assert gamma_test.shape == (n, K_test)

Не стесняйтесь сначала попробовать его сами!

Попытка №1

С первой попытки мы напишем все с помощью циклов for; при операциях с векторами / матрицами - только скаляры.

Посмотрев на уравнения, мы можем увидеть, что есть 3 цикла; По одному для каждого примера N, по одному для каждого кластера K и по одному для каждого измерения каждого объекта D, и мы собираемся выполнить цикл в этом порядке. Итак, мы собираемся заполнять матрицу γ по одному элементу за раз.

def E_step(X, mu, pi):
    N, D = X.shape
    K = pi.shape[0]
    gamma = np.zeros((N, K))
    for n in range(N):
        for k in range(K):
            m = 1
            for i in range(D):
                m *= mu[k][i]**X[n][i] * (1-mu[k][i])**(1-X[n][i])
            gamma[n][k] = m * pi[k]
        gamma[n] /= gamma[n].sum()
    return gamma

И наши результаты видны на графике ниже.

Мы определенно можем добиться большего!

Попытка №2

Хорошо начинать с внутреннего цикла и продвигаться к внешнему. И это именно то, что мы собираемся делать!

Мы хотим избавиться от цикла for D. Таким образом, каждый член, зависящий от D, должен теперь стать вектором. Внутри этого цикла for у нас есть две переменные; μ и x (см. Уравнение(2)). Таким образом, x и μ → векторы. Проблема; Это μ**x, вектор в степени другого вектора, и его трудно вычислить. Если бы мы только могли обойти это ...

Есть функция, которая превращает степенную операцию в умножение. Правильно, это логарифм! Итак, давайте логарифмируем наше выражение, а затем возьмем показатель степени, если результат!

Предпочтительны операции над логарифмическими вероятностями, поскольку они обеспечивают численную стабильность.

Несмотря на то, что в нашем случае это никак не повлияет, каждый раз, когда вы используете журналы, используйте константу epsilon внутри выражения для стабильности (чтобы не обнуляться, с is -inf).

Таким образом, нам предстоит поэлементное умножение векторов. Легкий ;)

def E_step(X, mu, pi):
    N, D = X.shape
    K = pi.shape[0]
    gamma = np.zeros((N, K))
    for n in range(N):
        for k in range(K):
            log_gamma = np.log(pi[k]) + (X[n] * np.log(mu[k]) \
                        + (1 - X[n])*np.log(1 - mu[k])).sum()
            gamma[n][k] = np.exp(log_gamma)
        gamma[n] /= gamma[n].sum()
    return gamma

И наши результаты…

Это ОГРОМНАЯ победа! Похоже, что ось x конкурирует с Algor. 1! Но мы можем сделать лучше;)

Попытка №3

По одной петле за раз: K ход!

В процессе векторизации мы движемся следующим образом:

скаляр → вектор → матрица

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

Мы берем нашу предыдущую реализацию и хотим удалить цикл K for. Таким образом, каждый скаляр, который зависит от K, превратится в vector, а каждый vector в matrix. Это означает, что X останется прежним, а μ превратится в матрицу, а π и γ в векторы. Обратите внимание на последний; как γ поле строка за строкой, результат нашего выражения теперь должен быть вектором! Таким образом, операция μ и X должна приводить к 1xK вектору, и быстрые индикаторы: (i) они должны суммироваться с вектором π, который также 1xK (ii) результатом является строка матрицы γ, также 1xK вектор.

Кодирование приводит к:

def E_step(X, mu, pi):
    N, D = X.shape
    K = pi.shape[0]
    gamma = np.zeros((N, K))
    for n in range(N):
        log_gamma = np.log(pi) + np.log(mu) @ X[n] \
                    + np.log(1 - mu) @ (1 - X[n])
        gamma[n] = np.exp(log_gamma)
        gamma[n] /= gamma[n].sum()
    return gamma

И вот результаты:

Удивительный! Нам удалось вдвое сократить время на n=1000! Нет никакого сравнения с Algor. 1. Но можем ли мы сделать лучше?

Попытка №4

У нас еще осталась одна петля. Можем ли мы иметь расчет без петель и питона? N, твое время вышло!

Поскольку мы собираемся преобразовать matrix * vector операций в matrix @ matrix операцию, нам нужно взять транспортную матрицу первой (@ - это обычное умножение матриц). Имейте в виду, что теперь на выходе должна быть вся матрица γ. Думаю, теперь вы уже поняли, как это происходит;).

Итак, наш код будет

def E_step(X, mu, pi):
    gamma = np.exp(np.log(pi) + X @ np.log(mu.T) \
            + (1 - X) @ np.log(1 - mu.T))    
    gamma /= gamma.sum(axis=1)[:, np.newaxis]
    return gamma

Ни одной петли! Код выглядит элегантно и состоит всего из трех строк! А теперь результаты, барабанная дробь ...

И все, от этого не может быть лучше! Для n=1000 мы перешли от среды выполнения 11.688 → 0.012, используя всего три строки кода!

Резюме

Итак, что вам нужно сделать, если вы хотите векторизовать выражение:

• Разберитесь с размерами матриц.
• Ручка и бумага: запишите формулу, перейдите от суммирования к суммированию и превратите ее в эквивалентную матричную операцию < br /> • Математика - ваш друг; всегда рассуждайте о количестве измерений, которые должно возвращать любое выражение; следите за операциями суммирования соседей, так как они имеют одинаковые размеры
• Цикл за циклом, шаг за шагом: скаляр → вектор → матрица
• Возьмите логарифм и убедитесь, что введена постоянная нормализации epsilon
• Закодируйте свою векторизованную версию вашего подхода и сияйте: D