Метод окаймления. Популярный алгоритм обращения симметричных матриц методом окаймления состоит в следующем.

A =
a
b
b'
c
A-1 =
a-1-pp'/d
p/d
p'/d
-1/d

где p=a-1b, d=b'p-c.


A=[ 1 2 3; 2 1 5; 3 5 1], 
B=invs(A), I={A*B}, mesh(I), B=?,

function: invs(A),
var i n b c d p, n=rowm(A), 
B=ones(1), B[0][0]=1.0/A[0][0], 
for i=1:n, 
c=A[i][i], b=block(A[i] 0 i-1),  
p={B*b},   d={b'*p}, d=d-c, c={p/d},  
B=res(B-p*c), p={[p -1]}, p={p/d}, 
B={[B p]}, tr(B), B={[B p]}, 
end, return B,
end,

Алгоритм Гревилля. Сходные идеи используют для обобщенного обращения (псевдоинверсии) матриц, наращиванием их боковой каймой.

A =
a
b
A+ =
a+-pc
c

где p=a+b, на старте d=a (если a=0, то a+=0), далее d=b-ap, c=d'/d'd и если d=0, то c=p'a+/(1+p'p).


A=[ 1 1 0; 0 1 0; 0 0 1], B=pinvs(A), I={A*B}, mesh(I), 

function: pinvs(A),
var i n m b c d p, n=colm(A), 
m=rowm(A), d={A[0]}, c={d'*d}, 
if c<>0, B={d/c}, else, B=zero(d), end, 
for i=1:n, b={A[i]}, p={B'*b},
c=blocks(A 0 m 0 i-1), 
d={c*p}, d={b-d},  c={d'*d}, 
if c=0,  c={p'*p}, c=1+c, d={B*p}, end,
c={d/c}, p={c*p}, B={B-p}, B={[B c]}, 
end, tr(B), 
return B,
end,

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

В свое время это породило обширную литературу с аналитическими выкладками решения разнообразных прикладных задач, которые мало оправданы с точки зрения возможностей реализации. Этот исторический аспект до сих пор не преодолен, порождая последующие цитирования.

Алгоритм Гревилля, приведенный выше, интересен тем, что на промежуточных этапах вычислений возникает вектор невязки текущего столбца матрицы A и линейной комбинации предыдущих столбцов d=b-ap, взятый в максимальном приближении выбором соответствующих коэффициентов p.

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

При ортогонализации системы векторов методом Грама-Шмидта такая стратегия отвечала бы выбору в первую очередь наиболее растворенных векторов.



Rambler's Top100