Первая страница / Теория / Компьютерная математика /

Матрицы и действия с ними

Голосование: 35, 9

Введение

Матрицы часто встречаются в научных расчетах, поэтому важно уметь эффективно с ними работать. Данная статья содержит краткое введение в теорию матриц и операций над ними. Особое внимание уделяется умножению матриц и решению систем линейных уравнений. Мы введем основные понятия и обозначения, затем изложим алгоритм Штрассена, позволяющий умножить две матрицы размера n×n за время Θ(n2.81), появление которого стало в свое время неожиданностью. Далее будет рассказано, как решать системы линейных уравнений с помощью так называемого LUP-разложения, и обсуждена задача обращения матриц.

Определения

Матрица элементов из некоторого пространства S размера m×n — это объект из пространства Sm×n, координаты которого упорядочены по строкам и столбцам. Далее мы будем рассматривать вещественные матрицы, которые можно считать прямоугольными таблицами чисел.

При транспонировании матрицы A ее строки становятся столбцами и наоборот. Матрица, получаемая из A транспонированием, обозначается AT, например для матрицы

Formula

получаем

Formula

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

Formula

является вектором из трех элементов. Стандартной формой вектора мы будем считать вектор-столбец, т. е. n×1 матрицу; при его транспонировании получается вектор-строка

xT = (2 3 5).

Вектор, i-й элемент которого равен 1, а все остальные элементы равны 0, называют единичным вектором и обозначают ei. Количество элементов единичного вектора обычно определяется из контекста.

Нулевой матрицей называется матрица, все элементы которой равны 0. Такая матрица обычно обозначается O. Ее размер также обычно определяется из контекста.

Часто встречаются квадратные матрицы — матрицы размера n×n. Рассмотрим некоторые виды таких матриц:

  1. У диагональной матрицы все внедиагональные элементы равны нулю (aij = 0 при i <> j), поэтому она может быть задана перечислением элементов, стоящих на диагонали.
    Formula
  2. Единичной матрицей называется диагональная матрица, диагональ которой заполнена единицами. При этом столбцами такой матрицы служат векторы e1e2, …, en.
    Formula
  3. У верхне-треугольной матрицы все элементы под главной диагональю равны нулю (uij = 0 при i > j):
    Formula
  4. А у нижне-треугольной матрицы все элементы над главной диагональю равны нулю (lij = 0 при i < j):
    Formula
  5. Матрица перестановки имеет в точности одну единицу в каждой строке и в каждом столбце; на всех прочих местах у нее стоят нули. Пример матрицы перестановки:
    Formula
    Умножение вектора x на матрицу перестановки приводит к перестановке его элементов.
  6. Симметрическая матрица удовлетворяет условию A = AT. Такая матрица, например:
    Formula

Действия с матрицами

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

Пусть даны (m×n)-матрицы A = (aij) и B = (bij). Назовем их суммой (m×n)-матрицу C = (cij) = A + B с элементами

cij = aij + bij,

где i = 1, 2, …, m и j = 1, 2, …, n. Другими словами, сложение матриц осуществляется покомпонентно. Нулевая матрица является нейтральной матрицей для операции сложения матриц:

A + O = A = O + A.

Пусть λ — число, а A = (aij) — матрица. Можно умножить матрицу A на число λ, умножив каждый элемент матрицы A на число λ, в результате получится матрица λA = (λaij). Особо отметим матрицу −A = (−1)A, называемую противоположной к A матрицей. Элемент с индексами ij в матрице −A равен −aij, поэтому

A + (−A) = O = (−A) + A.

Вычитание матрицы мы теперь можем определить как прибавление противоположной матрицы: A − B = A + (−B).

Умножение матрицы A на матрицу B осуществимо, лишь если они имеют согласованные размеры, то есть число столбцов A совпадет с числом строк B. Если A = (aij) является (m×n)-матрицей, а B = (bij) является (n×p)-матрицей, то их произведением C = AB называется (m×p)-матрица C = (cij), в которой

Formula

для i = 1, 2, …, m и k = 1, 2, …, p. Матрицы обладают многими (хотя и не всеми) свойствами чисел. Единичная матрица является нейтральным элементом для умножения: для любой (m×n)-матрицы A:

ImA = AIn = A.

Умножение на нулевую матрицу дает нулевую матрицу:

AO = O.

Умножение матриц ассоциативно:

A(BC) = (AB)C

для любых матриц AB и C согласованных размеров. Умножение матриц дистрибутивно относительно сложения:

A(B + C) = AB + AC,
(B + C)D = BD + CD.

При n <> 1 умножение (n×n)-матриц, вообще говоря, не коммутативно, что легко проверяется непосредственно. Например, для

Formula

имеем

Formula

Рассматривая вектор-столбец как (n×1)-матрицу, а вектор-строку как (1×n)-матрицу, мы можем перемножать векторы и матрицы. Скалярным произведением векторов x и y называют число, представляющее собой

Formula

Тензорным произведением этих же векторов называется матрица Z = xyT размера (n×n) с элементами zij = xiyj. Евклидова норма ||x|| вектора x определяется равенством

Formula

в евклидовом пространстве.

Способы умножения матриц. Метод Штрассена

Если перед нами стоит задача получить произведение C двух матриц A и B размера n×n, то это можно сделать несколькими способами. Лобовой алгоритм, умножающий по формуле cik = Σaijbjk, работает за время Θ(n3) = Θ(nlog 8). Другой простой алгоритм работает рекурсивно и разбивает каждую из умножаемых матриц на 4 части так, что произведение C = AB представимо следующим образом:

Формула простого рекурсивного умножения матриц

Он производит 8 умножений матриц размером n⁄2, поэтому также работает за время Θ(nlog 8). Алгоритм Штрассена, реализованный в визуализаторе, обозначает схему, при которой в рекурсивном алгоритме можно обойтись семью умножениями матриц размера n⁄2, тем самым сократив трудоемкость до Θ(nlog 7) = O(n2.81), что дает заметный выигрыш на больших плотных матрицах (начиная, примерно, от 64×64).

Описание алгоритма

Вначале мы проверяем, является ли размер умножаемых матриц n натуральной степенью двойки. Если нет, мы дополняем исходные матрицы дополнительными нулевыми строками и столбцами, оставляя на главной диагонали единицы. При этом мы получаем удобные для рекурсивного умножения размеры, но теряем в эффективности за счет дополнительных ненужных умножений. Алгоритм Штрассена умножает две (n×n)-матрицы A и B следующим образом:

  1. Каждая из матриц A и B разбивается на 4 блока по схеме, приведенной выше.
  2. Строятся 14 матриц A1, B1, A2, B2, …, A7, B7 размера (n/2×n⁄2) (для чего нужно Θ(n2) операций сложения/вычитания чисел).
  3. Рекурсивно вычисляются 7 произведений матриц меньшего размера Pi = AiBi (i = 1, …, 7).
  4. Вычисляются части r, s, t, u искомой матрицы C. Они являются линейными комбинациями матриц Pi с коэффициентами из множества {−1, 0, 1}, и вычисление их требует Θ(n2) операций сложения/вычитания чисел.

Формулы, по которым происходит умножение

iAiBiPi = AiBi
1ag − hag − ah
2a + bhah + bh
3c + dece + de
4df − edf − de
5a + de + hae + ah + de + dh
6b − df + hbf + bh − df − dh
7a − ce + gae + ag − ce − cg
r = P5 + P4 − P2 + P6 = ae + bf
s = P1 + P2 = ag + bh
t = P3 + P4 = ce + df
u = P5 + P1 − P3 − P7 = cg + dh

Видно, что алгоритм Штрассена действительно работает! Также доступен визуализатор [2], с помощью которого можно посмотреть работу алгоритма пошагово. Тем не менее, самый асимптотически быстрый из известных на сегодняшний день способ умножения матриц был опубликован Копперсмитом и Виноградом в 1987-м году и позволяет перемножать квадратные матрицы за время Θ(n2.38). На практике же метод Штрассена гораздо проще программируется и имеет меньшую константу в оценке трудоемкости.

Обратная матрица, ранг и детерминант

Матрицей, обратной к (n×n)-матрице A, называется матрица A−1, для которой AA−1 = In = A−1A (если таковая существует). Например,

Formula

Многие ненулевые (n×n)-матрицы не имеют обратных — такие матрицы называются необратимыми или вырожденными. Пример ненулевой вырожденной матрицы:

Formula

Матрица, имеющая обратную, называется обратимой или невырожденной. Если обратная матрица существует, то она единственна. Отметим еще два интересных свойства обратных матриц: если (n×n)-матрицы A и B обратимы, то

(BA)−1 = A−1B−1;

операция обращения матрицы перестановочна с операцией транспонирования —

(A−1)T = (AT)−1.

Говорят, что векторы x1x2, …, xn линейно зависимы, если найдется набор коэффициентов c1c2, …, cn, не все из которых равны нулю, для которого c1x1 + c2x2 + … + cnxn = 0. Например, векторы (1, 2, 3)T, (2, 6, 4)T и (4, 11, 9)T линейно зависимы, поскольку 2x1 + 3x2 − 2x3 = 0. Векторы, не являющиеся линейно зависимыми, называются линейно независимыми. Таковы, например, столбцы единичной матрицы.

Столбцовым рангом ненулевой (m×n)-матрицы A называется наибольшее число линейно независимых столбцов A. Строчным рангом называется наибольшее число линейно независимых строк. Для любой матрицы A эти два числа совпадают, так что их общее значение называется просто рангом. Ранг (m×n)-матрицы представляет собой целое число в пределах от 0 до min(nm). Ранг нулевой матрицы равен нулю, а ранг единичной матрицы In равен n. Квадратная (n×n)-матрица ранга n называется матрицей полного ранга. Основное свойство рангов, доказываемое в любом курсе линейной алгебры, таково: квадратная матрица имеет полный ранг тогда и только тогда, когда она невырождена.

Минором элемента aij матрицы A размера (n×n) называется ((n−1)×(n−1))-матрица A[ij], получаемая вычеркиванием i-й строки и j-го столбца в матрице A. Один из способов ввести понятие определителя — это рекурсивная формула

Formula

Множитель (−1)i+jdet(A[ij]) называется алгебраическим дополнением элемента aij. Отметим следующие важные свойства определителя квадратной матрицы A:

  • Если в какой-либо строке или каком-либо столбце матрицы стоят одни нули, то ее определитель равен 0;
  • Если умножить все элементы какой-либо строки матрицы на некоторое число λ, то ее определитель также умножится на λ (аналогично для столбцов);
  • Если прибавить к элементам одной строки соответствующие элементы другой, то определитель не изменится (аналогично для столбцов);
  • Определители матриц A и AT равны;
  • При перестановке двух строк или столбцов матрицы ее определитель меняет знак.

При этом квадратная матрица A вырождена тогда и только тогда, когда det(A) = 0. И если B — квадратная матрица совпадающего с A размера, то det(AB) = det(A)det(B).

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

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

Рассмотрим систему n линейных уравнений с n неизвестными x1x2, …, xn (число уравнений равно числу неизвестных):

Formula

Набор значений x1x2, …, xn, обращающий в равенство каждое из уравнений системы (*), называется решением системы (*). Перепишем систему (*) в виде матричного уравнения

Formula

или, введя обозначения A = (aij), x = (xj), b = (bi), в виде

Ax = b  (**)

Если матрица A невырождена, то для нее найдется обратная матрица A−1 и вектор

x = A−1b

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

Итак, у нас есть система Ax = b из n уравнений с n неизвестными. Можно было бы попытаться вычислять обратную матрицу A−1 и умножать ее на вектор b. Однако на практике такой подход часто сталкивается с вычислительной неустойчивостью: вещественные числа хранятся в памяти лишь приближенно и ошибки округления могут накапливаться. Другой подход — так называемое LUP-разложение — в меньшей степени подвержен этой опасности, к тому же, работает он раза в три быстрее. По сути метод LUP-разложения является более формализованным широко известным методом Гаусса, решающим системы линейных уравнений за Θ(n3), но более удобен, так как для систем с одинаковыми матрицами, но разными векторами b, достаточно один раз выполнить LUP-разложение матрицы A с кубической трудоемкостью, и впоследствии решать все системы уже с квадратичной трудоемкостью.

Общая схема LUP-разложения

Мы строим три (n×n)-матрицы LU и P, образующие LUP-разложение матрицы A в следующем смысле:

PA = LU,

причем

  • L является нижне-треугольной матрицей, и ее главная диагональ заполнена единицами,
  • U является верхне-треугольной матрицей, и
  • P является матрицей перестановки.

Всякая невырожденная матрица A обладает LUP-разложением.

Зачем нужно LUP-разложение? Зная его, мы сводим систему (**) к двум системам с треугольными матрицами, а такие системы решаются просто. Вот как происходит сведение. Умножим обе части уравнения Ax = b слева на матрицу P, получаем уравнение PAx = Pb (умножение слева на P соответствует перестановке уравнений). Имеем:

LUx = Pb.

Остается решить две системы с треугольными матрицами. Сначала из нижнетреугольной системы

Ly = Pb

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

Ux = y

найдем x с помощью "обратной подстановки". Вектор x будет искомым решением: матрица P обратима, поэтому

Ax = P−1LUx = P−1Ly = P−1Pb = b.

Начнем с прямой и обратной подстановок.

Прямая и обратная подстановки

Прямая подстановка позволяет решить нижнетреугольную систему Ly = Pb за время Θ(n2). Через π обозначим перестановку, соответствующую матрице P: Pi,π[i] = 1, а Pij = 0 при j <> π[i]. Тогда в матрице PA на пересечении i-й строки и j-го столбца стоит число aπ[i],j, а i-й элемент вектора Pb равен bπ[i]. Матрица L является нижне-треугольной с единицами на главной диагонали, поэтому система Ly = Pb имеет вид:

Formula

Первое уравнение дает нам значение

y1 = bπ[1],

подставив его во второе уравнение, получим

y2 = bπ[2] − l21y1.

Теперь, подставив y1 и y2 в третье уравнение, получим

y3 = bπ[3] − (l31y1 + l32y2).

И так далее: для получения yi подставим y1y2, …, yi−1 в i-е уравнение и получим

Formula

Обратная подстановка для уравнения Ux = y осуществляется аналогично, начиная с последнего уравнения, и также требует времени Θ(n2). Матрица U — верхне-треугольная, и система Ux = y имеет вид

Formula

Мы последовательно вычисляем xnxn−1, …, x1 по формулам

Formula

или, в общем виде:

Formula

Процедура LUP-Solve находит x по данным PL, U и b, делая прямую, а затем и обратную подстановку. В тексте процедуры n выражено как rows[L], а матрица перестановки задана массивом π.

LUP-Solve(L, U, n, b)
1 n ← rows[L]
2 for i ← 1 to n
3       do yi ← bπ[i] − Σij=11(lijyj)
4 for i ← n downto 1
5       do xi ← (yi − Σjn = i + 1(uijxj)) / uii
6 return x

В строках 2-3 проводится прямая подстановка для определения y, а в строках 4-5 — обратная (для x). Внутри каждого из двух циклов for неявно присутствуют циклы суммирования, поэтому время работы процедуры составляет Θ(n2).

Вычисление LUP-разложения

Остается научиться строить LUP-разложение невырожденной (n×n)-матрицы A. Наш способ построения LUP-разложения называется методом Гаусса исключения неизвестных. Сначала мы исключим из всех уравнений, кроме первого, неизвестное x1. Для этого вычтем из всех остальных уравнений первое, умноженное на подходящее число (разное для разных уравнений). Далее устраним из третьего и последующих уравнений переменную x2, вычитая второе с подходящими множителями. Продолжая эту процедуру, мы приведем систему к верхне-треугольному виду. Матрица полученной системы и будет искомой матрицей U. Матрицу L мы составим из множителей, участвовавших в процессе исключения неизвестных. Будем называть элементы, на которые мы делим в процессе разложения, ведущими элементами. Матрица P соответствует перестановке уравнений, которая предотвращает деление на 0 (или на очень маленькое число) в описанном процессе. Этот прием называется выбором ведущего элемента — мы будем выбирать наибольший по модулю элемент.

Опишем рекурсивный метод построения LUP-разложения. На первом шаге сначала переместим ненулевой элемент из первого столбца матрицы A — пусть это ak1 — в левый верхний угол матрицы. Переставим строки 1 и k в матрице A (что соответствует перестановке уравнений), умножив A слева на матрицу перестановки Q. Тогда

Formula

где v = (a21a31, …, ak−1,1a11ak+1,1, …,  an1)T (k-й элемент заменен на первый), wT = (ak2, …, akn), и матрица A′ имеет размер (n−1)×(n−1). Теперь ak1 <> 0 и мы можем представить следующим образом матрицу A, не опасаясь деления на 0:

Formula

Это равенство проверяется непосредственным умножением. Матрица A′ − vwT / ak1 называется дополнением Шура элемента a11 матрицы A и является невырожденной матрицей. По индуктивному предположению построим для ((n−1)×(n−1))-матрицы A′ − vwT / ak1 LUP-разложение, найдя нижне-треугольную матрицу L′ с единицами на главной диагонали, верхне-треугольную матрицу U′ и матрицу перестановки P′, для которых

P′(A′ − vwT / ak1) = LU′.

Матрица

Formula

будет матрицей перестановки, как произведение двух матриц перестановки. Искомым LUP-разложением будет

Formula

Матрица L′ — нижне-треугольная с единицами на главной диагонали, а U′ — верхне-треугольная, и потому матрицы L и U также обладают этим свойством.

В процедуре LUP-Decomposition, осуществляющей LUP-разложение, мы заменим рекурсию циклом. Матрицу перестановки P мы представляем массивом π, считая что в i-й строке P единица стоит на π[i]-м месте. Матрицы L и U будут вычисляться на месте матрицы A, и после работы процедуры

Formula
LUP-Decomposition(A)
1 n ← rows[A]
2 for i ← 1 to n
3       do π[i] ← i
4 for k ← 1 to n
5       do p ← 0
6             for i ← k to n
7                   do if |aik| > p
8                         then p ← |aik|
9                              k′ ← i
10            if p = 0
11                  then error "вырожденная матрица"
12            π[k] ↔ π[k′]
13            for i ← 1 to n
14                  do aki ↔ ak′i
15            for i ← k + 1 to n
16                  do aik ← aik / akk
17                        for i ← k + 1 to n
18                              do aij ← aij − aikakj

Вначале (строки 2-3) в массив π помещается тождественная перестановка. Внешний цикл, начинающийся в строке 4, заменяет рекурсию. При каждом значении k в строках 5-9 определяется наибольший по модулю элемент первого столбца текущей матрицы (размера  (n − k + 1) × (n − k + 1)), им будет akk. Если там одни нули, то матрица вырождена, и программа сообщает об ошибке (строки 10-11). В строках 12-14 мы переставляем наибольший элемент в угол (модифицируем массив π, меняем местами строки k и k′). Ведущим элементом становится akk. (Поскольку на P′ умножается не только дополнение Шура A′ − vwT / ak1, но и вектор v / ak1, строки нужно переставлять целиком.) Наконец, в строках 15-18 вычисляется дополнение Шура. Максимальное число вложенных циклов равно 3, так что время работы процедуры LUP-Decomposition есть Θ(n3). Выбор ведущего элемента увеличивает время работы лишь в ограниченное число раз.

Обращение матриц

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

Пусть нам дано LUP-разложение PA = LU матрицы A размера (n×n). Зная его, можно решить систему вида Ax = b с помощью процедуры LUP-Solve за время Θ(n2). Чтобы решить другую систему Ax = b′ (с той же матрицей, но другой правой частью), нам понадобится еще столько же времени. Таким образом, зная LUP-разложение матрицы A, можно решить k систем с матрицей A за время Θ(kn2).

Матричное уравнение AX = In для обратной матрицы X можно рассматривать как совокупность n систем вида Ax = b. Обозначим i-й столбец матрицы X через Xi; тогда AXi = eii = 1, …, n, поскольку i-м столбцом матрицы In является единичный вектор ei. Другими словами, нахождение обратной матрицы сводится к решению n уравнений с одной матрицей и разными правыми частями. После выполнения LUP-разложения (время Θ(n3)) на решение каждого из n уравнений нужно время Θ(n2), так что и эта часть работы требует времени Θ(n3).

Литература

  1. Кормен Т., Лейзерсон Ч., Ривест Р. Алгоритмы: построение и анализ. — М.: МЦНМО, 1999.

Визуализаторы

  1. Котов А. Алгоритм Штрассена (2004)

Александр Котов


Долгачев Георгий / 2008-10-07 11:02:26

Можно ли получить ответ на следующие вопросы?

Пусть дана квадратная матрица A. Из столбцов матрицы A создана матрица B, в которой столбцы располагаются в обратном порядке:

b[i,j] = a[i, N-j+1],

i = 1..N;

j = 1..N;

где: i - индекс строки;

j - индекс столбца;

N - размерность матриц.

Вопросы:

- имеет ли операция получения матрицы B общепринятое наименование, аналогичное операции "транспонирование"?

Возможно, "зеркальное отражение ..."?

- имеет ли сама матрица "B" общепринятое наименование?

Тогда, соответственно, "зеркально-отраженная ...".

Тупик Н. В. / 2009-04-23 19:00:42

А есть ли алгоритм перемножения двух матриц 2х2 с 7 умножениями и 15 сложениями/вычитаниями. Стандартный алгоритм Штрассена даёт 18 сложений/вычитаний

Евгений / 2010-02-07 15:36:43

Очень понравилось. Всё чётко и точно.

Алекс / 2010-12-17 13:31:12

Тоже очень понравилось, только поздно вас нашел..(

Что ж, есть разные алгоритмы поиска. Видимо, Ваш работал медленно. ;)

Ваше имя
Email
Текущий день недели (строчными буквами)
Комментарий