close

Вход

Забыли?

вход по аккаунту

?

Вычисление матричной степени и матричных функций.

код для вставкиСкачать
Для этого нужно показать, что матрица R невырожденная, H ? Cn , E = EH , R ? IЇEH ? FI,IEH ,
RA = H .
БЛАГОДАРНОСТИ: Работа выполнена при поддержке РФФИ, проект 04-07-90268.
Литература
1. Strassen V. Gaussian Elimination is not optimal. Numerische Mathematik. 1969. 13, 354-356.
2. Малашонок Г.И. О решении систем линейных уравнений р-адическим методом. Программирование, 2003, 2, с.8-22.
3. Малашонок Г.И. Матричные методы вычислений в коммутативных кольцах. Тамбов: Изд-во
ТГУ 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. 3155.
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. 737756.
6. Ахо А., Хопкрофт Дж., Ульман Дж. Построение и анализ вычислительных алгоритмов.
М.: Мир, 1979.
Поступила в редакцию 15 декабря 2007г.
Вычисление матричной степени
и матричных функций
c
1
c М.В. Старов
Г.И. Малашонок , Введение
Обычно для вычисления матричных степеней и матричных функций используют корни характеристического полинома матрицы. Однако корни полинома вычислить точно обычно не удается.
И даже задача приближенного вычисления корней полинома является трудоемкой задачей вычислительной математики, трудно решаемой для матриц большого размера.
Недавно появилась работа Виктора Робука, в которой предлагается способ вычисления матричной степени [1], в котором используются только коэффициенты характеристического полинома и не
требуется вычисления его корней. Этот алгоритм позволил перевести задачу вычисления матричной
степени из области вычислительной математики в область компьютерной алгебры.
Алгоритм вычисления матричной степени в работе [1] сводится к следующему. Пусть имеется
матрица M порядка N и m ? N , тогда вычисление матричной степени M m можно выполнить по
следующему алгоритму:
N
?1
X
Mm =
?i,m M i ,
(1)
i=0
где
?p,m = Zp,m ?
NX
?p?1
?k Zp+k,m ,
(2)
k=1
Zp,m =
X (k1 + . . . + kN )!
k1 ! . . . kN !
kN
?1k1 . . . ?N
(3)
где сумма пробегает все ki которые удовлетворяют уравнению:
k1 + 2k2 + . . . + N kN = m ? p
(4)
при условии p ? m, и ?q = (?1)q+1 ?q , где ?q коэффициенты характеристического полинома матрицы M .
2
Оценка сложности алгоритма
Чтобы получить представление о сложности вычислений, будем искать оценку снизу числа мультипликативных операций, которые требуются в этом алгоритме.
Сначала получим оценку снизу количества решений F (m, N ) диофантового уравнения (4) при
p = 0:
m
N
X
m?N kN
N ?1
X
m?N kN ?...?3k3
2
X
...
kN =0 kN ?1 =0
1
(5)
k2 =0
В этой формуле стоит N ?1 кратное суммирование целых положительных чисел, а верхние пределы
это целые части стоящих там дробей. Можно на эту сумму смотреть как на N ? 1 кратный объем,
составленный из единичных кубиков. Чтобы оценить этот объем снизу, будем вычислять объем
фигуры, которая в него целиком погружена, для этого уменьшим верхний предел на 1, а целое
деление заменим на рациональное деление. Получим выражение:
kN
Z m?N
Z m?N kN2?...?3k3 ?1
Z m
?1
N ?1
N ?1
...
dx2 . . . dxN ?1 dxN
(6)
0
0
0
Последовательное прямое вычисление интегралов позволяет получить простую нижнюю оценку,
если пренебречь меньшим из двух слагаемых, которое равно
k k?2
e?2?(k) k ? 2 e2 k
.
(7)
k(k
?
2)
=
(k!)2
2?
k2
k
Здесь k ? 3 нижний индекс переменной xk , по которой ведется интегрирование. В итоге получим
значение интеграла:
e?2?(N ) m ? 2N + 2 me2 N
mN ?2
.
N (m ? 2N + 2) =
Fm,N =
2
(N !)
2?
m2
N2
Здесь использовалась формула Стирлинга:
n n
?
1
e?(n) , где |?(n)| <
n! = 2?n
.
e
12n
(8)
Для вычисления всех требуемых значений ?i,m из (1) нужно знать все значения Zi,m , i = 0, N ? 1
(3). Если даже предполагать, что в (3) все факториалы и необходимые степени коэффициентов
характеристического полинома заранее вычислены и записаны в память, то тогда для вычисления
одно такого слагаемого Zi,m потребуется не меньше, чем 2N операций умножения. Следовательно,
для нахождения всех значений Zi,m требуется
Rp = 2N
m
X
Fk,N
k=m?N +1
операций умножения.
Вычисление всех ?i,m требует дополнительно N (N2?1) операций умножения, вычисление всех
матричных степеней M k , k < N требует N 3 (N ? 2) операций умножения и вычисления, в соответствии с (1) требуют N 2 (N ? 1) операций умножения. Всего дополнительно потребуется
N 4 ? N 3 ? N 2 /2 ? N/2 операций.
Следовательно, основной вклад в число мультипликативных операций дает сумма (6) со слагаемыми (5), которые являются показательными функциями от N . Такие алгоритмы, у которых
сложность превышает полиномиальную, считаются очень сложными. И их стараются не применять.
Поэтому в следующем параграфе мы предлагаем улучшенный алгоритм, который имеет полиномиальную сложность.
3
Алгоритм с полиномиальной сложностью
Пусть матрица M имеет характеристический полином
h(?) = ?N + h1 ?N ?1 + . . . + hN .
(9)
Пусть m > N , рассмотрим деление с остатком полинома ?m на h(?):
?m = f (?)h(?) + r(?),
(10)
где f (?) = f0 ?m?N + f1 ?m?N ?1 + . . . + fm?N и r(?) = r1 ?N ?1 + . . . + rN .
По теореме Гамильтона-Кэли h(M ) = 0. Следовательно:
M m = r(M ).
(11)
Таким образом, задача нахождения матричной степени сводится к нахождению коэффициентов полинома r(?), степень которого меньше степени характеристического полинома. Чтобы найти эти
коэффициенты, приравняем коэффициенты при равных степенях в левой и правой части равенства (10). Получим систему, состоящую из m + 1 уравнений и содержащую m + 1 неизвестных
f0 , .., fm?N , r0 , .., rN ?1 :
?
f0 = 1 ,
?
?
?
fi = ?(h1 fi?1 + h2 fi?2 + . . . + hi f0 ) ,
i = 1, N ? 1
(12)
fi = ?(h1 fi?1 + h2 fi?2 + . . . + hN fi?N ) , i = N, m ? N + 1
?
?
?
ri = ?(hi fm?N + . . . + hN fm?2N +i ) ,
i = 1, N
Матрица коэффициентов системы имеет треугольный вид, у которой все элементы, за исключением
главной диагонали и N прилегающих к ней снизу диагоналей, равны нулю. Поэтому для решения
этой системы потребуется всего N (m ? N ) операций умножения. Учитывая то, что для вычисления
всех степеней матрицы M k , k < N , требуется N 3 (N ? 2) операций умножения, а для вычисления произведений найденных коэффициентов на матрицы в (8) потребуется N 2 (N ? 1) операций
умножения, всего получим C = N (N 3 + m) ? N 2 (N + 2) операций умножения во всем алгоритме.
Приведем соответствующий алгоритм.
Algorithm MatrixPower
Input: M , m;
Output: M m ;
N = length(M );
Coef CharP ol[] = CoeCharPol(M );
n = m ? N + 1;
f [0] = 1;
for i = 1 to n ? 1 do
k = N ? 1;
for j = i ? N to i ? 1 do
if j ? 0 then
f [i] = f [i] + CoefCharPol[k] Ч (f [j]);
k = k ? 1;
f [i] = ?f [i];
for i = 0 to N ? 1 do
k = m ? N + 1;
for j = i to N ? 1 do
k = k ? 1;
if k ? 0 then
r[i] = r[i] + CoefCharPol[j] Ч f [k];
r[i] = ?r[i];
M atrP [0] = I;
M atrP [1] = M ;
for i = 2 to N ? 1 do
M atrP [i] = M Ч M atrP [i ? 1];
for i = 0 to N ? 1 do
M atrP [i] = M atrP [i] Ч R[N ? i ? 1];
M atrP ower = M atrP ower + M atrP [i];
return M atrP ower;
4
Улучшенный алгоритм с экспоненциальной сложностью
Исследование решения системы (12) позволяет по-новому взглянуть на алгоритм с экспоненциальной сложностью.
Рассмотрим решение системы (12) методом Крамера. Определитель матрицы коэффициентов
системы будет иметь вид:
1
0
... 0
0
...
0
0 0 . . . 0 0 h1
1
... 0
0
...
0
0 0 . . . 0 0 ..
..
.
.
.
.. .. . .
. . .
.
..
..
..
..
..
.
. .. .. .
. .
hN hN ?1 . . . 1
0
...
0
0 0 . . . 0 0 0
hN
. . . h1
1
...
0
0 0 . . . 0 0 ..
..
..
..
..
.. .. . .
. . ..
..
.
.
. .. .. .
.
.
.
. .
?= .
0
0
. . . hN hN ?1 . . .
1
0 0 . . . 0 0 0
0
... 0
hN
...
h1
1 0 . . . 0 0 0
0
... 0
0
...
h2
0 1 . . . 0 0 .
.
.
.
.
.. .. . .
. . .
.
..
..
..
..
..
..
..
. .. .. . .
0
0
... 0
0
. . . hN ?1 0 0 . . . 1 0 0
0
... 0
0
...
hN
0 0 ... 0 1 Так как этот определитель равен 1, то последние N неизвестных r1 , .., rN равны соответствующим
определителям, полученным заменой одного из последних столбцов на единичный столбец. Раскладывая их по последним N столбцам, найдем ri (i=1,..,N):
h1
1
... 0
0
...
0
...
0
0 ..
..
..
..
..
..
.. ..
..
..
.
.
.
.
.
.
.
.
. .
hN hN ?1 . . . 1
0
...
0
...
0
0 0
hN
. . . h1
1
...
0
...
0
0 ..
..
..
..
..
..
.. ..
..
..
.
.
.
.
.
.
.
.
. ri = .
0
1
...
0
0 0
. . . hN hN ?1 . . .
0
0
... 0
hN
...
h1
...
0
0 .
.
.
.
.
.
.. .
.
.
..
..
..
..
..
..
..
..
..
. 0
.
.
.
h
1
0
.
.
.
0
0
.
.
.
h
N ?i
1
0
0
... 0
0
. . . hN
. . . hi+1 hi Эти определители отличаются только последней строкой: определитель ri в последней строке имеет
N ? i + 1 коэффициент характеристического полинома hN , .., hi , а их первые m ? N строки берутся
из определителя ?. Разложим ri по последней строке. Получим:
ri =
N
X
?k Zm?N +i?k ,
k=i
где ?t = (?1)t+1 ht и Zt обозначает минор, образованный строками 1, 2, .., m ? N и столбцами
1, 2, .., t ? 1, t + 1, .., m ? N + 1 определителя ?.
Отметим, что Z1 = 1 и Zk , k > 1, равен значению верхнего левого углового минора порядка k ? 1
определителя ?. Докажем по индукции, что:
Zk =
X (k1 + . . . + kN )!
k1 ! · · · kN !
k
kN
?1k1 . . . ?N
для k = 2, .., m ? N . Здесь сумма пробегает все ki , которые удовлетворяют уравнению
k1 + 2k2 + . . . + N kN = k ? 1 .
Предположим, что это утверждение верно для всех Zi , i < k . Рассматривая Zk как угловой минор
порядка k ? 1, разложим его по последней строке. Получим:
Zk =
N
X
?i Zk?i
i=1
Здесь предполагается, что Zt = 0 при t < 1.
N
X
?i Zk?i =
i=1
=
N
X
i=1
?i
X (k1 + . . . + kN )!
k?i
k1 ! · · · kN !
N X
X
(k1 + . . . + ki ? 1 + . . . + kN )!
i=1
k1 ! · · · (ki ? 1)! · · · kN !
k
=
X (k1 + . . . + kN )!
k
k1 ! · · · kN !
kN
?1k1 . . . ?N
=
kN
?1k1 . . . ?N
=
kN
?1k1 . . . ?N
Здесь мы воспользовались тождеством
N
X
(k1 + . . . + ki ? 1 + . . . + kN )!
i=1
k1 ! · · · (ki ? 1)! · · · kN !
=
(k1 + . . . + kN )!
k1 ! · · · kN !
Таким образом мы получили еще один алгоритм вычисления матричной степени
Mm =
N X
N
X
?k Zm?N +i?k M N ?i .
i=1 k=i
Этот алгоритм напоминает алгоритм (1), однако тут нет лишних слагаемых, которые имеются в
сумме (1). Соответствие этих алгоритмов устанавливается равенством Zm,p = Zm?p . Этот алгоритм
также имеет экспоненциальную сложность.
5
Численные эксперименты
Были реализованы программы для приведенных выше алгоритмов с полиномиальной и экспоненциальной сложностью вычисляющих матричную степень. Приведем результаты экспериментов по
вычислению матричных степеней для матриц пятого порядка.
Rt
Степень Rt(мс) N t(мс) t = N
Rp
N p p = Rp
t
Nt
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
N p число операций умножения в алгоритме с полиномиальной сложностью,
Rp нижняя оценка числа операций умножения в алгоритме с экспоненциальной сложностью,
N t время вычисления матричной степени в алгоритме с полиномиальной сложностью,
Rt время вычисления матричной степени при использовании алгоритма с экспоненциальной сложностью,
N Rt время вычисления матричной степени при использовании алгоритма с полиномиальной
сложностью,
N Rp оценка числа операций умножения при использовании алгоритма с полиномиальной сложностью.
Литература
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, c. 280.
3. Математическая энциклопедия/ Гл. ред. И.М. Виноградов. М.: Сов. Энциклопедия, 1979.
Т.5.
Поступила в редакцию 15 декабря 2007г.
Использование жесткого диска в матричных вычислениях
c
М.С.Зуев
Вычисления с использованием жесткого диска реализованы, например, в известном пакете численных методов POOCLAPACK (Parallel out-of-core Linear algebra package, см. [5]).
В данном пакете реализованы два способа хранения матриц на диске, один из которых повторяет
способ хранения матриц в ОЗУ по строкам или столбцам, а второй предполагает хранение матрицы
по блокам. Во втором случае требуется использование блочных матричных алгоритмов.
В современной литературе блочные и блочно-рекурсивные матричные методы рассматриваются
как одни из самых эффективных. В отличие от стандартных строчно-ориентированных методов,
такие методы более эффективно используют процессорный кэш, и поэтому позволяют достигать
производительности процессора, близкой к пиковой (см. [2-4]). Например, в работе [2] указывается
достижение 92 % производительности процессора IBM Power 3 с блочной процедурой разложения
Холецкого и 76 % с соответствующей поэлементной процедурой.
Если процедура использует жесткий диск, то она может эффективно использовать оперативную
память. Считывание блоков из файла целиком и возможность варьировать их размер позволяет
определить баланс между экономией оперативной памяти и минимизацией количества обращений
к диску. При считывании блока из файла требуется одно обращение к таблице размещения файлов
(суперблоку) и одно обращение к файлу. Если бы блоки хранились в виртуальной памяти, то для
их считывания требовалось бы большее количество обращений к диску.
Суммируя все сказанное выше, можно предположить, что блочные матричные алгоритмы с
использованием жесткого диска могут быть эффективными.
Для исследования эффективности таких алгоритмов были проведены вычислительные эксперименты на однопроцессорной машине. Сравнивались алгоритмы умножения и обращения матриц
(одностороннее обращение, см. [1]) размера 4096 Ч 4096 в конечном поле по 28-битному простому
модулю. Как и в POOCLAPACK, матрица разбивается на блоки фиксированного размера, которые
записываются в отдельные файлы. Матрице сопоставляется папка на диске. Если разбить матрицу
на четыре равных блока, то возможно их хранение в отдельных папках.
Для определения лучшего алгоритма умножения были рассмотрены четыре алгоритма: стандартный блочный алгоритм (алгоритм 0), стандартный блочный алгоритм с записью блоков на
Документ
Категория
Без категории
Просмотров
4
Размер файла
189 Кб
Теги
степени, вычисления, матричный, функции, матричное
1/--страниц
Пожаловаться на содержимое документа