Блок факторизации ЛПНП для эрмитовских неопределенных матриц
L = ldl(A)
[L,D] = ldl(A)
[L,D,P] = ldl(A)
[L,D,p] = ldl(A,'vector')
[U,D,P] = ldl(A,'upper')
[U,D,p] = ldl(A,'upper','vector')
[L,D,P,S] = ldl(A)
[L,D,P,S] = LDL(A,THRESH)
[U,D,p,S] = LDL(A,THRESH,'upper','vector')
L = ldl(A) возвращает только переставленную нижнюю треугольную матрицу L как в форме с двумя выходами. Информация о перестановках теряется, как и коэффициент диагонали блока D. По умолчанию ldl ссылается только на диагональ и нижний треугольник Aи предполагает, что верхний треугольник является комплексно сопряженным транспонированием нижнего треугольника. Поэтому [L,D,P] = ldl(TRIL(A)) и [L,D,P] = ldl(A)оба возвращают одинаковые коэффициенты. Обратите внимание, что этот синтаксис недопустим для разреженного A.
[L,D] = ldl(A) сохраняет блок-диагональную матрицу D и переставленная нижняя треугольная матрица в L такой, что A = L*D*L'. Блок-диагональная матрица D имеет блоки 1 на 1 и 2 на 2 по диагонали. Обратите внимание, что этот синтаксис недопустим для разреженного A.
[L,D,P] = ldl(A) возвращает единицу нижней треугольной матрицы L, блок-диагональ Dи матрица перестановок P такой, что P'*A*P = L*D*L'. Это эквивалентно [L,D,P] = ldl(A,'matrix').
[L,D,p] = ldl(A,'vector') возвращает информацию о перестановке в виде вектора, pвместо матрицы. p output - вектор строки, такой, что A(p,p) = L*D*L'.
[U,D,P] = ldl(A,'upper') ссылается только на диагональ и верхний треугольник A и предполагает, что нижний треугольник является комплексно сопряженным транспонированием верхнего треугольника. Этот синтаксис возвращает верхнюю треугольную матрицу единиц измерения U такой, что P'*A*P = U'*D*U (при условии, что A эрмитова, а не только верхняя треугольная). Аналогично, [L,D,P] = ldl(A,'lower') задает поведение по умолчанию.
[U,D,p] = ldl(A,'upper','vector') возвращает информацию о перестановке в виде вектора, p, как и [L,D,p] = ldl(A,'lower','vector'). A должна быть полной матрицей.
[L,D,P,S] = ldl(A) возвращает единицу нижней треугольной матрицы L, блок-диагональ D, матрица перестановок Pи матрица масштабирования S такой, что P'*S*A*S*P = L*D*L'. Этот синтаксис доступен только для вещественных разреженных матриц, и только нижний треугольник A используется ссылка.
[L,D,P,S] = LDL(A,THRESH) использование THRESH в качестве проворотного допуска в алгоритме. THRESH должен быть двойным скаляром, лежащим в интервале [0, 0.5]. Значение по умолчанию для THRESH является 0.01. Использование меньших значений THRESH может дать более быстрое время факторизации и меньшее количество записей, но может также привести к менее стабильной факторизации. Этот синтаксис доступен только для вещественных разреженных матриц.
[U,D,p,S] = LDL(A,THRESH,'upper','vector') задает допуск поворота и возвращает верхний треугольник U и вектор перестановки p как описано выше.
Эти примеры иллюстрируют использование различных форм ldl функции, включая форму с одним, двумя и тремя выходами, и использование vector и upper варианты. Рассматриваются следующие темы:
Перед выполнением любого из этих примеров необходимо создать следующие положительные определенные и неопределенные эрмитовы матрицы:
A = full(delsq(numgrid('L', 10)));
B = gallery('uniformdata',10,0);
M = [eye(10) B; B' zeros(10)]; Структура M здесь очень часто возникают проблемы оптимизации и потока жидкости, и M фактически является неопределенным. Обратите внимание, что положительная определенная матрица A должен быть заполнен, как ldl не принимает разреженные аргументы.
Форма с двумя выходами ldl прибыль L и D такой, что A-(L*D*L') мал, L перестановка единицы нижнего треугольника, и D является диагональю 2 на 2 блока. Отметим также, что, потому что A является положительным определенным, диагональ D является положительным:
[LA,DA] = ldl(A); fprintf(1, ... 'The factorization error ||A - LA*DA*LA''|| is %g\n', ... norm(A - LA*DA*LA')); neginds = find(diag(DA) < 0)
Данный a b, решить Ax=b использование LA, DA:
bA = sum(A,2); x = LA'\(DA\(LA\bA)); fprintf(... 'The absolute error norm ||x - ones(size(bA))|| is %g\n', ... norm(x - ones(size(bA))));
Три выходные формы также возвращают матрицу перестановки, так что L на самом деле единица нижняя треугольная:
[Lm, Dm, Pm] = ldl(M); fprintf(1, ... 'The error norm ||Pm''*M*Pm - Lm*Dm*Lm''|| is %g\n', ... norm(Pm'*M*Pm - Lm*Dm*Lm')); fprintf(1, ... 'The difference between Lm and tril(Lm) is %g\n', ... norm(Lm - tril(Lm)));
Данный b, решить Mx=b использование Lm, Dm, и Pm:
bM = sum(M,2); x = Pm*(Lm'\(Dm\(Lm\(Pm'*bM)))); fprintf(... 'The absolute error norm ||x - ones(size(b))|| is %g\n', ... norm(x - ones(size(bM))));
D представляет собой блок-диагональную матрицу с блоками 1 на 1 и 2 на 2. Это делает его частным случаем тридиагональной матрицы. Когда входная матрица является положительной определенной, D почти всегда диагональна (в зависимости от того, насколько определенна матрица). Однако, когда матрица является неопределенной, D может быть диагональной или может выражать структуру блока. Например, с A как указано выше, DA диагональ. Но если вы сдвинетесь A немного, в итоге получается неопределенная матрица, а затем можно вычислить D имеет структуру блока.
figure; spy(DA); title('Structure of D from ldl(A)');
[Las, Das] = ldl(A - 4*eye(size(A)));
figure; spy(Das);
title('Structure of D from ldl(A - 4*eye(size(A)))');
Как и lu функция, ldl принимает аргумент, который определяет, возвращает ли функция вектор перестановки или матрицу перестановки. ldl возвращает последний по умолчанию. При выборе 'vector', функция выполняется быстрее и использует меньше памяти. По этой причине, указав 'vector' рекомендуется вариант. Другое дело, что индексирование обычно происходит быстрее, чем умножение для этого вида операций:
[Lm, Dm, pm] = ldl(M, 'vector');
fprintf(1, 'The error norm ||M(pm,pm) - Lm*Dm*Lm''|| is %g\n', ...
norm(M(pm,pm) - Lm*Dm*Lm'));
% Solve a system with this kind of factorization.
clear x;
x(pm,:) = Lm'\(Dm\(Lm\(bM(pm,:))));
fprintf('The absolute error norm ||x - ones(size(b))|| is %g\n', ...
norm(x - ones(size(bM))));Как и chol функция, ldl принимает аргумент, который определяет, на какой треугольник входной матрицы ссылается, а также ldl возвращает более низкий (L) или верхний (L') треугольный коэффициент. Для плотных матриц нет реальной экономии при использовании верхнего треугольного варианта вместо нижнего треугольного варианта:
Ml = tril(M); [Lml, Dml, Pml] = ldl(Ml, 'lower'); % 'lower' is default behavior. fprintf(1, ... 'The difference between Lml and Lm is %g\n', norm(Lml - Lm)); [Umu, Dmu, pmu] = ldl(triu(M), 'upper', 'vector'); fprintf(1, ... 'The difference between Umu and Lm'' is %g\n', norm(Umu - Lm')); % Solve a system using this factorization. clear x; x(pm,:) = Umu\(Dmu\(Umu'\(bM(pmu,:)))); fprintf(... 'The absolute error norm ||x - ones(size(b))|| is %g\n', ... norm(x - ones(size(bM))));
При указании обоих 'upper' и 'vector' опции, 'upper' должен предшествовать 'vector' в списке аргументов.
При использовании linsolve функция, вы можете испытать лучшую производительность, используя знания, что система имеет симметричную матрицу. Матрицы, используемые в примерах выше, немного малы, чтобы видеть это, поэтому для этого примера генерируется матрица большего размера. Матрица здесь симметричная положительная определенная, и ниже мы увидим, что с каждым битом знаний о матрице, есть соответствующее ускорение. То есть симметричный решатель быстрее, чем общий решатель, в то время как симметричный положительный определенный решатель быстрее, чем симметричный решатель:
Abig = full(delsq(numgrid('L', 30)));
bbig = sum(Abig, 2);
LSopts.POSDEF = false;
LSopts.SYM = false;
tic; linsolve(Abig, bbig, LSopts); toc;
LSopts.SYM = true;
tic; linsolve(Abig, bbig, LSopts); toc;
LSopts.POSDEF = true;
tic; linsolve(Abig, bbig, LSopts); toc;[1] Эшкрафт, К., Р. Г. Граймс и Дж. Г. Льюис. «Точные решатели симметричных неопределенных линейных уравнений». SIAM J. Matrix Anal. Appl. Том 20. Номер 2, 1998, стр. 513-561.
[2] Дафф, И. С. «MA57 - Новый код решения разреженных симметричных определённых и неопределённых систем». Технический отчет RAL-TR-2002-024, Лаборатория Резерфорда Эпплтона, 2002 год.