Операции разреженной матрицы

Эффективность операций

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

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

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

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

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

  • Функции, которые принимают матрицу и возвращают скаляр или вектор постоянного размера всегда, производят выход в полном формате устройства хранения данных. Например, size функционируйте всегда возвращает полный вектор, полон ли его вход или разрежен.

  • Функции, которые принимают скаляры или векторы и возвращают матрицы, такие как zerosединицыrand, и eye, всегда возвращайте полные результаты. Это необходимо, чтобы не вводить разреженность неожиданно. Разреженный аналог zeros(m,n) просто sparse(m,n). Разреженные аналоги rand и eye sprand и speye, соответственно. Нет никакого разреженного аналога для функционального ones.

  • Унарные функции, которые принимают матрицу и возвращают матрицу или вектор, сохраняют класс памяти операнда. Если S разреженная матрица, затем chol(S) также разреженная матрица и diag(S) разреженный вектор. Постолбцовые функции, такие как max и sum также возвратите разреженные векторы, даже при том, что эти векторы могут быть совершенно ненулевыми. Важными исключениями к этому правилу является sparse и full функции.

  • Бинарные операторы приводят к разреженным результатам, если оба операнда являются разреженными, и полными результатами, если оба полны. Для смешанных операндов результат полон, если операция не сохраняет разреженность. Если S разреженно и F полно, затем S+F, S*F, и F\S полны, в то время как S.*F и S&F разреженны. В некоторых случаях результат может быть разреженным даже при том, что матрица имеет немного нулевых элементов.

  • Конкатенация матриц с помощью любого cat функциональные или квадратные скобки приводят к разреженным результатам для смешанных операндов.

Перестановки и переупорядочение

Сочетание строк и столбцов разреженной матрицы S может быть представлен двумя способами:

  • Матрица перестановок P действия на строках S как P*S или на столбцах как S*P'.

  • Вектор сочетания p, который является полным вектором, содержащим сочетание 1:n, действия на строках S как S(p,:), или на столбцах как S(:,p).

Например:

p = [1 3 4 2 5]
I = eye(5,5);
P = I(p,:)
e = ones(4,1);
S = diag(11:11:55) + diag(e,1) + diag(e,-1)
p =

     1     3     4     2     5


P =

     1     0     0     0     0
     0     0     1     0     0
     0     0     0     1     0
     0     1     0     0     0
     0     0     0     0     1


S =

    11     1     0     0     0
     1    22     1     0     0
     0     1    33     1     0
     0     0     1    44     1
     0     0     0     1    55

Можно теперь попробовать некоторые сочетания с помощью вектора сочетания p и матрица перестановок P. Например, операторы S(p,:) и P*S возвратите ту же матрицу.

S(p,:)
ans =

    11     1     0     0     0
     0     1    33     1     0
     0     0     1    44     1
     1    22     1     0     0
     0     0     0     1    55
P*S
ans =

    11     1     0     0     0
     0     1    33     1     0
     0     0     1    44     1
     1    22     1     0     0
     0     0     0     1    55

Точно так же S(:,p) и S*P' оба производят

S*P'
ans =

    11     0     0     1     0
     1     1     0    22     0
     0    33     1     1     0
     0     1    44     0     1
     0     0     1     0    55

Если P разреженная матрица, затем оба представления используют устройство хранения данных, пропорциональное n и можно применить любого к S вовремя пропорциональный nnz(S). Векторное представление немного более компактно и эффективно, таким образом, различные стандартные программы сочетания разреженной матрицы, все возвращают полные векторы-строки за исключением вертящегося сочетания в LU (треугольная) факторизация, которая возвращает матрицу, совместимую с полной LU-факторизацией.

Преобразовывать между двумя представлениями сочетания:

n = 5;
I = speye(n);
Pr = I(p,:);
Pc = I(:,p);
pc = (1:n)*Pc
pc =

     1     3     4     2     5
pr = (Pr*(1:n)')'
pr =

     1     3     4     2     5

Инверсия P просто R = P'. Можно вычислить инверсию p с   r(p) = 1:n.

r(p) = 1:5
r =

     1     4     2     3     5

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

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

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

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

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

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

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

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

  • symamd — Используйте с симметричными матрицами.

  • colamd — Используйте с несимметричными матрицами и симметричными матрицами формы A*A' или A'*A.

Смотрите Переупорядочение и Факторизацию Разреженных матриц для примера с помощью 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.

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

  • lu(B) (Исходный): 1022

  • lu(B(r,r)) (Обратный алгоритм Катхилла-Макки): 968

  • lu(B(m,m)) (Аппроксимируйте минимальную степень): 636

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

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 позволяет решение разреженных проблем наименьших квадратов

minimizeAx-b2

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

[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
показывает тот A имеет 195 228 ненулей, и его полная факторизация Холесского без переупорядочения имеет 7 762 589 ненулей. В отличие от этого:
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 страница с описанием для получения дополнительной информации.

Собственные значения и сингулярные значения

Две функции доступны, которые вычисляют несколько заданных собственных значений или сингулярных значений. svds основан на eigs.

Функции, чтобы вычислить несколько собственных значений или сингулярных значений

Функция

Описание

eigs

Немного собственных значений

svds

Немного сингулярных значений

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

Оператор

[V,lambda] = eigs(A,k,sigma)

находит k собственные значения и соответствующие собственные вектора матричного A это является самым близким “сдвиг” sigma. Если sigma не использован, собственные значения, самые большие в величине, найдены. Если sigma нуль, собственные значения, самые маленькие в величине, найдены. Вторая матрица, B, может быть включен для обобщенной задачи о собственных значениях: Aυ = λBυ.

Оператор

[U,S,V] = svds(A,k)

находит k самые большие сингулярные значения A и

[U,S,V] = svds(A,k,'smallest')

находит k самые маленькие сингулярные значения.

Числовые методы используются в eigs и svds описаны в [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

Ссылки

[1] Amestoy, P. R., Т. А. Дэвис и я. S. Подновите, “Аппроксимированный Минимальный Алгоритм Упорядоченного расположения Степени”, SIAM Journal на Анализе матрицы и Приложения, Издание 17, № 4, октябрь 1996, стр 886-905.

[2] Барретт, R., М. Берри, Т. Ф. Чан, и др., Шаблоны для Решения Линейных систем: Базовые блоки для Итерационных методов, SIAM, Филадельфия, 1994.

[3] Дэвис, T.A., Гильберт, J. R. Larimore, S.I., Ын, E., Пейтон, B., “Столбец Аппроксимированный Минимальный Алгоритм Упорядоченного расположения Степени”, Proc. Конференция SIAM по Прикладной Линейной алгебре, октябрь 1997, p. 29.

[4] Гильберт, Джон Р., Клив Молер и Роберт Шрайбер, “Разреженные матрицы в MATLAB: Разработка и реализация”, SIAM J. Matrix Anal. Appl., Издание 13, № 1, январь 1992, стр 333-356.

[5] Larimore, S. I. аппроксимированный минимальный столбец степени, заказывая алгоритм, тезис MS, отдел информатики и вычислительной техники и разработки, Флоридского университета, Гейнсвилл, FL, 1998.

[6] Саад, Yousef, итерационные методы для разреженных линейных уравнений. PWS Publishing Company, 1996.

[7] Karypis, Джордж и Випин Кумар. "Быстрая и Высококачественная Многоуровневая Схема Разделения Неправильных Графиков". SIAM Journal на Научных вычислениях. Издание 20, Номер 1, 1999, стр 359–392.

Похожие темы