CC BY
34
Для этого нужно показать, что матрица R невырожденная, H G Cn, E = EH, R — IEh € Fi,ie > RA = H. H
БЛАГОДАРНОСТИ: Работа выполнена при поддержке РФФИ, проект 04-07-90268.
Литература
1. Strassen V. Gaussian Elimination is not optimal. Numerische Mathematik. 1969.13, 354-356.
2. Малашонок Г. 11. О решении систем линейных уравнений р-адическим методом. Программирование, 2003, 2, с.8-22.
3. Малашонок Г. 11. Матричные методы вычислений в коммутативных кольцах. Тамбов: Изд-во ТГУ 2002.
4. Gustavson F.G. High-performance linear algebra algorithms using new generalized data structures for matrices. IBM Journal for Research and development. Vol. 41, Issue 1, January, 2003, P. 31-55.
5. Gustavson F.G. Recursion leads to automatic variable blocking for dense linear-algebra algorithms. IBM Journal for Research and development. Vol. 41, Issue 6, November, 1997, P. 737-756.
6. Axo А., Хопкрофт Дж., Ульман Дж. Построение и анализ вычислительных алгоритмов. М.: Мир, 1979.
Поступила в редакцию 15 декабря 2007г.
Вычисление матричной степени и матричных функций
© Г.И. Малашонок , © М.В. Старов
1 Введение
Обычно для вычисления матричных степеней и матричных функций используют корни характеристического полинома матрицы. Однако корни полинома вычислить точно обычно не удается. И даже задача приближенного вычисления корней полинома является трудоемкой задачей вычислительной математики, трудно решаемой для матриц большого размера.
Недавно появилась работа Виктора Робука, в которой предлагается способ вычисления матричной степени [1], в котором используются только коэффициенты характеристического полинома и не требуется вычисления его корней. Этот алгоритм позволил перевести задачу вычисления матричной степени из области вычислительной математики в область компьютерной алгебры.
Алгоритм вычисления матричной степени в работе [1] сводится к следующему. Пусть имеется матрица М порядка N и т > N, тогда вычисление матричной степени Мт можно выполнить по следующему алгоритму:
N-1
Mm = ^ Ф;,тМ* , (1)
где
N-p-1
фр,т Zp,m 'У ^ akZp+k,m •> (2)
к=1
^ (ki + ... + kN)! ki kn ro\
ZP,m = b ki! ...kN! ai ..^N (3)
где сумма пробегает все к* которые удовлетворяют уравнению:
к\ + 2к2 + ... + Мкы = т — р (4)
при условии р < т, и ад = ( — 1)9+1стд, где ад коэффициенты характеристического полинома матрицы М.
2 Оценка сложности алгоритма
Чтобы получить представление о сложности вычислений, будем искать оценку снизу числа мультипликативных операций, которые требуются в этом алгоритме.
Сначала получим оценку снизу количества решений Г(т, N) диофантового уравнения (4) при р = 0:
т т — №кМ т—Ыкн —...—3кз
N N—1 2
ЕЕ ... Е 1 (5)
кн=0 кн —1=0 к2 = 0
В этой формуле стоит N — 1 кратное суммирование целых положительных чисел, а верхние пределы - это целые части стоящих там дробей. Можно на эту сумму смотреть как на N — 1 кратный объем,
составленный из единичных кубиков. Чтобы оценить этот объем снизу, будем вычислять объем
фигуры, которая в него целиком погружена, для этого уменьшим верхний предел на 1, а целое деление заменим на рациональное деление. Получим выражение:
т л т—Ыкы л т — №кы — ... — 3кз ..
/*т— 1 г 1 С 2 1
/ / ... с1х2 ..Лхх-^ (6)
0 0 0
Последовательное прямое вычисление интегралов позволяет получить простую нижнюю оценку, если пренебречь меньшим из двух слагаемых, которое равно
кк-2 . е-20(к) к — 2 (е2\к ..
Щрк(к — 2> = 2п (Т). (7)
Здесь к > 3 - нижний индекс переменной хк, по которой ведется интегрирование. В итоге получим значение интеграла:
тм-2 е-2в(м) т — 2N + 2 ( те2
= («п2N(т — 2N+ 2> = ^---------------т— (-ж) .
Здесь использовалась формула Стирлинга:
—! = /2пп(—') ев(п\ где |0(—)| < ----. (8)
\е/ 12—
Для вычисления всех требуемых значений Ф^т из (1) нужно знать все значения Zim, i = 0, N — 1 (3). Если даже предполагать, что в (3) все факториалы и необходимые степени коэффициентов характеристического полинома заранее вычислены и записаны в память, то тогда для вычисления одно такого слагаемого Zim потребуется не меньше, чем 2N операций умножения. Следовательно, для нахождения всех значений Zim требуется
m
Rp = 2N ^ Fk,N
k=m-N +1
операций умножения.
Вычисление всех Ф^т требует дополнительно N(N2—~ операций умножения, вычисление всех матричных степеней Mk, k < N требует N3(N — 2) операций умножения и вычисления, в соответствии с (1) требуют N2 (N — 1) операций умножения. Всего дополнительно потребуется N4 — N3 — N2/2 — N/2 операций.
Следовательно, основной вклад в число мультипликативных операций дает сумма (6) со сла-
N
сложность превышает полиномиальную, считаются очень сложными. И их стараются не применять.
Поэтому в следующем параграфе мы предлагаем улучшенный алгоритм, который имеет полиномиальную сложность.
3 Алгоритм с полиномиальной сложностью
М
^(Л) = А" + ^Л"-1 + ... + Нм . (9)
Пусть т > N рассмотрим деление с остатком полинома Ат на Л.(Л):
Лт = / (Л)Л-(Л) + г(Л), (10)
где /(Л) = /оЛт-" + /1Лт-"-1 + ... + /т_" и г(Л) = г 1 Л"-1 + ... + Г".
По теореме Гамильтона-Кэли Л.(М) = 0. Следовательно:
Мт = г(М). (11)
Таким образом, задача нахождения матричной степени сводится к нахождению коэффициентов по-Г(Л)
коэффициенты, приравняем коэффициенты при равных степенях в левой и правой части равенства (10). Получим систему, состоящую из т + 1 уравнений и содержащую т + 1 неизвестных
/о,.., /т_м, го,.., Г"_1:
/о = 1 ,
fi = — (h1fi —1 + h2fi —2 + ... + hi/o) , i = 1 N — 1 (12)
fi = —(h1fi —1 + h2fi —2 + ... + hN fi — N) , i = N m — N +1
ri = (hi fm—N + ... + hNfm —2N+i) , i = 1, N
Матрица коэффициентов системы имеет треугольный вид, у которой все элементы, за исключением N
этой системы потребуется всего N(m — N) операций умножения. Учитывая то, что для вычисления всех степеней матрицы Mk, k < N, требуется N3(N — 2) операций умножения, а для вычнсле-
N2 ( N — 1)
умножения, всего получим C = N(N3 + m) — N2(N + 2) операций умножения во всем алгоритме. Приведем соответствующий алгоритм.
Algorithm MatrixPower
Mm
Mm
N = length(M);
C oef C harPol [] = CoeffCharPol(M); n = m — N + 1; f [0] = 1;
i = 1 n — 1 k = N — 1;
for j = i — N to i — 1 do if j > 0 then
f[i] = f [i] + CoefCharPol[k] x (f [j]); k = k — 1; f [i] = —f [i];
i = 0 N — 1 k = m — N +1; j = i N— 1 k = k — 1; if k > 0 then
r[i] = r[i] + CoefCharPol[j] x f [k]; r[i] = —r[i];
MatrP [0] = I;
MatrP [1] = M;
i = 2 N — 1 MatrP [i] = M x MatrP [i — 1];
i = 0 N — 1 MatrP[i] = MatrP[i] x R[N — i — 1]; MatrPower = MatrPower + MatrP [i]; return MatrPower;
4 Улучшенный алгоритм с экспоненциальной сложностью
Исследование решения системы (12) позволяет по-новому взглянуть на алгоритм с экспоненциальной сложностью.
Рассмотрим решение системы (12) методом Крамера. Определитель матрицы коэффициентов системы будет иметь вид:
д
і G. .G G. .. G G G. .G G
hl і. .G G. .. G G G. .G G
hN hN — l . .і G. .. G G G. .G G
G hN . . hl і. .. G G G. .G G
G G. . hN hN — l . .. і G G. .G G
G G. .G hN . . hl і G. .G G
G G. .G G. . . h2 G і. .G G
G G. .G G. . . hN—l G G. .і G
G G. .G G. . . hN G G. .G 1
Так как этот определитель равен 1, то последние N неизвестных Г1, ..,Г" равны соответствующим определителям, полученным заменой одного из последних столбцов на единичный столбец. Раскладывая их по последним N столбцам, найдем г* (1=1,..,1М):
hl і. .. G G. .. G . . . G G
hN hN — l . .. і G. .. G . . . G G
G hN . . . hl і. .. G . . . G G
G G. . . hN hN — l . .. і . . . G G
G G. .. G hN . . . hl . . . G G
G G. .. G G. . hN—i . . . hl 1
G G. .. G G. . . hN . . hi+1 hi
Эти определители отличаются только последней строкой: определитель г* в последней строке имеет N — I + 1 коэффициент характеристического полинома ^,.., а их первые т — N строки берутся
из определителя Д. Разложим г* по последней строке. Получим:
N
ri ^ ^ ak Zm—N+i — ^ k=i
где а4 = ( —1)4+1 ^ и ^ обозначает минор, образованный строками 1, 2,.., т — N и столбцами
1, 2,.., t — 1,4 + 1,.., т — N + 1 определителя Д.
r
Отметим, что = 1 и Zk, к > 1, равен значению верхнего детого углового минора порядка к — 1 определителя Д. Докажем по индукции, что:
(к1 + ... + )! к1 к
— ' ------------------------«1 ...а
к1! • • • !
N
N
для к = 2,.., т — N Здесь сумма пробегает все к*, которые удовлетворяют уравнению
к1 + 2к2 + ... + Жк^ = к — 1 .
Предположим, что это утверждение верно для всех Zj, г < к. Рассматривая Zk как угловой минор к—1
N
Zfc ^ ^ aiZfc—г
і=1
Здесь предполагается, что ^ = ^и і < 1.
N N
^ ^ «і^к-і ^ ^ «і ^ ^
і=1 і=1 к—і
(к1 + ... + ^ )! кі
к^ ••• kN!
а ... а
N
ЕЕ
і=1 к
(к1 + ... + кі — 1 + ... + kN)! к1
----------------------------------а 1
к1! ••• (кі - 1)! ••• kN! 1
N
к N ■ «/ =
^ (к1 + ... + ^)! к1 кн ^-кї-ттгк^-^ ...«
N
Здесь мы воспользовались тождеством
N
Е
і=1
(к1 + ... + кі — 1 + ... + kN)! (к1 + ... + kN)!
к1! ••• (кі — 1)! ••• kN!
к1! • • • kN!
Таким образом мы получили еще один алгоритм вычисления матричной степени
N N
Мт = +і-кмN-і .
і=1 к=і
Этот алгоритм напоминает алгоритм (1), однако тут нет лишних слагаемых, которые имеются в сумме (1). Соответствие этих алгоритмов устанавливается равенством = Ят-р. Этот алгоритм также имеет экспоненциальную сложность.
к
N
5 Численные эксперименты
Были реализованы программы для приведенных выше алгоритмов с полиномиальной и экспоненциальной сложностью вычисляющих матричную степень. Приведем результаты экспериментов по вычислению матричных степеней для матриц пятого порядка.
Степень Ді(мс) Жі(мс) і Л£ 1 ~ /V* Др Жр р N *
25 16.7 0.6 27.8 3745 575 6
50 388 0.7 554 7.8 х 104 700 111
75 2847 0.9 3163 4.4 х 105 825 533
100 11536 1.2 9613 1.5 х 106 950 1550
125 34599 1.5 23066 3.7 х 106 1075 3459
150 88902 2.1 42334 7.8 х 106 1200 6570
Np — число операций умножения в алгоритме с полиномиальной сложностью,
Rp — нижняя оценка числа операций умножения в алгоритме с экспоненциальной сложностью,
Nt — время вычисления матричной степени в алгоритме с полиномиальной сложностью,
Rt — время вычисления матричной степени при использовании алгоритма с экспоненциальной сложностью,
NRt — время вычисления матричной степени при использовании алгоритма с полиномиальной сложностью,
NRp
ностью.
Литература
1. Robuk V.N. A constructive formula for function of matrix. Alternative to the Lagrange-Silvester formula. Nuclear Instruments and Methods in Physics Research A 534 (2004) 319-323.
2. Ланкастер П. Теория матриц. М.: Наука, 1973, с. 280.
3. Математическая энциклопедия/ Гл. ред. И.М. Виноградов. — М.: Сов. Энциклопедия, 1979. Т.5.
Поступила в редакцию 15 декабря 2007г.
Использование жесткого диска в матричных вычислениях
© М.С.Зуев
Вычисления с использованием жесткого диска реализованы, например, в известном пакете численных методов POOCLAPACK (Parallel out-of-core Linear algebra package, см. [5]).
В данном пакете реализованы два способа хранения матриц на диске, один из которых повторяет способ хранения матриц в ОЗУ по строкам или столбцам, а второй предполагает хранение матрицы по блокам. Во втором случае требуется использование блочных матричных алгоритмов.
В современной литературе блочные и блочно-рекурсивные матричные методы рассматриваются как одни из самых эффективных. В отличие от стандартных строчно-ориентированных методов, такие методы более эффективно используют процессорный кэш, и поэтому позволяют достигать производительности процессора, близкой к пиковой (см. [2-4]). Например, в работе [2] указывается достижение 92 % производительности процессора IBM Power 3 с блочной процедурой разложения Холецкого и 76 % с соответствующей поэлементной процедурой.
Если процедура использует жесткий диск, то она может эффективно использовать оперативную память. Считывание блоков из файла целиком и возможность варьировать их размер позволяет определить баланс между экономией оперативной памяти и минимизацией количества обращений к диску. При считывании блока из файла требуется одно обращение к таблице размещения файлов (суперблоку) и одно обращение к файлу. Если бы блоки хранились в виртуальной памяти, то для их считывания требовалось бы большее количество обращений к диску.
Суммируя все сказанное выше, можно предположить, что блочные матричные алгоритмы с использованием жесткого диска могут быть эффективными.
Для исследования эффективности таких алгоритмов были проведены вычислительные эксперименты на однопроцессорной машине. Сравнивались алгоритмы умножения и обращения матриц (одностороннее обращение, см. [1]) размера 4096 х 4096 в конечном поле по 28-битному простому модулю. Как и в POOCLAPACK, матрица разбивается на блоки фиксированного размера, которые записываются в отдельные файлы. Матрице сопоставляется папка на диске. Если разбить матрицу на четыре равных блока, то возможно их хранение в отдельных папках.
Для определения лучшего алгоритма умножения были рассмотрены четыре алгоритма: стандартный блочный алгоритм (алгоритм 0), стандартный блочный алгоритм с записью блоков на
