[Перевод] Семь грехов численной линейной алгебры22.10.2022 16:02

В численной линейной алгебре нас интересуют точное и эффективное решение задач и понимание чувствительности задач к возмущениям. К старту флагманского курса по Data Science делимся материалом от профессора Ника Хигэма о семи грехах линейной алгебры, из-за которых теряется точность/эффективность или информация о чувствительности [к возмущениям] оказывается недостоверной.

1. Обращение матрицы

На курсе линейной алгебры мы узнаём, что решение линейной системы $Ax = b$ , состоящей из $n$ уравнений, содержащих $n$ неизвестных, можно записать в виде $x = A^{-1}b$ , где $A^{-1}$ — обратная матрица. При этом не всегда подчёркивается, что случаев, когда следует вычислять $A^{-1}$ , очень мало. На самом деле, ( $n=1$ ) скалярную систему $7x = 21$ нельзя решить, вычислив $x = 7^{-1} \times 21$ . Вместо этого лучше выполнить деление $x = 21/7$ [прим. ред. — мы помним, что 7 в минус первой степени — это одна седьмая единицы, и, конечно, помним о свойствах умножения и деления дробей. Возможно, здесь имеется в виду разница в подходе к вычислениям в контексте вычислительной машины]. В случае $n\times n$ быстрее и точнее будет решение линейной системы (методом Гаусса) путём выбора ведущего элемента столбца, а не обращения $A$ , которое в любом случае потребует LU-разложения матрицы).

Случаи, когда требуется использовать $A^{-1}$ , редки, но возможны, когда диагональные элементы обратной ковариационной матрицы являются значимыми величинами, а также в некоторых алгоритмах вычисления матричных функций.

2. Формирование перекрёстного произведения матриц

Решение линейной задачи наименьших квадратов $\min_x\| b - Ax \|_2$ , где $A$ — матрица полного ранга $m\times n$ , где $m\ge n$ удовлетворяет нормальным уравнениям $A^T\!A x = A^Tb$ . Поэтому естественно сформировать симметричную положительно определённую матрицу $A^T\!A$ и решить нормальные уравнения с помощью разложения Холецкого. Хотя этот метод быстрый, он численно неустойчив при слабообусловленной $A$ . В противоположность этому решение задачи наименьших квадратов с разложением QR матрицы всегда численно устойчиво.

Что не так с перекрёстным произведением матриц $A^T\!A$ (известным как матрица Грама)? Она возводит данные в квадрат, что может привести к потере информации в арифметике с плавающей точкой. Например, при

$A = \begin{bmatrix} 1 & 1 \\ \epsilon & 0 \end{bmatrix}, \quad 0 < \epsilon < \sqrt{u},$

где $u$ — единичная ошибка округления для плавающей запятой, то

$A^T\!A = \begin{bmatrix} 1 + \epsilon^2 & 1 \\ 1 & 1 \end{bmatrix}$

является положительно определённой, но, так как $\epsilon^2<u$ в арифметике с плавающей точкой $1+\epsilon^2$ округляется до $1$ , получается вырожденная

$\mathrm{f\kern.2ptl}( A^T\!A) = \begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix}$ ,

при этом информация в $\epsilon$ потеряна.

Другая проблема с перекрёстным произведением матриц заключается в том, что $2$ — нормальное число обусловленности $A^T\!A$ — является квадратом такового для $A$ . Это приводит к численной нестабильности алгоритмов, которые работают с $A^T\!A$ , если число обусловленности велико.

3. Неэффективный порядок действий при определении произведения матриц

Число действий при определении произведения матриц зависит от порядка его определения (если принять, что не все матрицы $n\times n$ ). Иными словами, умножение матриц ассоциативно, поэтому $A(BC) = (AB)C$ . В общем случае число действий при определении произведения матриц зависит от того, где поставить скобки. Один порядок может быть намного лучше других, поэтому не следует просто оценивать произведение в фиксированном порядке слева направо или справа налево. Например, если $x$ , $y$ и $z$ — $n$ — векторы, то $xy^Tz$ можно определить так:

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

Словом, определение места скобок в матричном произведении $A_1A_2\dots A_k$ для минимизации числа действий — это сложная задача, но для многих случаев на практике хороший порядок определить легко.

4. Предположение, что матрица является положительно определённой

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

Определённость подразумевает, что:

диагональные элементы матрицы положительны;
$$\det(A) > 0$» />; </li> <li><img src=$ для всех $i \ne j$ .

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

Лучший способ проверить определённость, который зачастую может потребоваться сам по себе — рассчёт разложения Холецкого. функция MATLAB chol возвращает сообщение об ошибке, если разложение не удалось. При этом может быть запрошен второй выходной аргумент, в этом случае ему задаётся номер этапа, на котором случилась неудача факторизации, или 0, если факторизация прошла успешно. В случае неудачи частично вычисленный фактор $R$ возвращается в первом аргументе, и его можно использовать, например, для вычисления направления отрицательной кривизны, как это необходимо для оптимизации.

Этот грех занимает первое место в «Семи грехах оптимизации портфеля» (Seven Sins in Portfolio Optimization) Шмельцера и Хаузера, поскольку в этом случае отрицательное собственное значение в ковариационной матрице может определять портфель с отрицательной дисперсией, обещая при этом условно высокие инвестиции без риска!

5. Неиспользование структуры матрицы

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

Для задач нахождения седловой точки матрицы являются симметрично неопределёнными и имеют вид:

$\notag C = \begin{bmatrix} A & B^T \\ B & 0 \end{bmatrix},$

где $A$ является симметричной положительно определённой. Разработка числовых методов для решения $Cx = b$ с использовании блочной структуры и возможной разреженности $A$ и $B$ оказалась весьма трудоёмкой задачей. Ещё один пример — циркулянтная матрица:

$\notag C = \begin{bmatrix} c_1 & c_2 & \dots & c_n \\ c_n & c_1 & \dots & \vdots \\ \vdots & \ddots & \ddots & c_2 \\ c_2 & \dots & c_n & c_1 \\ \end{bmatrix}.$

Важное свойство циркулянтных матриц — диагонализируемость унитарной матрицей дискретного преобразования Фурье. При помощи этого свойства $Cx = v$ можно решить за $O(n \log_2n)$ , а не $O(n^3)$ операций, что потребовалось бы при игнорировании циркулянтной структуры.

В идеале программное обеспечение линейной алгебры должно обнаруживать структуру в матрице и вызывать алгоритм, использующий эту структуру. Ярким примером такого мета-алгоритма в MATLAB является функция x = A\b с обратным слешем для решения $Ax = b$ . Функция с обратным слешем проверяет, является ли данная матрица треугольной (или преобразованной треугольной) матрицей, верхней матрицей Хессенберга, симметричной или симметричной положительно определённой. После такой проверки применяется подходящий метод. При этом $A$ также может быть прямоугольной. При этом задача наименьших квадратов решается, если строк больше, чем столбцов, или недоопределённая система, если столбцов больше, чем строк.

6. Определение приближения вырожденности с помощью определителя

$n\times n$ матрица $A$ является невырожденной только при ненулевом детерминанте. Поэтому можно ожидать, что малое значение $\det(A)$ указывает на почти вырожденную матрицу. Однако размер $\det(A)$ ничего не говорит о вырожденности. Действительно, поскольку $\det(\alpha A) = \alpha^n \det(A)$ , мы можем получить любое значение определителя, умножив его на скаляр $\alpha$ , однако $\alpha A$ имеет приближение вырожденности не больше и не меньше, чем $A$ для $\alpha \ne 0$ .

Ещё одно ограничение для определителя показано на примере двух матриц:

$\notag T = \begin{bmatrix} 1 & -1 & -1 & \dots & -1\\ & 1 & -1 & \dots & -1\\ & & 1 & \dots & \vdots\\ & & & \ddots & -1 \\ & & & & 1 \end{bmatrix}, \quad U = \begin{bmatrix} 1 & 1 & 1 & \dots & 1\\ & 1 & 1 & \dots & 1\\ & & 1 & \dots & \vdots\\ & & & \ddots & 1 \\ & & & & 1 \end{bmatrix} \qquad (1)$

В обеих матрицах есть единично диагональные и внедиагональные элементы, ограниченные по модулю $1$ . Поэтому $\det(T) = \det(U) = 1$ , но

$\notag T^{-1} = \begin{bmatrix} 1 & 1 & 2 & \dots & 2^{n-2}\\ & 1 & 1 & \dots & \vdots\\ & & 1 & \ddots & 2\\ & & & \ddots & 1 \\ & & & & 1 \end{bmatrix}, \quad U^{-1} = \begin{bmatrix} 1 & -1 & & & \\ & 1 & -1 & & \\ & & 1 & \ddots & \\ & & & \ddots & -1 \\ & & & & 1 \end{bmatrix}.$
Таким образом, $T$ является слабообусловленной для большого множества $n$ . Фактически при изменении элемента $(n,1)$ матрицы $T$ до $-2^{n-2}$ матрица становится вырожденной! В противоположность этому $U$ всегда хорошо обусловлена. Определитель не позволяет отличить слабообусловленную $T$ от хорошо обусловленной $U$ .

7. Определение обусловленности при помощи собственных значений

Для любой $n\times n$ матрицы $A$ и любой нормы согласованной матрицы справедливо, что $\|A\| \ge |\lambda_i|$ для всех $i$ , где $\lambda_i$ являются собственными значениями $A$ . Поскольку собственные значения $A^{-1}$ — это $\lambda^{-1}$ , число обусловленности матрицы $\kappa(A) = \|A\| \, \|A^{-1}\|$ ограничено снизу отношением наибольшего собственному значению к наименьшему по абсолютной величине, т. е.:

$\notag \kappa(A) \ge \displaystyle\frac{ \max_i |\lambda_i| } { \min_i |\lambda_i| }.$

Но, как показала матрица $T$ в (1), эта связь может быть очень слабой.

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

$\notag \kappa_2(A) = \displaystyle\frac{\sigma_1}{\sigma_n},$

где $A = U\Sigma V^T$ — сингулярное разложение матрицы с $U$ и $V$ ортогональностью и $\Sigma = \mathrm{diag}(\sigma_i)$ , $\sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_n \ge 0$ . Если $A$ симметрична, то, например, множества $\{ |\lambda_i| \}$ и $\{\sigma_i \}$ идентичны, но в целом собственные $\lambda_i$ и вырожденные $\sigma_i$ значения могут заметно различаться.