Блокируйте 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 ссылается.
[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 функция, включая одну - 2D, и форма с тремя выходами и использование vector и upper опции. Затронутые темы:
Прежде, чем запустить любой из этих примеров, необходимо будет сгенерировать следующие положительные определенные и неопределенные Эрмитовы матрицы:
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 использование LAdA :
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] Ashcraft, C., Р.Г. Граймс и Дж.Г. Льюис. “Точные Симметричные Неопределенные Решатели Линейных уравнений”. Издание 20 SIAM J. Matrix Anal. Appl. Номер 2, 1998, стр 513–561.
[2] Подновите, я. S. "MA57 — новый код для решения разреженных симметричных определенных и неопределенных систем". Технический отчет RAL TR-2002-024, Лаборатория Резерфорда Эпплтона, 2002.