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

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

Вычислительная сложность разреженных операций пропорциональна 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, разреженная единичная матрица, для вставки -3 с на диагональ B.

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

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

Вот ненулевые счетчики для трёх сочетаний.

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

The 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».

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

R = qr(S)

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

R'*R = S'*S

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

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

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

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

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

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

[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)). The ILU(0) никогда не поворачивается, и результирующие множители имеют только ненули в положениях, где входная матрица имела ненули. Оба ILUC(tau) и ILUTP(tau)однако выполняйте отбрасывание на основе порога с пользовательским допуском отбрасывания tau.

Для примера:

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

ans =

        7840

nnz(lu(A))

ans =

      126478

показывает, что A имеет 7840 ненулевых, и его полная LU-факторизация имеет 126478 ненулевых. С другой стороны, следующий код показывает различные выходы 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

Эти вычисления показывают, что коэффициенты нулевого заполнения имеют 7840 ненулевых, ILUTP(1e-4) факторы имеют 31147 ненули, и ILUC(1e-4) факторы имеют 31083 ненулевых. Кроме того, относительная погрешность продукта коэффициентов нулевого заполнения по существу равна нулю на шаблоне 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 и L*L' находится в том же порядке, что и допуск. Важно отметить, что в отличие от факторов, обеспечиваемых chol, коэффициенты по умолчанию, предоставляемые ichol являются нижними треугольными. См. ichol Страница с описанием для получения дополнительной информации.

Наименьшее Собственное значение Разреженной Матрицы

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

Установите пятиточечный Дифференциального оператора Лапласа на сетке 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