Язык 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.