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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Если S разреженная матрица, следующая команда возвращает три разреженных матрицы LU, и 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].

Самое маленькое собственное значение разреженной матрицы

В этом примере показано, как найти самое маленькое собственное значение и собственный вектор разреженной матрицы.

Настройте Дифференциальный оператор Лапласа с пятью точками на 65 65 сетка в L-образной, двумерной области.

L = numgrid('L',65);
A = delsq(L);

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

size(A)

ans = 1×2

        2945        2945

nnz(A)

ans = 14473

A матрица порядка 2945 с 14 473 ненулевыми элементами.

Вычислите самое маленькое собственное значение и собственный вектор.

[v,d] = eigs(A,1,'smallestabs');

Распределите компоненты собственного вектора по соответствующим узлам решетки и произведите контурный график результата.

L(L>0) = full(v(L(L>0)));
x = -1:1/32:1;
contour(x,x,L,15)
axis square