Тройная факторизация
5.3 Тройная факторизация.
Часто вместо разложения Холецкого используют разложение 
, в котором 
‑ диагональная матрица, а треугольные матрицы 
и 
имеют единичную диагональ. В матричном виде тройная факторизация выглядит так
Чтобы найти коэффициенты матриц 
вначале вычислим произведение 
. Получим нижнюю треугольную матрицу 
с элементами 
. Произведение 
дает
Отсюда при 
находим
, (6)
а при 
имеем
. (7)
Рекомендуемые материалы
Вычисления тройной факторизации по формулам (6), (7) проводятся в том же порядке, как и в разложении Холецкого. Вычисления ведутся по строкам 
(сверху-вниз). Для произвольной строки 
вначале по формуле (6) находим диагональ 
, вспоминаем, что 
, а затем для 
по формуле (7) вычисляем 
. Заметим, что при тройной факторизации не надо вычислять квадратные корни.
После того, как множители 
определены, вычисляем решение в два этапа. Сначала решаем задачу 
, или 
с нижней треугольной матрицей. Пользуясь формулой (3), получаем
Затем решаем задачу 
с верхней треугольной матрицей:
Пример 1. Решение с помощью разложения Холецкого.
Входные данные. matrix A. Количество ненулевых в треугольной части 7
4 1 2 0.5 2
1 0.5 0 0 0
2 0 3 0 0
0.5 0 0 0.625 0
2 0 0 0 16
vector f
7 3 7 -4 -4
Результат разложения – матрица . Ненулевых нет, заполнение 6.
2.00 0.50 1.00 0.25 1.00
0.00 0.50 -1.00 -0.25 -1.00
0.00 0.00 1.00 -0.50 -2.00
0.00 0.00 0.00 0.50 -3.00
0.00 0.00 0.00 0.00 1.00
Решение задачи 
. Вектор 
(3.50 2.50 6.00 -2.50 -0.50).
Решение задачи . Вектор 
(3.50 2.50 6.00 -2.50 -0.50).
------------------------------------
Код решения задачи.
procedure TMainForm.RunMClick(Sender: TObject);
begin
ReadAF('AandF.txt',N, A,f);
CholezkyFactor(N,A, U);
SaveMatrix('U.txt', N, U);
SolveL(N,U,f, y);
SaveVector('y.txt', N, y);
SolveU(N,U,y, x);
SaveVector('x.txt', N, x);
Close;
end; //Run
Коды процедуры разложения Холецкого.
Procedure CholezkyFactor(N:integer; A:matrix; var U:matrix);
var
i,j,k:integer;
s,d:float;
begin
For i:=1 to N do
begin
For j:=1 to i-1 do U[i,j]:=0;
s:=0;
for k:=1 to i-1 do s:=s + sqr(U[k,i]);
U[i,i]:=sqrt(A[i,i] - s);
For j:=i+1 to N do
begin
s:=0;
for k:=1 to i-1 do s:=s + U[k,i]*U[k,j];
U[i,j]:=(A[i,j] - s)/U[i,i];
end; //j
end; //i
end; //CholezkyFactor
Код решения систем с треугольными матрицами.
Procedure SolveL(N:integer;U:matrix;f:vector; var x:vector);
{L=UT}
var
i,j:integer;
s:float;
begin
For i:=1 to N do
begin
s:=0;
For j:=1 to i-1 do s:=s + U[j,i]*x[j];
x[i]:=(f[i]-s)/U[i,i];
end; //i
end; //SolveL
Procedure SolveU(N:integer;U:matrix;f:vector; var x:vector);
var
i,j:integer;
s:float;
begin
For i:=N downto 1 do
begin
s:=0;
For j:=i+1 to N do s:=s + U[i,j]*x[j];
x[i]:=(f[i]-s)/U[i,i];
end; //i
end; //SolveU
На рассмотренном примере мы видели, что при факторизации Холецкого происходит полное заполнение треугольных множителей. Этот нежелательный эффект можно уменьшить, перенумеровав переменные. Если сделать перенумерацию по правилу
то последнее уравнений станет первым, второе снизу – вторым и т.д., наконец бывшее первое уравнение станет последним. В результате получим эквивалентную систему уравнений 
, где
matrix 
16 0 0 0 2
0 0.625 0 0 0.5
0 0 3 0 2
0 0 0 0.5 1
2 0.5 2 1 4
vector 
=(-4 -4 7 3 7).
В результате разложения получается матрица 
4.00 0.00 0.00 0.00 0.50
0.00 0.79 0.00 0.00 0.63
0.00 0.00 1.73 0.00 1.15
0.00 0.00 0.00 0.71 1.41
0.00 0.00 0.00 0.00 0.13
У нее заполнение равно нулю, т.е. новых ненулевых элементов не появилось. Решение задачи, разумеется, то же самое, с точностью до наоборот:
( -0.50 -8.00 1.00 2.00 2.00).
В реальной жизни такую перенумерацию придумать не удается, однако существуют алгоритмы, уменьшающие заполнение.
Пример 2. Тройная факторизация.
Матрица возникла в ходе конечноэлементной сборки на нерегулярной сетке, построенной с помощью разбиения Делоне (алгоритм изучали ранее).
Портрет матрицы показан на рисунке. Черным отмечены ненулевые элементы (их 35 985 из 28 206 721 в матрице ), а красным – заполнение (1 078 803).
Рекомендуем посмотреть лекцию "15. 2. Процессуальные характеристики конфликта".
После перенумерации по алгоритму Катхилла-МакКи, уменьшающему ширину ленты, получается матрица с портретом
У этой матрицы, разумеется, столько же ненулевых элементов, однако заполнение при факторизации «всего» 349 611, т.е. в 3 раза меньше, чем у исходной. Однако это не самый оптимальный алгоритм перенумерации. Например, алгоритм минимальной степени, локально минимизирующий именно заполнение, а не ширину ленты, дает матрицу с неленточным портретом, однако заполнение при факторизации такой матрицы равно 87 923, т.е. в 4 раза меньше, чем в алгоритме Катхилла-МакКи. На больших задачах это преимущество еще более возрастает.
Задание. Написать программу тройной факторизации, проверить заполнение на примере матрицы, для которой проводилось разложение Холецкого, решить задачу, убедиться, что решения совпадают. .
Заключение по прямым методам решения сеточных СЛАУ. Прямые методы, основанные на тругольных разложениях исходной матрицы, с успехом используются в ВГ для задач малой и средней размерности. Достоинство этих методов – возможность получения точного (ошибки округления могут его искажать) решения за конечное число операций. Особенно привлекательны прямые методы в нестационарных задачах с оператором, не зависящим от времени. Тогда факторизацию матрицы можно провести один раз, а затем на каждом временном слое использовать треугольные множители для вычисления решения по меняющейся правой части. В этом случае единовременные накладные расходы на факторизацию становятся малыми по сравнению с общим временем счета. В таких задачах прямые методы выглядят даже предпочтительнее многих итерационных. Хотя и они сходятся довольно быстро, поскольку начальное приближение берется с предыдущего слоя, и при гладких зависимостях правой части от времени и достаточно малых это начальное приближение оказывается очень хорошим. Главный недостаток факторизованных методов – заполнение, которое 1) увеличивает до непомерных размеров объем памяти для хранения множителей; 2) увеличивает количество операций при умножении заполненного треугольника на вектор. Преодоление этих трудностей достигается за счет 1) перенумерации узлов сетки (= переупорядочивание неизвестных); 2) применение схем хранения разреженных векторов и матриц и применение таких методов их обработки, которые исключают действия с нулями.
