Документация

Вычислительная сложность

Вычислительная сложность разреженных операций пропорциональна nnz, количеству ненулевых элементов в матрице. Вычислительная сложность также зависит линейно от размера строки m и размер столбца n матрицы, но независима от продукта m*n, общее количество нулевых и ненулевых элементов.

Сложность довольно сложных операций, таких как решение разреженных линейных уравнений, включает факторы как упорядоченное расположение и временная замена, которые обсуждены в предыдущем разделе. В целом, однако, машинное время, требуемое для операции разреженной матрицы, пропорционально количеству арифметических операций на ненулевых количествах.

Алгоритмические детали

Разреженные матрицы распространяют посредством вычислений согласно этим правилам:

Переупорядочение для разреженности

Переупорядочение столбцов матрицы может часто делать свой LU или факторы QR более разреженными. Переупорядочение строк и столбцов может часто делать свои Факторы Холесского более разреженными. Самое простое такое переупорядочение должно отсортировать столбцы по ненулевому количеству. Это иногда - хорошее переупорядочение для матриц с очень неправильными структурами, особенно если существует большое изменение в ненулевых количествах строк или столбцов.

colperm вычисляет перестановку, которая заказывает столбцы матрицы количеством ненулей в каждом столбце от самого маленького до самого большого.

Переупорядочение для уменьшения пропускной способности

Упорядоченное расположение обратного алгоритма Катхилла-Макки предназначается, чтобы уменьшать профиль или пропускную способность матрицы. Это, как гарантируют, не найдет самую маленькую пропускную способность, но это обычно делает. Функция symrcm на самом деле управляет на ненулевой структуре симметрической матрицы A + A', но результат также полезен для несимметричных матриц. Это упорядоченное расположение полезно для матриц, которые прибывают из одномерных проблем или проблем, которые находятся в некотором смысле, длинном и тонком.

Приближение минимального упорядоченного расположения градуса

Степень узла в графике является количеством связей с тем узлом. Это совпадает с количеством недиагональных ненулевых элементов в соответствующей строке матрицы смежности. Аппроксимированный минимальный алгоритм степени генерирует упорядоченное расположение на основе того, как эти степени изменены во время Исключения Гаусса или факторизации Холесского. Это - сложный и мощный алгоритм, который обычно приводит к более разреженным факторам, чем большинство других упорядоченных расположений, включая количество столбцов и обратный алгоритм Катхилла-Макки. Поскольку отслеживание степени каждого узла очень длительно, аппроксимированный минимальный алгоритм степени использует приближение для степени, а не точной степени.

Эти функции MATLAB^® реализуют аппроксимированный минимальный алгоритм степени:

Смотрите Переупорядочение и Факторизацию для примера с помощью symamd.

Можно изменить различные параметры, сопоставленные с деталями алгоритмов с помощью функции spparms.

Для получения дополнительной информации на алгоритмах, используемых colamd и symamd, см. [5]. Аппроксимированная степень использование алгоритмов основана [1].

Вложенное упорядоченное расположение рассечения

Как аппроксимированное минимальное упорядоченное расположение степени, вложенный алгоритм упорядоченного расположения рассечения, реализованный функцией dissect, переупорядочивает строки матрицы и столбцы, полагая, что матрица матрица смежности графика. Алгоритм уменьшает проблему вниз до намного меньшей шкалы путем сворачивания вместе пар вершин в графике. После переупорядочения маленького графика алгоритм затем применяет проекцию, и улучшение продвигается, чтобы расширить график назад до первоначального размера.

Вложенный алгоритм рассечения производит высококачественные переупорядочения и выполняет особенно хорошо с матрицами конечного элемента по сравнению с другими методами переупорядочения. Для получения дополнительной информации о вложенном алгоритме упорядоченного расположения рассечения, см. [7].

LU-факторизация

Если S является разреженной матрицей, следующая команда возвращает три разреженных матрицы L, U и P, таким образом что P*S = L*U.

[L,U,P] = lu(S);

lu получает факторы Исключением Гаусса с частичным поворотом. Матрица перестановок P имеет только n ненулевые элементы. Как с плотными матрицами, оператор [L,U] = lu(S) возвращает переставленный модуль нижняя треугольная матрица и верхняя треугольная матрица, продуктом которой является S. Отдельно, lu(S) возвращает L и U в одной матрице без информации о центре.

Синтаксис с тремя выводами [L,U,P] = lu(S) выбирает P через числовой частичный поворот, но не вертится, чтобы улучшить разреженность в факторах LU. С другой стороны, синтаксис с четырьмя выводами [L,U,P,Q] = lu(S) выбирает P через порог частичный поворот и выбирает P и Q, чтобы улучшить разреженность в факторах LU.

Можно управлять поворотом в использовании разреженных матриц

lu(S,thresh)

где thresh является порогом центра в [0,1]. Поворот происходит, когда диагональный элемент в столбце имеет значение меньше, чем времена thresh значение любого поддиагонального элемента в том столбце. thresh = 0 обеспечивает диагональный поворот.   thresh = 1 является значением по умолчанию. (Значением по умолчанию для thresh является 0.1 для синтаксиса с четырьмя выводами).

Когда вы вызываете lu с тремя или меньше выходными параметрами, MATLAB автоматически выделяет память, необходимую, чтобы содержать разреженный L и факторы U во время факторизации. За исключением синтаксиса с четырьмя выводами, MATLAB не использует символьной предварительной факторизации LU, чтобы определить требования к памяти и настроить структуры данных заранее.

Переупорядочение и Факторизация.  Этот пример показывает эффекты переупорядочения и факторизации на разреженных матрицах.

Если вы получаете хорошую перестановку столбца p, который уменьшает временную замену, возможно, от symrcm или colamd, то вычисление lu(S(:,p)) занимает меньше времени и устройства хранения данных, чем вычисление lu(S).

Создайте разреженную матрицу с помощью Маркерного примера шара.

B = bucky;

B имеет точно три ненулевых элемента в каждой строке и столбце.

Создайте две перестановки, r и m с помощью symrcm и symamd соответственно.

r = symrcm(B);
m = symamd(B);

Эти две перестановки являются симметричным упорядоченным расположением обратного алгоритма Катхилла-Макки и симметричным аппроксимированным минимальным упорядоченным расположением степени.

Создайте графики шпиона показать три матрицы смежности Маркерного графика Шара с этими тремя различными нумерациями. Локальная, основанная на пятиугольнике структура исходной нумерации не очевидна в других.

figure
subplot(1,3,1)
spy(B)
title('Original')

subplot(1,3,2)
spy(B(r,r))
title('Reverse Cuthill-McKee')

subplot(1,3,3)
spy(B(m,m))
title('Min Degree')

Упорядоченное расположение обратного алгоритма Катхилла-Макки, r, уменьшает пропускную способность и концентрирует все ненулевые элементы около диагонали. Аппроксимированное минимальное упорядоченное расположение степени, m, производит подобную фракталу структуру с большими блоками нулей.

Чтобы видеть временную замену, сгенерированную в LU-факторизации Маркерного шара, используйте speye, разреженную единичную матрицу, чтобы вставить-3s на диагонали B.

B = B - 3*speye(size(B));

Поскольку каждая сумма строки является теперь нулем, этот новый B на самом деле сингулярен, но это все еще поучительно, чтобы вычислить его LU-факторизацию. Когда названо только одним выходным аргументом, lu возвращает два треугольных фактора, L и U, в одной разреженной матрице. Количество ненулей в той матрице является мерой времени и устройства хранения данных, требуемого решить линейные системы, включающие B.

Вот ненулевые счета для этих трех рассматриваемых перестановок.

Даже при том, что это - небольшой пример, результаты типичны. Исходная схема нумерации приводит к большей части временной замены. Временная замена для упорядоченного расположения обратного алгоритма Катхилла-Макки сконцентрирована в полосе, но это почти столь же обширно как первые два упорядоченных расположения. Для аппроксимированного минимального упорядоченного расположения степени относительно большие блоки нулей сохраняются во время устранения, и сумма временной замены является значительно меньше, чем сгенерированный другими упорядоченными расположениями.

spy строит ниже отражения характеристики каждого переупорядочения.

figure
subplot(1,3,1)
spy(lu(B))
title('Original')

subplot(1,3,2)
spy(lu(B(r,r)))
title('Reverse Cuthill-McKee')

subplot(1,3,3)
spy(lu(B(m,m)))
title('Min Degree')

Факторизация Холесского

Если S является симметричным (или Эрмитов), положительная определенная, разреженная матрица, оператор ниже возвратов разреженная, верхняя треугольная матрица R так, чтобы R'*R = S.

R = chol(S)

chol автоматически не вертится для разреженности, но можно вычислить аппроксимированную минимальную степень и ограничивающие перестановки профиля для использования с chol(S(p,p)).

Поскольку алгоритм Холесского не использует поворот для разреженности и не требует поворота для числовой устойчивости, chol делает быстрое вычисление требуемого объема памяти и выделяет всю память в начале факторизации. Можно использовать symbfact, который использует тот же алгоритм в качестве chol, чтобы вычислить, сколько памяти выделяется.

QR-факторизация

MATLAB вычисляет полную QR-факторизацию разреженной матрицы S с

 [Q,R] = qr(S)

[Q,R,E] = qr(S)

но это часто непрактично. Унитарной матрице Q часто не удается иметь высокий процент нулевых элементов. Более практическая альтернатива, иногда известная как “QR-факторизацию Q-less”, доступна.

С одним разреженным входным параметром и одним выходным аргументом

R = qr(S)

возвращает только верхний треугольный фрагмент QR-факторизации. Матричный R обеспечивает факторизацию Холесского для матрицы, сопоставленной с нормальными уравнениями:

R'*R = S'*S

Однако потери числовой информации, свойственной от вычисления S'*S, избегают.

С двумя входными параметрами, имеющими одинаковое число строк, и два выходных аргумента, оператор

[C,R] = qr(S,B)

применяет ортогональные преобразования к B, производя C = Q'*B, не вычисляя Q.

QR-факторизация Q-less позволяет решение разреженных проблем наименьших квадратов

с двумя шагами:

[c,R] = qr(A,b);
x = R\c

Если A является разреженным, но не квадратным, MATLAB использует эти шаги для линейного уравнения, решая оператор наклонной черты влево:

x = A\b

Или, можно сделать факторизацию сами и исследовать R на дефицит ранга.

Также возможно решить последовательность линейных систем наименьших квадратов с различными правыми сторонами, b, которые не обязательно известны, когда R = qr(A) вычисляется. Подход решает “полунормальные уравнения R'*R*x = A'*b с

x = R\(R'\(A'*b))

и затем использует один шаг итеративного улучшения, чтобы уменьшать, округляют ошибку:

r = b - A*x;
e = R\(R'\(A'*r));
x = x + e

Неполные факторизации

ilu и функции ichol обеспечивают аппроксимированные, неполные факторизации, которые полезны как предварительные формирователи для разреженных итерационных методов.

Функция ilu производит три incomplete lower-upper (ILU) факторизации: zero-fill ILU (ILU(0)), версия Crout ILU (ILUC(tau)) и ILU с порогом, понижающимся и вертящимся (ILUTP(tau)). ILU(0) никогда не вертится, и получившиеся факторы только имеют ненули в положениях, где входная матрица имела ненули. И ILUC(tau) и ILUTP(tau), однако, делают пороговое отбрасывание с пользовательским допуском отбрасывания tau.

Например:

A = gallery('neumann', 1600) + speye(1600);

ans =

        7840

nnz(lu(A))

ans =

      126478

показывает, что A имеет 7 840 ненулей, и его полная LU-факторизация имеет 126 478 ненулей. С другой стороны, следующий код показывает различному ILU выходные параметры:

[L,U] = ilu(A);
nnz(L)+nnz(U)-size(A,1)

ans =

        7840

norm(A-(L*U).*spones(A),'fro')./norm(A,'fro')

ans =

   4.8874e-17

opts.type = 'ilutp';
opts.droptol = 1e-4;
[L,U,P] = ilu(A, opts);
nnz(L)+nnz(U)-size(A,1)

ans =

       31147

norm(P*A - L*U,'fro')./norm(A,'fro')

ans =

   9.9224e-05

opts.type = 'crout';
[L,U,P] = ilu(A, opts);
nnz(L)+nnz(U)-size(A,1)

ans =

       31083

norm(P*A-L*U,'fro')./norm(A,'fro')

ans =

   9.7344e-05

Эти вычисления показывают, что заполнить нулями факторы имеют 7 840 ненулей, факторы ILUTP(1e-4) имеют 31 147 ненулей, и факторы ILUC(1e-4) имеют 31 083 ненуля. Кроме того, относительная погрешность продукта заполнить нулями факторов является по существу нулем на шаблоне A. Наконец, относительная погрешность в факторизациях, произведенных с пороговым отбрасыванием, находится на том же порядке допуска отбрасывания, несмотря на то, что это, как гарантируют, не произойдет. Смотрите страницу с описанием ilu для большего количества опций и деталей.

Функция ichol обеспечивает zero-fill incomplete Cholesky factorizations (IC(0)), а также пороговые понижающиеся неполные факторизации Холесского (ICT(tau)) симметричных, положительных определенных разреженных матриц. Эти факторизации являются аналогами неполных LU-факторизаций выше и имеют многие из тех же характеристик. Например:

A = delsq(numgrid('S',200));
nnz(A)

ans =

      195228

nnz(chol(A,'lower'))

ans =

     7762589

L = ichol(A);
nnz(L)

ans =

      117216

norm(A-(L*L').*spones(A),'fro')./norm(A,'fro')

ans =

   3.5805e-17

opts.type = 'ict';
opts.droptol = 1e-4;
L = ichol(A,opts);
nnz(L)

ans =

     1166754

norm(A-L*L','fro')./norm(A,'fro')

ans =

   2.3997e-04

IC(0) имеет ненули только в шаблоне более низкого треугольника A, и на шаблоне A, продукте соответствий факторов. Кроме того, факторы ICT(1e-4) значительно более разреженны, чем полный Фактор Холесского, и относительная погрешность между A and L*L' находится на том же порядке допуска отбрасывания. Важно отметить, что различающийся факторы, обеспеченные chol, факторы по умолчанию, обеспеченные ichol, являются нижними треугольными. Смотрите страницу с описанием ichol для получения дополнительной информации.

Прямые методы

Прямые методы обычно быстрее и в более общем плане применимы, чем косвенные методы, если существует достаточно устройства хранения данных, доступного, чтобы выполнить их. Итерационные методы обычно применимы к ограниченным случаям уравнений и зависят от свойств как диагональное преобладание или существование базового дифференциального оператора. Прямые методы реализованы в ядре программного обеспечения MATLAB и сделаны максимально эффективными для общих классов матриц. Итерационные методы обычно реализуются в файлах языка MATLAB и могут использовать прямое решение подпроблем или предварительных формирователей.

Используя Различное Предварительное упорядоченное расположение.  Если A не является диагональным, не соединен, треугольным, или перестановка треугольной матрицы, наклонная черта влево (\) переупорядочивает индексы A, чтобы уменьшать сумму временной замены — то есть, количество ненулевых записей, которые добавляются к разреженным матрицам факторизации. Новое упорядоченное расположение, названное предварительным упорядоченным расположением, выполняется перед факторизацией A. В некоторых случаях вы можете смочь обеспечить лучшее предварительное упорядоченное расположение, чем то, используемое алгоритмом наклонной черты влево.

Чтобы использовать различное предварительное упорядоченное расположение, сначала выключите оба из автоматических предварительных упорядоченных расположений, которые наклонная черта влево может выполнить по умолчанию, с помощью функционального spparms можно следующим образом:

defaultParms = spparms;
spparms('autoamd',0);
spparms('autommd',0);

Теперь, принятие вас создало вектор перестановки p, который задает предварительное упорядоченное расположение индексов A, примените наклонную черту влево к матричному A(:,p), столбцы которого являются столбцами A, переставленного согласно векторному p.

x = A(:,p) \ b;
x(p) = x;
spparms(currentParms);

spparms(defaultParms) команды восстанавливает средства управления к их предшествующему состоянию, в случае, если вы используете A\b позже, не задавая соответствующее предварительное упорядоченное расположение.

Итеративные методы

Одиннадцать функций доступны, которые реализуют итерационные методы для разреженных систем одновременных линейных систем.

Эти методы разработаны, чтобы решить A x = b или минимизировать норму b – A x. Для Предобусловленного Метода сопряженных градиентов, pcg, A должен быть симметричной, положительной определенной матрицей. minres и symmlq могут использоваться на симметричных неопределенных матрицах. Для lsqr, матричная потребность не быть квадратным. Другие семь могут обработать несимметричные, квадратные матрицы, и каждый метод обладает отличным преимуществом.

Все одиннадцать методов могут использовать предварительные формирователи. Линейная система

заменяется эквивалентной системой

Предварительный формирователь M выбран, чтобы ускорить сходимость итерационного метода. Во многих случаях предварительные формирователи происходят естественно в математической модели. Дифференциальное уравнение с частными производными с переменными коэффициентами может быть аппроксимировано одним с постоянными коэффициентами, например. Неполные матричные факторизации могут использоваться в отсутствие естественных предварительных формирователей.

Приближение конечной разности с пятью точками к уравнению Лапласа на квадратной, двумерной области обеспечивает пример. Следующие операторы используют предобусловленный предварительный формирователь метода сопряженных градиентов M = L *L', где L является заполнением нулями неполного Фактора Холесского A.

A = delsq(numgrid('S',50));
b = ones(size(A,1),1);
tol = 1e-3;
maxit = 100;
L = ichol(A);
[x,flag,err,iter,res] = pcg(A,b,tol,maxit,L,L');

Двадцать одна итерация требуется, чтобы достигать предписанной точности. С другой стороны, использование различного предварительного формирователя может привести к лучшим результатам. Например, с помощью ichol, чтобы создать измененный неполный Холесский, предписанной точности встречают только после 15 итераций:

L = ichol(A,struct('type','nofill','michol','on'));
[x,flag,err,iter,res] = pcg(A,b,tol,maxit,L,L');

Справочная информация об этих итерационных методах и неполных факторизациях доступна в [2] и [6].