Matrices#
Табличный способ хранения данных чрезвычайно широко распространён, и если эти данные числовые, то с математической точки зрения они образуют матрицу — двумерный массив из чисел. Для эффективной обработки и анализа таких данных, а также обучения моделей машинного обучения на них, важно хорошо ориентироваться в матричных операциях и разложениях и их свойствах.
Матрицы принято обозначать заглавными жирными буквами: \(\boldsymbol A\), \(\boldsymbol B\), \(\boldsymbol C\), \(\boldsymbol X\), \(\boldsymbol Y\) и т.п. Множество матриц, имеющих \(m\) строк и \(n\) столбцов обозначается через \(\mathbb R^{m\times n}\). Элемент \(i\)-й строки и \(j\)-го столбца матрицы \(\boldsymbol A \in \mathbb R^{m\times n}\) будем обозначать \(A_{ij}\). Допустимо также обозначение \(a_{ij}\), особенно когда дело доходит до подробной записи
Эту же матрицу можно представить как набор строк или столбцов:
Если \(m=n\), то матрица \(\boldsymbol A\) называется квадратной, в противном случае — прямоугольной. Вектор \(\boldsymbol x \in \mathbb R^n\) можно считать матрицей размера \(n\times 1\) (столбец) или \(1\times n\) (строка).
Транспонированная матрица \(\boldsymbol A^\mathsf{T} \in \mathbb R^{n\times m}\) в качестве строк содержит столбцы матрицы \(\boldsymbol A\): \(A^\mathsf{T}_{ij} = A_{ji}\). В частности, вектор-столбец получается транспонированием из вектора-строки, и наоборот, чем и объясняется обозначение \(\boldsymbol x^\mathsf{T}\) для векторов-строк.
Feature matrix
In machine learning a numeric dataset is usually represented by the feature matrix \(\boldsymbol X\) in which each training sample is represented by a row vector \(\boldsymbol x_i^\mathsf{T}\). Hence,
np.array
can handle multidimenstional arrays including matrices. For example, let’s generate a random matrix of shape \(2\times 3\)
import numpy as np
A = np.random.rand(2, 3)
A
array([[0.4495285 , 0.3634043 , 0.67298861],
[0.67521062, 0.10220722, 0.69235149]])
For transposing use attribute .T
:
A.T
array([[0.89698755, 0.52984442],
[0.59051524, 0.76303108],
[0.22936288, 0.33085304]])
Square matrices#
A square matrix \(\boldsymbol A \in \mathbb R^{n\times n}\) has the same number of rows and columns. Each square matrix has two diagonals: main (primaty, principle) and secondary (side, antidiagonal).
Trace of a matrix is equal to the sum of the elements on its main diagonal: \(\mathrm{tr}(\boldsymbol A) = \sum\limits_{i=1}^n a_{ii}\).
Identity matrix#
All ones on the main diagonal, all zeros outside it. Identity matrix is denoted as \(\boldsymbol I\) or \(\boldsymbol I_n\) if we want to underscore its shape:
Identity matrix in NumPy is created by np.eye
:
I = np.eye(4, dtype=np.int32)
I
array([[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]], dtype=int32)
By default the underlying type is float:
np.eye(3)
array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])
Diagonal matrix#
Can have nonzero elements only on the main diagonal:
Identity matrix is a special case of a diagonal matrix:
np.diag
can either create a diagonal matrix or extract the main diagonal of a matrix. Here we create a diagonal matrix:
D = np.diag([1, 2, 3])
D
array([[1, 0, 0],
[0, 2, 0],
[0, 0, 3]])
A = np.random.normal(size=(3, 3))
A
array([[-0.44717254, 1.4078295 , -0.6438964 ],
[ 0.30502598, 1.0711746 , 1.15591706],
[ 1.51407551, 1.83549146, -1.7051814 ]])
And here the extraction of the main diagonal happens:
np.diag(A)
array([-0.44717254, 1.0711746 , -1.7051814 ])
Triangular matrices#
An upper triangular matrix \(\boldsymbol U\) has all zeros above the main diagonal: \(U_{ij} = 0\) if \(i > j\),
A lower triangular matrix \(\boldsymbol L\), on the other hand, has all zeros below the main diagonal: \(L_{ij} = 0\) if \(i < j\),
Note that \(\boldsymbol U\) is upper triangular iff \(\boldsymbol U^\mathsf{T}\) is lower triangular.
Symmetric matrices#
A square matrix \(\boldsymbol A\) is called
symmetric if \(\boldsymbol A^\mathsf{T} = \boldsymbol A\);
skew-symmetric if \(\boldsymbol A^\mathsf{T} = -\boldsymbol A\).
Any square matrix \(\boldsymbol A\) can be represented as a sum of a symmetric and a skew-symmetric matrix:
Orthogonal matrices#
A square matrix \(\boldsymbol Q = [\boldsymbol q_1 \ldots \boldsymbol q_n]\) is orthogonal if its columns are orthonormal, i.e., \(\langle \boldsymbol q_i, \boldsymbol q_j\rangle = \delta_{ij}\), where \(\delta_{ij}\) is Kronecker delta.
An example of orthogonal matrix is given by the rotation matrix
Matrix operations#
Перечислим поэлементные операции с матрицами. Пусть матрицы \(\boldsymbol A, \boldsymbol B \in \mathbb R^{m\times n}\), \(\alpha \in \mathbb R\), \(1 \leqslant i \leqslant m\), \(1\leqslant j \leqslant n\).
Сложение матриц:
Умножение матрицы на скаляр:
Поэлементное умножение (произведение Адамара) и деление матриц:
\[ \boldsymbol C = \boldsymbol A \odot \boldsymbol B \iff C_{ij} = A_{ij} B_{ij},\quad \boldsymbol D = \boldsymbol A \oslash \boldsymbol B \iff D_{ij} = \frac{A_{ij}}{B_{ij}} \](последнее возможно, разумеется, только если все элементы матрицы \(\boldsymbol B\) не равны нулю). Стоит отметить, что под произведением матриц обычно понимают совсем иную операцию, нежели поэлементное умножение, и о ней предстоит отдельный разговор.
Первые два свойства указывают на то, что матрицы можно рассматривать как вектора особого вида. Можно даже пойти дальше и записать матрицу размера \(m\times n\) в виде одномерного массива длины \(mn\) (например, построчно), и получится вектор в чистом виде. Такая операция применяется, к примеру, в свёрточных нейронных сетях, когда на определённом этапе представленная матрицей или тензором картинка выпрямляется (flatten) в вектор для дальнейшего прохождения через один или несколько полносвязных слоёв.
Matrix operations in NumPy are quite natural:
A = np.random.randint(-2, 3, size=(2, 3))
B = np.arange(1, 7).reshape(2, 3)
print(A)
print(B)
[[ 0 -2 0]
[ 1 1 -1]]
[[1 2 3]
[4 5 6]]
A + B
array([[1, 0, 3],
[5, 6, 5]])
A - B
array([[-1, -4, -3],
[-3, -4, -7]])
A * B
array([[ 0, -4, 0],
[ 4, 5, -6]])
A / B
array([[ 0. , -1. , 0. ],
[ 0.25 , 0.2 , -0.16666667]])
Matrix-vector product#
Matrix times column = column#
If \(\boldsymbol A \in \mathbb R^{m\times n}\), \(\boldsymbol x \in \mathbb R^n\), then \(\boldsymbol A\boldsymbol x = \boldsymbol b \in \mathbb R^m\), and its coordinates are
Observe that \(b_i = \boldsymbol a_i^\mathsf{T} \boldsymbol x\) where \(\boldsymbol a_i^\mathsf{T}\) — \(i\)-th row of \(\boldsymbol A\). Alternatively, if we denote \(j\)-th column of \(\boldsymbol A\) as \(\boldsymbol a_j\), then
Hence, \(\boldsymbol b = \boldsymbol A\boldsymbol x\) is a linear combination of columns of \(\boldsymbol A\).
Since matrix times vector operation is quite similar to inner products, the Python syntax is the same as here:
A = np.arange(10).reshape(2, -1)
b = np.array([2, 1, -1, 0, -2])
A
array([[0, 1, 2, 3, 4],
[5, 6, 7, 8, 9]])
A.dot(b), np.dot(A, b), A @ b
(array([-9, -9]), array([-9, -9]), array([-9, -9]))
If the number of columns of \(\boldsymbol A\) is different form the size of \(\boldsymbol b\), the try to multiply them will result in an error:
A @ np.ones(4)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[21], line 1
----> 1 A @ np.ones(4)
ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 4 is different from 5)
Row times matrix = row#
Умножая матрицу на вектор-столбец, мы получаем снова вектор-столбец; умножение матрицы на вектор-строку записывают в обратном порядке:
Результат такого умножения представляет собой линейную комбинацию строк матрицы \(\boldsymbol A\): \(\boldsymbol c^\top = \sum\limits_{i=1}^m y_i \boldsymbol a_i^\top\).
A = np.arange(-2, 4).reshape(3, 2)
b = np.array([2, 1, -1])
A
array([[-2, -1],
[ 0, 1],
[ 2, 3]])
b.T @ A
array([-6, -4])
Without transposition \(\boldsymbol b\) is treated as a column vector, and attempt to multiply it by \(\boldsymbol A\) fails:
b @ A
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[22], line 1
----> 1 b @ A
ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 2 is different from 5)
Row times column#
Заметим, что матрицу \(\boldsymbol A\) размера \(1 \times n\) можно представить в виде вектора-строки \(\boldsymbol a^\mathsf{T}\). С другой стороны, согласно определению результат её умножения на вектор \(\boldsymbol x \in \mathbb R^n\) равен числу
а это в точности скалярное произведение \(\langle \boldsymbol a, \boldsymbol x\rangle\). Отсюда и проистекает его весьма популярное обозначение через матрично-векторное произведение \(\boldsymbol a^\mathsf{T} \boldsymbol x\), которым мы также будем активно пользоваться в дальнейшем.
Matrix product#
Матрицу \(\boldsymbol A\) можно умножить на матрицу \(\boldsymbol B\), если количество столбцов матрицы \(\boldsymbol A\) равно числу строк матрицы \(\boldsymbol B\). А именно, произведением матриц \(\boldsymbol A \in\mathbb R^{m\times n}\) и \(\boldsymbol B \in\mathbb R^{n\times p}\) является матрица \(\boldsymbol C = \boldsymbol A \boldsymbol B \in\mathbb R^{m\times p}\), каждый элемент которой равен
Элементы матрицы \(\boldsymbol C\) можно записать в виде \(C_{ik} = \boldsymbol a_i^\mathsf{T} \boldsymbol b_k\), где \(\boldsymbol a_i^\mathsf{T}\) — строки матрицы \(\boldsymbol A\), а \(\boldsymbol b_k\) — столбцы матрицы \(\boldsymbol B\). Вычисление каждого элемента \(C_{ik}\) требует \(O(n)\) арифметических операций, поэтому общая сложность подсчёта произведения матриц составляет \(O(mnp)\). Если все матрицы квадратные размера \(n\times n\), то сложность будет \(O(n^3)\).
А быстрее можно?
Алгоритм Штрассена позволяет перемножать квадратные матрицы чуть быстрее, а именно за \(O(n^{\log_2 7})\), что примерно составляет \(O(n^{2.81})\). Алгоритм Штрассена был разработан ещё в 1969 году, и с тех пор было предложено множество алгоритмов с лучшей асимптотикой. Подробнее про эту гонку можно почитать в этой статье, где текущая SoTA указана как \(O(n^{2.3728596})\). Математики не оставляют надежд добиться асимптотики \(O(n^2)\) или хотя бы \(O(n^{2+\varepsilon})\) для любого сколь угодно малого \(\varepsilon > 0\), однако, все эти изыскания носят сугубо теоретический характер из-за гигантских констант в O-большом и сложности реализации подобных алгоритмов.
Впрочем, сегодня новые алгоритмы перемножения матриц придумывают не только математики, но и нейронные сети. Например, компания DeepMind с помощью глубокого обучения подкрепления разработала алгоритм AlphaTensor способный перемножать матрицы размера \(5\times 5\) эффективнее, чем алгоритм Штрассена (\(76\) умножений против \(80\)).
Произведение матриц \(\boldsymbol A \boldsymbol B\) можно представить также в виде набора столбцов, получающихся умножением матрицы \(\boldsymbol A\) на столбцы матрицы \(\boldsymbol B\): если \(\boldsymbol B = [\boldsymbol b_1 \ldots \boldsymbol b_p]\), то
Или же можно строки матрицы \(\boldsymbol A\) умножать на матрицу \(\boldsymbol B\): если
Matrix product in NumPy is similar to that of vectors:
A = np.array([[1, 2, 3], [4, 5, 6]])
B = np.arange(12).reshape(3, 4)
print(A.shape, B.shape)
(2, 3) (3, 4)
A.dot(B)
array([[ 32, 38, 44, 50],
[ 68, 83, 98, 113]])
Properties of matrix product#
\((\boldsymbol{AB})\boldsymbol{C} = \boldsymbol{A}(\boldsymbol{BC})\) (ассоциативность)
\((\boldsymbol A + \boldsymbol B)\boldsymbol{C} = \boldsymbol{AC} + \boldsymbol{BC}\), \(\boldsymbol{C}(\boldsymbol A + \boldsymbol B) = \boldsymbol{CA} + \boldsymbol{CB}\) (дистрибутивность)
\((\boldsymbol{AB})^\mathsf{T} = \boldsymbol B^\mathsf{T} \boldsymbol A^\mathsf{T}\)
\(\mathrm{tr}(\boldsymbol{AB}) = \mathrm{tr}(\boldsymbol{BA})\)
Последнее свойство обобщается до равенства \(\mathrm{tr}(\boldsymbol{ABC}) = \mathrm{tr}(\boldsymbol{CAB})\), которое справедливо для любого числа сомножителей с согласованными размерами. Это так называемое циклическое свойство следа матрицы.
Let’s check this properties for random matrices. Start with associativity:
m, n, p, q = 2, 3, 4, 5
A = np.random.randn(m, n)
B = np.random.randn(n, p)
C = np.random.randn(p, q)
np.allclose((A @ B).dot(C), A.dot(B @ C))
True
Distributivity:
np.allclose((A + B) @ C, A @ C + B @ C)
True
Transposing:
A = np.random.randn(m, n)
B = np.random.randn(n, p)
np.allclose((A @ B).T, B.T @ A.T)
True
Trace:
A = np.random.randn(m, n)
B = np.random.randn(n, m)
np.trace(A @ B), np.trace(B @ A)
(0.6484650575676159, 0.6484650575676162)
Несколько дополнительных свойств для квадратных матриц:
\(\boldsymbol{AI} = \boldsymbol{IA} = \boldsymbol A\);
если \(\boldsymbol P\) — матрица перестановки, то матрица \(\boldsymbol{PA}\) получается из матрицы \(\boldsymbol A\) перестановкой строк, а матрица \(\boldsymbol{AP}\) — перестановкой столбцов;
произведение двух верхних (нижних) треугольных матриц тоже верхняя (нижняя) треугольная матрица;
если матрица \(\boldsymbol Q\) ортогональна, то \(\boldsymbol{QQ}^\mathsf{T} = \boldsymbol Q^\mathsf{T} \boldsymbol Q = \boldsymbol I\).
А вот коммутировать квадратные матрицы не обязаны: в общем случае \(\boldsymbol{AB} \ne \boldsymbol{BA}\).
A = np.random.randn(n, n)
B = np.random.randn(n, n)
np.allclose(A @ B, B @ A)
False
Powers of matrices#
Квадратные матрицы можно возводить в натуральную степень так же, как и обычные числа. Вот, например, индуктивное определение степени матрицы:
Справедливы равенства
Можно даже пойти дальше, и определить взятие экспоненты или синуса от матрицы через ряд Тейлора:
(можно доказать, что эти ряды сходятся).
Warning
В машинном (и особенно глубинном) обучении часто берут числовую функцию от матрицы, например \(\log(\boldsymbol A)\), \(\tanh(\boldsymbol A)\) или \(\sigma(\boldsymbol A)\), где \(\sigma(x) = \frac 1{1+e^{-x}}\) — сигмоида. Такая запись почти наверное подразумевает поэлементное применение функции к каждой ячейке матрицы \(\boldsymbol A\), а вовсе не матричные ряды! Тем более что матрицы в машинном обучении чаще всего прямоугольные, а для них нельзя определить ни степень, ни ряд.
Exponent in numpy
is elementwise:
np.exp(B)
array([[1.00000000e+00, 2.71828183e+00, 7.38905610e+00, 2.00855369e+01],
[5.45981500e+01, 1.48413159e+02, 4.03428793e+02, 1.09663316e+03],
[2.98095799e+03, 8.10308393e+03, 2.20264658e+04, 5.98741417e+04]])
One-rank matrix#
Одноранговая матрица задаётся произведением столбца на строку: \(\boldsymbol A = \boldsymbol {uv}^\mathsf{T}\), \(\boldsymbol u \in \mathbb R^m\), \(\boldsymbol v \in \mathbb R^n\). Элементы этой прямоугольной матрицы равны \(A_{ij} = u_i v_j\), а все её строки/столбцы пропорциональны. Отметим также, что если \(m=n\), то
Произведение матриц \(\boldsymbol A \in\mathbb R^{m\times n}\) и \(\boldsymbol B \in\mathbb R^{n\times p}\) можно записать в виде суммы одноранговых матриц:
В англоязычной литературе умножение столбца на строку \(\boldsymbol {uv}^\mathsf{T}\) называется outer product («внешнее произведение»), по аналогии со скалярным произведением \(\boldsymbol u^\mathsf{T}\boldsymbol v\), называемым inner product («внутреннее произведение»).
Block matrices#
Блочная матрица имеет вид
Блочные матрицы могут иметь произвольное число матричных блоков по каждому измерению. Интересно, что перемножать блочные матрицы можно по тем же правилам, что и матрицы из чисел, например
(при условии, что размеры матриц позволяют корректно произвести все матричные умножения).
Exercises#
Calculate \(\mathrm{tr}(\boldsymbol I_n)\).
Prove that the main diagonal of a skew-symmetric matrix contains only zeros.
Prove that the rotation matrix (45) is orthogonal.
Prove that any reflection matrix \(\boldsymbol I - 2\boldsymbol u \boldsymbol u^\top\), \(\boldsymbol u \in\mathbb R^n\), \(\Vert \boldsymbol u \Vert_2 = 1\), is orthogonal.
How many arithmetic operations is required for calculation \(\boldsymbol{Ax}\) if \(\boldsymbol A \in \mathbb R^{m\times n}\), \(\boldsymbol x \in \mathbb R^n\)? What if \(m=n\)?
Give an example of two square matrices \(\boldsymbol A\) and \(\boldsymbol B\) such that \(\boldsymbol{AB} \ne \boldsymbol{BA}\).
Prove that \(\mathrm{tr}(\boldsymbol{A} + \boldsymbol{B}) = \mathrm{tr}(\boldsymbol{A}) + \mathrm{tr}(\boldsymbol{B})\), \(\mathrm{tr}(\lambda\boldsymbol{A}) = \lambda\mathrm{tr}(\boldsymbol{A})\), \(\mathrm{tr}(\boldsymbol{AB}) = \mathrm{tr}(\boldsymbol{BA})\).
Let \(\boldsymbol A\) and \(\boldsymbol B\) be two symmetric matrices of the same shape. Is their product is necessarily symmetric?
Find the transpose to the block matrix
where \(\boldsymbol A \in \mathbb R^{k\times m}\), \(\boldsymbol B \in \mathbb R^{k\times n}\), \(\boldsymbol C \in \mathbb R^{\ell\times m}\), \(\boldsymbol D \in \mathbb R^{\ell\times n}\).