ldl

Блокируйте 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 не принимает разреженные аргументы.

Пример 1 - двух-Выходная Форма 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)))); 

Пример 2 - Три Выходных Формы ldl

Три выходных формы возвращают матрицу перестановок также, так, чтобы 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)))); 

Пример 3 - структура D

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)))');

Пример 4 - Используя 'векторную' опцию

Как 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))));

Пример 5 - Используя 'верхнюю' опцию

Как 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' в списке аргументов.

Пример 6 - linsolve и Эрмитов неопределенный решатель

При использовании 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.

Расширенные возможности

Смотрите также

| |