Блокируйте LDL-разложение для Эрмитовых неопределенных матриц
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 как в 2D выходной форме. Информация о перестановке потеряна, как диагональ блока факторный 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 имеет блоки 2 на 2 и 1 на 1 на своей диагонали. Отметьте, этот синтаксис не допустим для разреженного 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 вывод является вектором - строкой, таким образом что 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 ссылаются. ldl использует MA57 для разреженного действительного симметричного A.
[L,D,P,S] = LDL(A,THRESH) использование THRESH как допуск центра в MA57. THRESH должен быть двойным скаляром, лежащим в интервале [0, 0.5]. Значением по умолчанию для THRESH является 0.01. Используя меньшие значения THRESH может дать более быстрые времена факторизации и меньше записей, но может также привести к менее стабильной факторизации. Этот синтаксис доступен только для действительных разреженных матриц.
[U,D,p,S] = LDL(A,THRESH,'upper','vector') устанавливает допуск центра и возвращает верхний треугольный U и вектор перестановки p, как описано выше.
Эти примеры иллюстрируют использование различных форм функции ldl, включая одну - 2D, и формы с тремя выводами и использования опций upper и vector. Затронутые темы:
Прежде, чем запустить любой из этих примеров, необходимо будет сгенерировать следующие положительные определенные и неопределенные Эрмитовы матрицы:
A = full(delsq(numgrid('L', 10)));
B = gallery('uniformdata',10,0);
M = [eye(10) B; B' zeros(10)]; Структура M здесь очень распространена в оптимизации и проблемах потока жидкости, и M на самом деле неопределенен. Обратите внимание на то, что положительный определенный матричный A должен быть полным, когда ldl не принимает разреженные аргументы.
2D выходная форма ldl возвращает L и D, таким образом, что A-(L*D*L') является маленьким, L является переставленным нижним треугольным модулем, и D является диагональю блока 2-by-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;ldl использует стандартные программы MA57 в Harwell Subroutine Library (HSL) для действительных разреженных матриц.
[1] Ashcraft, C., Р.Г. Граймс и Дж.Г. Льюис. “Точные Симметричные Неопределенные Решатели Линейных уравнений”. Издание 20 SIAM J. Matrix Anal. Appl. Номер 2, 1998, стр 513–561.
[2] Подновите, я. S. "MA57 — новый код для решения разреженных симметричных определенных и неопределенных систем". Технический отчет RAL TR-2002-024, Лаборатория Резерфорда Эпплтона, 2002.