C++ и Python

Язык Python

Вычисления с библиотекой numpy

Библиотека numpy является основой экосистемы python для научных вычислений. В основе numpy лежит код на языке C, обеспечивающий высокую производительность. Центральным объектом в numpy является многомерный массив numpy.ndarray. Разнообразие и удобство инструментов для работы с многомерными массивами позволяют очень быстро реализовывать разнообразные вычислительные процедуры с минимальным количеством кода и во многих случаях без использования циклов.

Создание и свойства объектов numpy.ndarray

Массив numpy можно создать из списка чисел:

import numpy as np
a = np.array([1, 2, 3])
a.ndim   # 1
a.shape  # (3,)
a.size   # 3
a.dtype  # np.int64
a[0]     # 1

Двумерный массив из списка списков:

import numpy as np
b = np.array([
    [1, 2, 3.],
    [4, 5, 6]
])
b.ndim   # 2
b.shape  # (3, 3)
b.size   # 6
b.dtype  # np.float64
b[0]     # [1., 2., 3.]
b[0, 0]  # 1.
b[:, 0]  # [1., 4.]

Количество измерений может быть любым. В numpy есть несколько функций для создания массивов:

# трехмерный массив с типом элементов np.float64, заполненный нулями
zeros = np.zeros(shape=(3, 3, 3), dtype=np.float64)
# двумерный массив целых чисел, заполненный единицами
ones = np.ones((3, 3), dtype=np.int64)
# не проинициализированный четырехмерный массив комплексных чисел
empty = np.empty((4, 4, 4, 4), dtype=np.complex128)
# единичная матрица 3x3
unit3 = np.eye(3)
# диагональная матрица 4x4
diag4 = np.diag([1, 2, 3, 4])
# равномерная сетка в диапазоне значений
np.linspace(0, 1, 101)  # [0.00, 0.01, ..., 1.00]
# равномерная логарифмическая сетка
np.logspace(1, 9, 9)  # [1.e+01, 1.e+02, ..., 1.e+09]
# матрица 4x4, заполненная последовательными целыми числами
mtx = np.arange(16).reshape(-1, 4)

В последнем примере мы использовали метод reshape, позволяющий изменить форму массива. Для успешного выполнения этой операции массив должен содержать подходящее количество элементов. В качестве первого аргумента мы передали -1. Это значит, что первая размерность будет вычислена из размера массива и значения второй размерности: 16 // 4 = 4.

Метод reshape не изменяет объект, у которого он вызван. Кроме того, если это возможно, метод reshape не выполняет копирование данных. Вместо этого он возвращает новое представление (view) массива. Такое поведение очень удобно при работе с большими массивами, но надо иметь его в виду при изменении объекта через разные представления:

a = np.arange(6).reshape(-1, 2)
# [[0, 1],
#  [2, 3],
#  [4, 5]]
b = a.reshape(-1, 3)
# [[0, 1, 2],
#  [3, 4, 5]]
b[0, 0] = 9
# [[9, 1, 2],
#  [3, 4, 5]]
print(a)
# [[9, 1],
#  [2, 3],
#  [4, 5]]

Метод ravel позволяет получить одномерное представления массива любой формы. Как и reshape, ravel не выполняет копирование, если это возможно. Действие ravel аналогично вызову reshape(-1).

Операции с массивами, broadcasting

С массивами numpy можно выполнять арифметические операции. Они при этом будут выполняться поэлементно:

a = np.arange(9).reshape(-1, 3)
# [[0, 1, 2],
#  [3, 4, 5],
#  [6, 7, 8]]
a + 1
# [[1, 2, 3],
#  [4, 5, 6],
#  [7, 8, 9]]
a * 2
# [[ 0,  2,  4],
#  [ 6,  8, 10],
#  [12, 14, 16]]
a**2
# [[ 0,  1,  4],
#  [ 9, 16, 25],
#  [36, 49, 64]]
2**a
# [[  1,   2,   4],
#  [  8,  16,  32],
#  [ 64, 128, 256]]
a < 5
# [[ True,  True,  True],
#  [ True,  True, False],
#  [False, False, False]]

Сравним для иллюстрации скорость применения операций к элементам массива numpy и к элементам list с помощью модуля timeit:

a = list(range(10000))
b = np.arange(10000)

list(map(lambda x: x**2, a))  # Способ 1
# 3.19 ms ± 50.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

[x**2 for x in a]  # Способ 2
# 2.74 ms ± 23.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

b**2  # Способ 3
# 5.33 µs ± 215 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Кроме более краткой записи, мы получили ускорение более чем в 500 раз.

В numpy определено множество функций, которые могут работать с массивами любой формы, например:

alpha = 2.*np.pi * np.linspace(0, 1, 12)
sin_alpha = np.sin(alpha)

Вот еще примеры арифметических операций с массивами:

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])
a + b  # [5, 7, 9]
a * b  # [ 4, 10, 18]
a**b   # [  1,  32, 729]
a.reshape(-1, 1) + b.reshape(1, -1)
# [[5, 6, 7],
#  [6, 7, 8],
#  [7, 8, 9]]

Последний пример требует пояснения. Поэлементные операции с массивами имеют очевидный смысл, когда формы массивов совпадают. В последнем примере мы сложили разные массивы: вектор-строку и вектор-столбец. Если форма массивов не совпадает, то numpy пытается применить процедуру broadcasting — копирование массива вдоль новой оси так, чтобы новая форма позволяла выполнять операцию. В данном случае строка была скопирована вдоль столбцов, а столбец — вдоль строк:

# [[1, 2, 3],     [[4, 4, 4],    [[5, 6, 7],
#  [1, 2, 3],  +   [5, 5, 5],  =  [6, 7, 8],
#  [1, 2, 3]]      [6, 6, 6]]     [7, 8, 9]]

Broadcasting во многих случаях упрощает работу с массивами, обеспечивая интуитивное обобщение операций на массивы разной формы.

Способы индексации

Массивы numpy поддерживают все приемы индексации python и предоставляют дополнительные возможности. Во-первых, индексы для нескольких измерений можно передавать через запятую:

a = np.arange(9).reshape(-1, 3)
# [[0, 1, 2],
#  [3, 4, 5],
#  [6, 7, 8]]
a[:-1,:-1]
# [[0, 1],
#  [3, 4]]

Одним из самых полезных приемов индексирования является использование булевых масок:

a = np.arange(9).reshape(-1, 3)
# [[0, 1, 2],
#  [3, 4, 5],
#  [6, 7, 8]]
a[a > 5]  # [6, 7, 8]

Логические операторы применяются поэлементно, и результат их применения может быть использован в качестве булевой маски: будут выбраны только элементы, для которых маска содержит значение True. Форма массива-маски должна соответствовать форме индексируемого массива.

Элементы массива можно выбирать с помощью списка индексов:

a = np.arange(9)
a[[5, 0, 1, 5]]  # [5, 0, 1, 5]

Такая индексация полезна при сортировке массива. Выполнить обычную сортировку можно с помощью функции np.sort. Функция np.argsort вместо сортировки массива возвращает перестановку индексов, которую нужно применить для получения отсортированного массива. Индексация по списку индексов позволяет отсортировать массив по заданной перестановке:

a = np.arange(9)
a[np.argsort(a % 3)]  # сортируем массив по остатку деления на три
# [0, 3, 6, 1, 4, 7, 2, 5, 8]

Случайные числа в numpy

Объект Generator из модуля numpy.random предоставляет интерфейс ко многим дискретным и непрерывным случайным распределениям. Создавать объект Generator рекомендуется следующим образом:

import numpy as np
rnd = np.random.default_rng(seed=None)
type(rnd)  # <class 'numpy.random._generator.Generator'>

Для данного значения параметра seed, например seed=5, будет получаться одна и та же последовательность случайных чисел:

rnd1 = np.random.default_rng(seed=1)
rnd1.random(5)
# [0.51182162, 0.9504637 , 0.14415961, 0.94864945, 0.31183145]
rnd2 = np.random.default_rng(seed=2)
rnd2.random(5)
# [0.26161213, 0.29849114, 0.81422574, 0.09191594, 0.60010053]
rnd3 = np.random.default_rng(seed=1)
rnd3.random(5)
# [0.51182162, 0.9504637 , 0.14415961, 0.94864945, 0.31183145]

Это удобно для воспроизводимости численных экспериментов. Значение по умолчанию seed=None означает, что генератор будет проинициализирован случайно. В последнем примере использована функция random, которая возвращает выборку из равномерного распределения в интервале [0.0, 1.0). Вот еще несколько примеров:

rnd = np.random.default_rng()
rnd.random(4, 2)  # массив формы (4, 2)
# выборка целых чисел от 0 до 10 включительно
rnd.integers(0, 10, (5, 5))  # массив формы (5, 5)
# выборка из нормального распределения
rnd.normal(loc=5, scale=2, size=(4, 6))  # массив формы (4, 6)

Функции Generator.shuffle и Generator.permutation позволяют получать случайные перестановки элементов массива. Функция Generator.shuffle модифицирует переданный массив, а функция Generator.permutation возвращает новый массив:

arr = np.arange(10)         # arr = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
rnd.shuffle(arr)            # arr = [8, 3, 5, 0, 4, 7, 6, 1, 2, 9]
per = rnd.permutation(arr)  # per = [3, 7, 1, 0, 4, 5, 9, 6, 8, 2]

Функция Generator.choice позволяет получать случайные выборки из заданного набора эдементов:

rnd.choice([0, 1], 10)  # [1, 1, 0, 0, 1, 1, 1, 1, 0, 1]

Линейная алгебра в numpy

Пакет numpy поддерживает операции линейной алгебры. В частности, произведение векторов и матриц по правилам линейной алгебры может быть выполнено с помощью оператора @ или функции numpy.dot:

a = np.array([1, 2, 3]).reshape(1, -1)  # вектор-строка
b = np.array([4, 5, 6]).reshape(1, -1)  # вектор-строка
mtx = np.arange(9).reshape(-1, 3)

Скалярное произведение:

a @ b.T  # (1, 3) x (3, 1) -> (1, 1)
# [[32]]

Для получения вектора-столбца мы выполнили транспонирование вектора b.

Внешнее произведение:

a.T @ b  # (3, 1) x (1, 3) -> (3, 3)
# [[ 4,  5,  6],
#  [ 8, 10, 12],
#  [12, 15, 18]]

Умножение матрицы на вектор, умножение матриц происходит аналогично:

mtx @ a.T  # (3, 3) x (3, 1) -> (3, 1)
# [[ 8],
#  [26],
#  [44]]

a @ mtx  # (1, 3) x (3, 3) -> (1, 3)
# [[24, 30, 36]]

mtx @ mtx  # (3, 3) x (3, 3) -> (3, 3)
# [[ 15,  18,  21],
#  [ 42,  54,  66],
#  [ 69,  90, 111]]

Модуль numpy.linalg содержит различные алгоритмы линейной алгебры, например:

  • numpy.linalg.inv возвращает обратную матрицу
  • numpy.linalg.det возвращает определитель матрицы
  • numpy.linalg.eig возвращает собственные значения и собственные векторы матрицы
  • numpy.linalg.svd выполняет сингулярное разложение матрицы

Полный набор доступных инструментов можно найти в документации.

Свертка тензоров

Операции умножения массивов по правилам линейной алгебры становятся гораздо менее очевидными, когда приходится иметь дело с размерностью больше, чем два. В этом случае мы имеем дело с тензорами и их сверткой по различным индексам. В таких случаях мы рекомендуем рассмотреть использование функции np.einsum. Предлагаем читателю самостоятельно изучить детали ее использования, приведем только несколько примеров:

def mtxab(A, B):
    """ A и B — трехмерные массивы с одинаковой первой размерностью.
        Интерпретируем A и B как списки матриц и выполняем матричное
        произведение каждой пары матриц
    """
    return np.einsum('kij, kjl -> kil', A, B)
def mtxabt(A, B):
    """ Аналогично mtxab, но дополнительно выполняем транспонирование
        матриц из массива B
    """
    return np.einsum('kij, klj -> kil', A, B)
def mtxabc(A, B, C):
    """ Выполняем произведение соответствующих трех матриц из
        трехмерных массивов A, B и C
    """
    return np.einsum('kij, kjl, klm -> kim', A, B, C)

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

Сериализация и десериализация массивов numpy

Зачастую вычисления занимают много времени. В таких случаях результат вычислений разумно сохранять в файл и при необходимости читать данные из файла вместо того, чтобы повторять вычисления. Функция np.save позволяет сохранять массив в бинарный файл с расширением .npy. Функция np.savetxt сохраняет массив в текстовом формате. Функция np.savez позволяет сохранить несколько массивов в один файл. Функции load, loadtxt выполняют чтение сохраненных массивов. Рассмотрим простой пример:

a = np.array([1, 2, 3])
np.save('array', a)

b = np.load('array.npy')  # b = [1, 2, 3]

Резюме

В этом разделе мы обсудили основы работы с массивами numpy, инструменты для генерирования случайных чисел с numpy, средства numpy для линейной алгебры, средства сериализации массивов numpy.

Источники