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

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

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

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

Алгоритмические сведения

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

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

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

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

Переупорядочивание для уменьшения полосы пропускания

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

Приблизительный порядок минимальных степеней

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

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

Пример использования см. в разделе Переупорядочивание и факторизация разреженных матриц symamd.

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

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

Порядок вложенных рассечений

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

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

Факторизация логической единицы

Если 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 не использует какую-либо символическую предварительную факторизацию логической единицы для определения требований к памяти и предварительной настройки структур данных.

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

В этом примере показано влияние переупорядочивания и факторизации на разреженные матрицы.

Если вы получите хорошую перестановку столбцов 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, создает фракталоподобную структуру с большими блоками нулей.

Чтобы увидеть заполнение, сгенерированное в факторизации логической единицы шарика Баки, используйте 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», доступна.

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

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 функция производит три неполные нижние-верхние (ILU) факторизации: ILU с нулевым заполнением (ILU(0)), версия ILU Crout (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 имеет 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 функция обеспечивает нулевое заполнение неполных факторизаций Холеского (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 = 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