Блокируйте 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 как в форме с двумя выходами. Информация о сочетании теряется, как и блочный диагональный коэффициент 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, вместо матрицы. The 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 функция, включая одно-, двух- и трех- выходную форму, и использование 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] Ashcraft, C., R.G. Grimes, and J.G. Lewis. «Точные симметричные неопределенные линейные Уравнения Решателей». SIAM J. Matrix Anal. Appl. vol. 20. № 2, 1998, с. 513-561.
[2] Duff, I. S «. MA57 - Новый код для решения разреженных симметричных определенных и неопределенных систем». Технический отчет RAL-TR-2002-024, Лаборатория Резерфорда Эпплтона, 2002.