Математика в библиотеке Numpy

Все курсы > Программирование на Питоне > Занятие 10

Библиотека Numpy — это, прежде всего, инструмент для проведения математических и статистических вычислений над массивами чисел. Рассмотрим этот функционал Numpy более подробно.

Откроем ноутбук к этому занятию⧉

Основные операции

По умолчанию математические операции в массиве Numpy выполняются поэлементно. Возьмем две матрицы 2 x 3.

Сложение и вычитание массивов

Рассмотрим сложение и вычитание массивов. Вначале выполним поэлементное сложение.

Мы также можем выполнить сложение двух элементов внешнего измерения одного массива.

Аналогично выполняется вычитание массивов.

Умножение и деление

Рассмотрим поэлементное умножение. Такая операция еще называется произведением Адамара (Hadamard product).

Аналогично выполняется операция деления.

Размерность массива должна позволять выполнить поэлементную операцию.

При этом как и в случае со сложением и вычитанием, Numpy позволяет умножать и делить на число.

Другие математические операции

Помимо базовых операций в Numpy есть множество других математических функций. Например, мы можем найти максимальное значение массива с помощью функции np.max().

Как и в случае np.sum(), мы можем указать ось, вдоль которой выполняется операция.

Большинство функций имеют соответствующий метод, который приводит к точно такому же результату.

Выбор между функцией и методом во многом определяется удобством прочтения кода. Например, метод транспонирования матрицы .T более похож на привычную математическую запись этой операции нежели функция np.transpose().

Дополнительно такая операция демонстрирует возможность поэлементного сравнения двух матриц. Про эту и другие логические операции мы поговорим в следующем разделе.

При этом некоторые функции отличаются от соответствующих им методов. В частности, как мы помним, функция np.sort(), в отличие от метода .sort(), не меняет исходный массив.

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

Логические операции

Сравнение массивов

Мы можем поэлементно проверить равенство этих массивов.

Кроме того, мы можем проверить, меньше ли элементы массива a элементов массива b.

Эту же задачу можно решить с помощью специальной функции np.less().

Мы можем сравнить элементы массива с числом.

Про эту операцию мы уже говорили в контексте логической маски массива Numpy, а также при рассмотрении порогового преобразования изображения.

Логические операторы

Логические операции И, ИЛИ, исключающее ИЛИ и отрицание (логическое НЕ) реализованы с помощью функций np.logical_and(), np.logical_or(), np.logical_xor() и np.logical_not() соответственно.

Для того чтобы понять как они работают, полезно вспомнить таблицы истинности (truth table) для этих операций. Начнем с логического И. Операция истинна, только если оба значения истинны (оба True или 1).

Применим функцию np.logical_and() и выведем True только если значение массива a больше нуля И меньше трех.

Рассмотрим логическое ИЛИ. Здесь хотя бы одно значение должно быть истинно.

Функция np.logical_or() выведет значение True либо если элемент массива b меньше трех, либо если элемент равен единице.

Логическая операция исключающего ИЛИ выдаст True только если одно из значений истинно.

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

Теперь применим функцию np.logical_xor(). Выведем True, если элемент массива a больше одного или (но не одновременно) элемент массива b меньше трех.

Логическое НЕ переводит истинное значение в ложное и ложное в истинное.

Применим функцию np.logical_not() к массиву a.

Обсудим подробнее, как мы получили такой результат. В Питоне все числа кроме нуля считаются истинными, ноль наоборот приравнивается к ложному значению. Убедимся в этом с помощью двух несложных примеров.

Как следствие, так как в массиве a ни одно из чисел не равно нулю, все они считаются истинными. Применив функцию np.logical_not(), все значения стали ложными.

Функции np.all() и np.any()

Создадим массив, состоящий из логических значений.

Функция np.all() позволяет проверить, все ли элементы массива оцениваются как истинные.

Функция np.any() наоборот проверяет истинно ли хотя бы одно значение.

Полный перечень логических операций также можно найти в документации Numpy⧉.

Функция np.where()

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

Наша задача заключается в том, чтобы заменить двойку на «Не сдал», а остальные оценки на «Сдал». Вначале создадим массив из логических значений True или False.

Теперь создадим списки с категориями, такой же длины, как и исходный массив scores.

Посмотрим на один из списков.

Передадим функции np.where() массив с логическими значениями и два списка.

Рассмотрим, как мы получили такой результат. Функция np.where() взяла логический массив, и там, где было значение True, заменила это значение на элемент из первого массива, а там где False — из второго.

Если не передавать списки с категориями, можно узнать, например, индексы двоечников ИЛИ отличников.

Обратите внимание, что в данном случае мы не можем использовать оператор or, потому что нам нужно применить условие к каждому элементу массива (побитовая операция). Возможным вариантом было бы использование функции np.logical_or().

Такой результат можно использовать в качестве фильтра исходного массива.

Константы в Numpy

Питон позволяет использовать несколько математических констант, в частности, число Пи и число Эйлера.

Векторизация и трансляция

Векторизация

Как мы увидели выше, при выполнении арифметических операций одновременно над всеми элементами массива мы не использовали цикла for. Такую запись принято называть векторизацией (vectorization). Одним из преимуществ является скорость исполнения кода.

Сравним время исполнения кода с циклами и без них на простом примере. Вначале создадим два длинных одномерных массива.

Используем цикл for и list comprehension для поэлементного умножения двух векторов.

В данном случае фактическое время исполнения кода на серверах Google (wall time) составило 355 миллисекунд (при повторном исполнении кода время будет немного отличаться).

Отдельно поговорим про так назывыемые магические команды (magic commands) %time и %%time в Google Colab.

• Команда %time позволяет замерить время исполнения отдельной строки кода.
• При этом, %%time с двумя символами процента замеряет время исполнения всей ячейки.

Теперь решим ту же задачу, но уже с помощью векторизованного кода.

Как мы видим, время исполнения сократилось более чем в 100 раз.

Отмечу, что циклы в данном случае все же применяются, но написаны они (как и почти весь Numpy) на более быстром с точки зрения исполнения кода языке C.

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

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

Лаконичный и эффективный код принято называть питоническим (pythonic).

Трансляция

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

В этом случае, для того чтобы сделать возможной векторизацию сложения, Numpy предварительно увеличивает меньший массив (число) до размеров большего (вектора) и затем выполняет сложение без использования цикла. Такой подход называется трансляцией (broadcasting).

Прежде чем выполнить векторизованную операцию (сложение), размерность числа увеличивается до размерности вектора.

Напоследок замечу, что далеко не всегда можно транслировать массивы таким образом, чтобы сделать их совместимыми для определенной операции. Существуют правила трансляции (broadcasting rules). Их рассмотрение однако находится за рамками этого занятия.

Умножение векторов, матриц и тензоров

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

Как мы узнали на двенадцатом занятии вводного курса, скалярное произведение (dot product) удобно использовать, например, в моделях машинного обучения для умножения вектора данных (признаков) на вектор весов.

Вначале повторим самый простой случай скалярного произведения двух векторов.

В библиотеке Numpy такая операция выполняется с помощью функции np.dot().

Уточнение терминов

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

В русском языке под скалярным произведением понимается поэлементное умножение и сложение получившихся произведений. Вероятно, такое название связано с тем, что результатом операции будет число (скаляр). При этом если мы просто хотим провести поэлементное (покомпонентное) умножение, речь идет произведении Адамара. Также возможно умножение вектора или матрицы на число.

В английском языке скалярное произведение принято называть dot product или scalar product, поэлементное умножение — element-wise product или Hadamard product, умножение на число — scalar multiplication.

Для наглядности соберем результат этих наблюдений в таблицу.

Произведение матрицы и вектора

Пример из области машинного обучения

В случае если наши данные содержат несколько наблюдений (строк), речь уже идет об умножении двумерного массива или матрицы на вектор весов.

В этом случае для каждой строки наблюдения мы получаем конкретный числовой результат (прогноз).

Давайте пошагово разберем, как были сделаны такие вычисления. На первом этапе умножаем первую строку на вектор-столбец весов и получаем 14 (как и в предыдущем примере).

Затем умножаем вторую строку матрицы на вектор-столбец весов и получаем число 32.

И наконец остается умножить третью строку на вектор весов, чтобы получить 50.

Обратите внимание, в получившемся векторе столько же строк, сколько было в матрице данных (то есть три), и столько же столбцов, сколько было в векторе весов (то есть один).

Это первое свойство скалярного произведения выполняется всегда для векторов, матриц и тензоров. Если вы вернетесь к умножению двух векторов, то увидите, что это правило также соблюдалось. Одна строка вектора данных и один столбец вектора весов в результате дают скаляр (число).

Кроме того, количество столцов матрицы совпадает с количеством строк вектор-столбца. В противном случае умножение было бы невозможно.

Это второе свойство скалярного произведения. Мы еще раз рассмотрим эти свойства, когда будем умножать матрицу на матрицу.

Реализуем умножение матрицы и вектора на Питоне.

Что интересно, если переставить матрицу и вектор местами и перемножить по тем же правилам (строку первого множителя на столбец второго), то результат будет иным.

Из этого можно сделать два важных вывода:

Вывод 1. Скалярное произведение матриц некоммутативно. Другими словами, если в скалярном произведении, хотя бы один из множителей — матрица, то от перемены мест множителей произведение меняется.

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

Пример для системы линейных уравнений

Помимо моделей машинного обучения, операцию умножения матрицы на вектор удобно использовать, например, при перемножении коэффициентов и неизвестных системы линейных уравнений (set of simultaneous equations).

$$ \begin{cases} 4x_1-2x_2+3x_3=1 \\ x_1+3x_2-4x_3=-7 \\ 3x_1+x_2+2x_3=5 \end{cases} $$

Перепишем эту систему как произведение матрицы и вектора.

$$ \begin{pmatrix} 4 & -2 & 3 \\ 1 & 3 & -4 \\ 3 & 1 & 2 \end{pmatrix} \cdot \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 1 \\ -7 \\ 5 \end{pmatrix} $$

Если обозначить матрицу через A, а векторы через X и B, то систему линейных уравнений можно записать следующим образом.

$$ AX = B $$

Мы займемся решением этой системы уравнений сразу после того, как познакомимся с понятием обратной матрицы.

Умножение двух матриц

Перейдем к более общему случаю умножения матрицы на матрицу. Для наглядности возьмем две матрицы размерностью 2 x 3 и 3 x 2 соответственно. Создадим первую матрицу.

Создадим вторую.

Пошаговый разбор

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

Шаг 1. Перемножим компоненты первой строки первой матрицы и первого столбца второй матрицы и сложим произведения. Получается число пять. Запишем этот результат в первую строку и первый столбец новой, результирующей матрицы.

Шаг 2. Теперь снова используем первую строку первой матрицы, при этом из второй матрицы возьмем уже второй столбец. Перемножим компоненты и сложим произведения. Результат запишем во второй столбец первой строки.

Шаги 3 и 4. Перемножим вторую строку первой матрицы сначала с первым, затем со вторым столбцом второй. Так мы получим результат для второй строки новой матрицы.

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

Выполним ту же операцию на Питоне.

Опять же обратите внимание, если множители поменять местами, произведение изменится.

Свойства умножения матриц

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

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

Другими словами, если у нас есть матрица A размерностью i x j и матрица B размерностью j x k, то результатом скалярного произведения станет матрица C размерностью i x k.

$$ A_{i \times j} \cdot B_{j \times k} = C_{i \times k}$$

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

Второе свойство легко продемонстрировать на примере. Попробуем умножить, например, две матрицы размерностью 2 x 3 каждая.

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

Как мы видим, Питон говорит о том, что второе изменение (столбцы) первой матрицы (dim 1) отличается по количеству элементов от первого (строки) измерения (dim 0) второй.

Умножение матриц с помощью циклов for

В учебных целях продемонстрируем умножение матриц a и b с помощью циклов for. Для этого обозначим оси через индексы i, j, k и будем умножать вдоль индекса j.

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

На втором этапе запишем результаты в нулевую матрицу 2 x 2, воспользовавшись соответствующими индексами.

Обратите внимание на три момента:

• Мы легко смогли заменить оси массива на индексы
• Умножение происходило вдоль «внутреннего» индекса j
• Индекс j исчез из результата вычислений

Эти наблюдения помогут нам в дальнейшем понять, что такое соглашение о суммировании Эйнштейна и как работает функия np.einsum().

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

Умножение матриц в нейронных сетях

Когда мы изучали основы нейронных сетей в рамках вводного курса, то при написании кода прямого распространения (forward propagation) использовали словарь для хранения данных и цикл for для умножения данных на соответствующие веса (weight) и прибавления смещений (bias).

Теперь мы можем решить эту же задачу гораздо эффективнее с помощью массива Numpy и скалярного произведения.

Архитектура нейронной сети

Возьмем несложную нейросеть, которая на входе будет принимать рост человека (X), пропускать эти данные через один скрытый слой (h₁, h₂) и на выходе предсказывать вес (o₁) и обхват шеи (o₂). Приведем схему нейросети.

Прямое распространение

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

Шаг 1. Представим, что у нас три наблюдения, обозначим их как x_a, x_b и x_c и поместим в матрицу размерностью 3 x 1. Умножим эти значения на матрицу весов w₁ и w₂ с размерностью 1 x 2. В результате получаем матрицу размерностью 3 x 2.

Шаг 2. Прибавляем к этой матрице смещения b₁ и b₂. Причем обратите внимание, смещение b₁ по правилам трансляции прибавляется к первому столбцу, а b₂ — ко второму.

Шаг 3. Результат сложения мы пропускаем через новую для нас функцию активации ReLU (о ней я расскажу чуть ниже).

С помощью выделенных красным ячеек можно проследить вычисления для одного компонента матрицы.

Теперь рассмотрим второй этап от скрытого до выходного слоя.

Шаг 4. Вначале умножаем матрицу скрытого слоя (3 x 2) на матрицу весов w₃, w₄, w₅ и w₆ (2 x 2). Получается матрица 3 x 2.

Шаг 5. Затем мы прибавляем смещения b₃ и b₄.

Обратите внимание, при такой архитектуре сети можно было бы подумать, что в скрытом слое нейрон h₁ отвечает за вес, а h₂ за обхват шеи. Однако если присмотреться, оба этих нейрона учавствуют в расчете как веса тела, так и обхвата шеи. При этом веса модели и смещения относятся строго каждый к «своему» показателю.

Шаг 6. Передаем результат в функцию активации для получения конечных значений веса тела (o₁) и обхвата шеи (o₂) для каждого из наблюдений.

Функция активации ReLU

Теперь несколько слов про функцию активации ReLU или Rectified Linear Unit. Ее принцип очень прост: если передаваемое ей значение больше нуля, то мы оставляем значение без изменений, в остальных случаях возвращем ноль.

$$ f(x) = max(0, x) $$

На Питоне можно использовать функцию np.maximum(), которой в качестве первого параметра мы передадим ноль.

Графически ReLU выглядит следующим образом.

Помимо встроенной в Numpy функции np.maximum(), ReLU можно реализовать самостоятельно, используя свойства трансляции.

В данном случае с помощью (x > 0) мы формируем массив из единиц (если x больше 0) и нулей (в остальных случаях). После этого мы поэлементно умножаем исходный массив на ноль или один. Умножение на единицу не изменяет значения, умножение на ноль дает ноль.

Реализация на Питоне

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

Теперь пропишем веса и смещения.

Объявим функцию прямого распространения.

Вызовем фунцию, передав ей матрицу наблюдений.

Подбор весов и смещений, то есть обратное распространение (back propagation), мы рассмотрим на курсе по оптимизации.

Тензорное произведение

Мы конечно помним, что массивы Numpy не обязательно должны быть двумерными, при этом функция np.dot() используется для вычисления скалярного произведения именно двумерных массивов.

Для того чтобы перемножить тензоры любой размерности существует функция np.tensordot(). Рассмотрим применение этой функции вначале на двумерных, затем на трехмерных массивах.

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

Функция np.tensordot() на примере двумерного тензора

Вновь возьмем матрицы а и b.

Повторим умножение матриц с помощью функции np.dot().

По умолчанию, функция np.dot() выполняет умножение по последней оси первого множителя (матрицы a) и по предпоследней оси второго множителя (матрицы b).

Идея np.tensordot() заключается в том, чтобы явно указать оси массива, по которым будет происходить покомпонентное умножение и сложение получившихся произведений (по-английски такое сокращение еще называют sum-reduction). Делается это с помощью параметра axes.

Повторим умножение матриц a и b, но уже через функцию np.tensordot().

В данном случае мы явно прописали, что при умножении используем последнюю ось (ось 1) матрицы a и предпоследнюю ось (ось 0) матрицы b.

Обратите внимание, при умножении исчезают (сокращаются) те оси, вдоль которых происходит операция. Например, здесь мы взяли матрицы размерностью 3 x 2 и 2 x 3, умножение происходило по «внутренним» осям (в которых по два элемента) и в результате мы получили матрицу 3 x 3 (то есть остались только «внешние» оси).

Теперь давайте для закрепления материала попробуем поиграть с параметром axes и посмотрим какие результаты мы получим.

Рассмотрим вот такой пример.

На схеме это выглядело бы следующим образом.

Как вы видите, мы транспонировали матрицу b, но при этом выполнили умножение вдоль ее последней оси (axis = 1). Правило сокращения измерений также действует. Мы умножили две матрицы 2 x 3 и 2 x 3 по вторым изменениям (в которых по три элемента) и у нас осталась матрица 2 x 2.

Теперь умножим а на b, но перевернем оси. В матрице a будем умножать по первой (предпоследней) оси, а в матрице b по второй (последней).

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

Например, проведем умножение по обеим осям обеих матриц.

По сути такая операция аналогична «вытягиванию» первого массива вдоль оси 1, второго массива вдоль оси 0 и скалярному произведению двух получившихся векторов.

Теперь рассмотрим особые случаи, когда в параметр axes мы передаем значения 1, 2 или 3.

Особые случаи axes = 0, axes = 1 и axes = 2

Параметр axes = 0 позволяет умножить каждый элемент первой матрицы, на элементы второй.

Размерность при этом всегда увеличивается, чтобы «уместить» увеличившееся количество матриц.

Параметр axes = 1 позволяет провести умножение по последней оси первой матрицы и первой оси второй. Для двумерных массивов речь идет об осях 1 и 0 соответственно.

Параметр axes = 2 для двумерных массивов выполняет функцию умножения вдоль обеих осей обоих массивов.

Функция np.tensordot() на примере трехмерного тензора

Перейдем к трехмерным тензорам.

Посмотрим на умножение вдоль последней оси первого массива и вдоль предпоследней оси второго.

Что интересно, именно так работает функция np.dot(), если ее применить к массивам, в которых больше двух измерений.

Особые случаи для трехмерного тензора

Рассмотрим особые случаи для трехмерного массива. Их будет четыре: axes = 0, axes = 1, axes = 2 и axes = 3.

Результат работы функции np.tensordot() c параметром axes = 0 аналогичен двумерному массиву. Каждый элемент первого тензора поочередно умножается на каждый элемент второго.

Параметр axes = 1 призводит умножение по последней оси первого тензора и по первой оси второго.

Мы конечно можем прописать оси вручную и получим точно такой же результат.

С параметром axes = 2 мы производим умножение по первым двум осям первого тензора и по последним двум осям второго.

И наконец параметр axes = 3 выполняет умножение вдоль всех трех осей тензора.

Соглашение Эйнштейна и функция np.einsum()

Соглашение о суммировании Эйнштейна (Einstein Summation Convention) упрощает запись операций над тензорами через их индексы (вместо осей функции np.tensordot()).

Сам Альберт Эйнштейн в письме другу говорил следующее: «Я сделал важное открытие в математике. Мне удалось отказаться от знака суммы в тех случаях, когда суммирование происходит по индексу, который появляется дважды» ^[1].

Рассмотрим это утверждение на примере умножения двух матриц.

Принцип соглашения о суммировании

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

$$ \sum^n_{j=1} a_{i, j} \cdot b_{j, k} = c_{i, k} $$

В данном случае мы берем две матрицы и умножаем их по индексу j (как мы уже делали выше), который появляется дважды.

Без знака суммирования эту запись можно представить следующим образом:

$$ a_{ij} b_{jk} = c_{ik}$$

Функция np.einsum() отражает такую запись.

Функция np.einsum()

Вначале мы передаем в np.einsum() исходные и желаемые индексы ('ij, jk -> ik'), а затем сами матрицы a и b.

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

Как вы помните, np.dot() выдавала аналогичный результат.

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

Транспонирование, сложение и поэлементное умножение в np.einsum()

Транспонировать матрицу можно просто поменяв индексы строк и столбцов местами.

Кроме того, мы можем сложить элементы вдоль определенного индекса. Для этого мы берем оба индекса матрицы и оставляем только тот, по которому хотим провести суммирование.

Например, если мы хотим сложить вдоль строк (то есть оси 1), то оставим индекс i.

И наоборот, если нужно просуммировать столбцы, оставим j.

Кроме того, мы можем просуммировать элементы по обоим индексам.

При поэлементном умножении размерность должна совпадать. Для этого мы можем в одной операции транспонировать матрицу b и затем, за счет использования одинаковых индексов i и j для обеих матриц провести поэлементное умножение.

Операции с векторами в функции np.einsum()

Посмотрим как использовать функцию np.einsum() при работе с одномерными массивами.

Найдем сумму элементов.

Выполним поэлементное умножение.

Найдем скалярное произведение.

Типы и свойства матриц

Создание специальных матриц

Изучим несколько особых видов матриц и способы их создания с помощью библиотеки Numpy.

Нулевая матрица

С точки зрения математики, нулевая матрица (zero matrix) — это матрица размера m x n, в которой все элементы равны нулю. С точки зрения Numpy — это двумерный массив, заполненный нулями.

$$ O_{3, 2} = \begin{pmatrix} 0 & 0 \\ 0 & 0 \\ 0 & 0 \end{pmatrix} $$

Его можно создать с помощью функции np.zeros().

Единичная матрица

Единичная матрица (identity matrix) представляет собой квадратную матрицу (square matrix) размера n x n, в которой по диагонали слева направо (ее называют главной диагональю, main diagonal) расположены единицы, а остальные элементы заполнены нулями.

$$ I = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} $$

В библиотеке Numpy такую матрицу можно создать с помощью функций np.identity() и np.eye().

Отличие заключается в том, что np.eye() позволяет сместить диагональ с помощью параметра k.

В линейной алгебре нулевая и единичная матрица выполняют ту же функцию, что и 0 и 1 в обычных вычислениях. В частности, произведение любой матрицы на нулевую матрицу равно нулевой матрице.

$$ A \cdot O = O $$

При умножении квадратной матрицы на единичную матрицу такого же размера получится исходная матрица. При этом обратите внимание, такая операция коммутативна.

$$ A_{m, n} \cdot I_n = I_m \cdot A_{m, n} = A_{m, n} $$

Диагональная матрица

Главная диагональ матрицы не обязательно должна состоять из единиц. В этом случае говорят про диагональную матрицу (diagonal matrix).

Диагональ также можно заполнить произвольными значениями.

Функции np.diagonal() или np.diag() возвращают значения диагонали.

При этом функция np.diag() может не только возвращать значения диагонали, но и использоваться для создания диагональной матрицы.

Ленточная матрица

Ленточная матрица (band или banded matrix) — это матрица, в которой ненулевыми являются также диагонали, примыкающие к главной.

$$ B = \begin{pmatrix} 3 & 1 & 0 \\ 2 & 2 & 2 \\ 0 & 1 & 1 \end{pmatrix} $$

В Питоне нет специальной функции для создания ленточной матрицы. Как поступить? Мы можем вначале создать диагональную матрицу.

Затем аналогичным способом создать еще две матрицы такого же размера со смещенными (offset) вниз и вверх диагоналями. Для этого можно использовать параметр k функции np.diag().

Останется лишь поэлементно сложить получившиеся матрицы.

Можно и так.

Более подробно про генерацию случайных чисел и, в частности, функцию np.random.randint() мы поговорим на следующем занятии.

Треугольная матрица

Треугольная матрица (triangular matrix) — это матрица, в которой все элементы ниже или выше главной диагонали равны нулю. Матрицу называют верхней треугольной матрицей (upper triangular), если все элементы ниже главной диагонали равны нулю и нижней треугольной (lower triangular), если выше.

Верхнюю треугольную матрицу можно построить из обычной матрицы с помощью функции np.triu().

Параметр k этой функции позволяет регулировать, ниже какой диагонали обнулять элементы.

Нижняя треугольная матрица создается с помощью функции np.tril().

Диагональ нижней треугольной матрицы также можно сдвигать вверх и вниз с помощью положительных и отрицательных значений параметра k.

Матрица Numpy

Помимо массива (ndarray) в библиотеке Numpy существует еще один объект, матрица Numpy (Numpy matrix). Рассмотрим отличия матрицы от массива:

• В матрице может быть только два измерения
• Если перемножить массивы с помощью символа звездочки *, мы выполним поэлементное умножение (произведение Адамара), при перемножении матриц — этот символ предполагает скалярное произведение.
• Символ ** возводит каждый элемент массива в соответствующую степень, в случае объекта matrix предполагается возведение в степень всей матрицы (другими словами ее скалярное умножение саму на себя заданное количество раз).

На курсе я буду стараться избегать применения matrix, поскольку разработчики Numpy рекомендуют отказаться от использования⧉ этого объекта. Вероятно это связано с тем, что весь функционал объекта matrix присутствует в объекте ndarray (массиве Numpy).

Транспонирование матрицы

Еще раз вернемся к транспонированной матрице.

Транспонированная матрица (the transpose of a matrix) — это матрица, в которой строки и столбцы исходной матрицы поменялись местами.

Как вы уже конечно знаете, в библиотеке Numpy для транспонирования матрицы можно использовать функцию np.transpose() или метод .T.

Соответственно, если размерность исходной матрицы была m x n, то у транспонированной матрицы размерность будет n x m.

Продемонстрируем с помощью Numpy некоторые свойства траспонированных матриц.

Дважды транспонированная матрица равна исходной.

$$ (A^T)^T = A $$

Транспонированная сумма двух матриц равна сумме транспонированных матриц.

$$ (A + B)^T = A^T + B^T $$

Транспонированное произведение двух матриц равно произведению транспонированных матриц, поставленных в обратном порядке.

$$ (AB)^T = B^TA^T $$

Нахождение обратной матрицы

Обратная матрица (inverse of a matrix) — это такая матрица A⁻¹, при умножении которой на матрицу А, получится единичная матрица I.

$$ AA^{-1} = A^{-1}A = I $$

Обратную матрицу можно сравнить с обратным числом. Например, для числа 5 обратным числом будет 1/5.

$$ 5 \times \frac{1}{5} = 1 $$

Для нахождения обратной матрицы можно использовать функцию np.linalg.inv() из модуля линейной алгебры библиотеки Numpy. Продемонстрируем использование этой функции на очень простом примере.

Приведем несколько полезных свойств и определений.

• Только квадратная матрица может иметь обратную матрицу
• Если квадратная матрица обратима, то говорят, что это невырожденная матрица (invertible или non-singular)
• Если квадратная матрица не обратима, то говорят, что это вырожденная или сингулярная матрица (non-invertible или singular)

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

Решение систем линейных уравнений

Теперь вернемся к системе линейных уравнений, которую, напомню, мы выразили с помощью следующего уравнения.

$$ AX = B $$

Оказывается, что если матрица A обратима, то мы можем умножить обе части уравнения на A⁻¹.

$$ A^{-1}AX = A^{-1}B $$

Это возможно, так как произведение матрицы на обратную ей матрицу дает единичную матрицу, а произведение матрицы на единичную матрицу равно этой матрице.

$$ A^{-1}A = I $$

$$ XI = X $$

Исходя из этого, мы можем выразить X через произведение A⁻¹ на B.

$$ X = A^{-1}B $$

Применим этот метод на практике и решим рассмотренную ранее систему уравнений.

$$ \begin{pmatrix} 4 & -2 & 3 \\ 1 & 3 & -4 \\ 3 & 1 & 2 \end{pmatrix} \cdot \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 1 \\ -7 \\ 5 \end{pmatrix} $$

Вначале создадим массивы Numpy.

Теперь у нас есть два способа для решения это задачи.

Способ 1. Воспользуемся функцией np.linalg.inv() для нахождения обратной матрицы и методом .dot() для умножения матриц.

Корнями уравнения будут x₁ = −1, x₂ = 2, x₃ = 3.

Способ 2. Используем функцию np.linalg.solve(), которой передадим матрицы A и B.

Проверим правильность найденного решения.

В каких ситуациях матрица A обратима, а у системы уравнений есть решение?

• Неизвестных должно быть столько же, сколько уравнений. Другими словами, матрица A должна быть квадратной.
• Одно уравнение нельзя превратить в другое через операции сложения уравнений и их умножения на число. В этих случаях говорят, что строки матрицы (уравнения) линейно независимы (linearly independent rows).

Статистика в Numpy

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

Меры среднего

Для измерения среднего значения (measures of central tendency) можно использовать функции np.mean(), np.average() и np.median().

Функция np.mean()

Функция np.mean() вычисляет среднее арифметическое значение данных вдоль заданной оси.

Функция np.median()

С помощью функции np.median() можно вычислить медиану.

Аналогичного результата можно добиться, применив функцию np.percentile() и передав ей параметр 50, ведь медиана — это 50-й процентиль (что такое процентили, мы узнаем на курсе по анализу и обработке данных).

Функция np.average()

Функция np.average() позволяет найти средневзвешенное значение (weighted average). При вычислении средневзвешенного значения (W) каждое наблюдение (X_i) умножается на соответствующий вес (w_i) и получившиеся произведения складываются.

$$ W = \sum_{n}^{i=1} w_i X_i $$

Рассмотрим на примере. Предположим, что студент получил оценки 3, 4 и 5 на экзаменах по линейной алгебре, математическому анализу и истории соответственно. При этом вес линейной алгебры в средней оценке за семестр составляет 50%, вес матана составляет 40%, а вес — истории только 10%.

Каким будет средневзвешенный результат за семестр?

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

Теперь в качестве упражнения выполним эти вычисления вручную.

Меры разброса

Давайте посмотрим, как с помощью Numpy можно посчитать некоторые уже знакомые нам меры разброса (measures of dispersion).

Функция np.ptp()

Функция np.ptp() (от английского peak-to-peak) находит разницу между максимальным и минимальным значениями (т.е. диапазон, range). Мы уже использовали ее, когда изучали основы нейронных сетей.

Функции np.std() и np.var()

Функция np.std() считает среднее квадратическое отклонение (СКО, standard deviation).

Одновременно функция np.var() считает дисперсию (т.е. квадрат СКО, variance).

Полный перечень статистических функций можно найти в документации Numpy⧉.

Более подробно с теорией и практикой описательной статистики мы опять же познакомимся на следующем курсе по анализу и обработке данных.

Подведем итог

На сегодняшнем занятии мы подробно рассмотрели математику в библиотеке Numpy.

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

Последующие же курсы анализа данных, оптимизации, линейной алгебры и статистики вывода будут в большей степени сосредоточены именно на изучении математических объектов.

Вопросы для закрепления

Вопрос. На примере умножения двух векторов объясните, чем произведение Адамара отличается от скалярного произведения?

Вопрос. Что такое векторизация кода?

Вопрос. В чем отличие функции np.dot() от np.tensordot()?

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

Перед этим рекомендую обратить внимание на дополнительные упражнения⧉.