22.1.4.3 Математические аспекты ¶

Была предпринята попытка сделать поведение разреженных матриц идентичным поведению их полных аналогов. Однако существуют определенные различия, особенно по сравнению с реализациями разреженных матриц в других продуктах.

Во-первых, операторы "./" и ".^" следует использовать с осторожностью. Подумайте, что дадут следующие примеры:

s = speye (4);
  a1 = s .^ 2;
  a2 = s .^ s;
  a3 = s .^ -2;
  a4 = s ./ 2;
  a5 = 2 ./ s;
  a6 = s ./ s;

Первый пример возведения s в степень 2 не вызывает проблем. Однако возведение s в степень, заданную элементом s, подразумевает большое количество слагаемых 0 .^ 0, что равно 1. В этом случае s .^ s получается полная матрица.

Аналогично, s .^ -2 включает в себя слагаемые типа 0 .^ -2, что равно бесконечности, и поэтому s .^ -2 также является полной матрицей.

Для оператора "./" s ./ 2 проблем нет, но 2 ./ s также включает в себя большое количество слагаемых, равных бесконечности, и также является полной матрицей. В случае s ./ s встречаются слагаемые типа 0 ./ 0, что является NaN, и, следовательно, это также полная матрица, где нулевые элементы s заполнены значениями NaN.

Это поведение соответствует поведению полных матриц, но не соответствует поведению разреженных матриц в других продуктах.

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

a = 0 ./ [-1, 1; 1, -1];
 b = 1 ./ a
 ⇒ -Inf            Inf
     Inf           -Inf
 c = 1 ./ sparse (a)
 ⇒  Inf            Inf
     Inf            Inf

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

В общем случае любая функция или оператор, используемый с разреженной матрицей, приведет к разреженной матрице с таким же или большим количеством ненулевых элементов, чем у исходной матрицы. Это особенно справедливо для важного случая факторизации разреженных матриц. Обычно это решается путем переупорядочения матрицы таким образом, чтобы ее факторизация была более разреженной, чем факторизация исходной матрицы. То есть факторизация L * U = P * S * Q имеет более разреженные члены L и U по сравнению с эквивалентной факторизацией L * U = S.

Доступно несколько функций для переупорядочения в зависимости от типа факторизуемой матрицы. Если матрица симметрично положительно определена, следует использовать symamd или csymamd. В противном случае следует использовать amd, colamd или ccolamd. Для полноты также доступны функции переупорядочения colperm и randperm.

См. Рисунок 22.3 для примера структуры простой положительно определенной матрицы.

Стандартную факторизацию Холецкого этой матрицы можно получить с помощью той же команды, что и для полной матрицы. Это можно визуализировать с помощью команды r = chol (A); spy (r);. См. Рисунок 22.4. Исходная матрица имела 598 ненулевых элементов, а эта факторизация Холецкого — 10200, при этом хранится только половина симметричной матрицы. Это значительный уровень заполнения, и хотя это не проблема для такого простого тестового случая, это может представлять значительную нагрузку при работе с другими разреженными матрицами.

Соответствующее сохраняющее разреженность переупорядочение исходной матрицы задается функцией symamd, а визуализацию факторизации с этим переупорядочением можно получить с помощью команды q = symamd (A); r = chol (A(q,q)); spy (r). Это даёт 399 ненулевых элементов, что является значительным улучшением.

Сама факторизация Холецкого может использоваться для определения соответствующего сохраняющего разреженность переупорядочения матрицы во время факторизации. В этом случае это можно получить с тремя возвращаемыми аргументами как [r, p, q] = chol (A); spy (r).

В случае несимметричной матрицы соответствующее сохраняющее разреженность переупорядочение — colamd, а визуализацию факторизации с этим переупорядочением можно получить с помощью команды q = colamd (A); [l, u, p] = lu (A(:,q)); spy (l+u).

Наконец, Octave неявно переупорядочивает матрицу при использовании операторов div (/) и ldiv (\), поэтому пользователю не нужно явно переупорядочивать матрицу для максимальной производительности.

: p = amd (S) ¶

: p = amd (S, opts) ¶

Возвращает приближенное переупорядочение по минимальной степени для матрицы.

Это переупорядочение таково, что факторизация Холецкого S (p, p) имеет тенденцию быть более разреженной, чем факторизация Холецкого самой S. amd обычно быстрее, чем symamd, но служит аналогичной цели.

Необязательный параметр opts — это структура, которая управляет поведением amd. Поля структуры:

opts.dense

Определяет, что amd считает плотной строкой или столбцом входной матрицы. Строки или столбцы с более чем max (16, (dense * sqrt (n))) элементами, где n — порядок матрицы S, игнорируются amd при вычислении перестановки. Значение dense должно быть положительным скаляром, а значение по умолчанию равно 10,0.

opts.aggressive

Если это значение — ненулевой скаляр, amd выполняет агрессивное поглощение. По умолчанию агрессивное поглощение не выполняется.

Автор кода — Timothy A. Davis (см. http://faculty.cse.tamu.edu/davis/suitesparse.html).

См. также: symamd, colamd.

: p = ccolamd (S) ¶

: p = ccolamd (S, knobs) ¶

: p = ccolamd (S, knobs, cmember) ¶

: [p, stats] = ccolamd (…) ¶

Перестановка приближённой минимальной степени столбцов с ограничениями.

p = ccolamd (S) возвращает вектор перестановки приближённой минимальной степени столбцов для разреженной матрицы S. Для несимметричной матрицы S, S(:, p) имеет тенденцию к более разреженным факторам LU, чем S. chol (S(:, p)' * S(:, p)) также имеет тенденцию к большей разреженности, чем chol (S' * S). p = ccolamd (S, 1) оптимизирует порядок для lu (S(:, p)). Порядок следует за пост-упорядочиванием дерева исключения столбцов.

knobs — это необязательный вектор входных данных от 1 до 5 элементов, со значением по умолчанию [0 10 10 1 0] если отсутствует или пуст. Отсутствующие элементы устанавливаются по умолчанию.

knobs(1)

если ненулевое, порядок оптимизирован для lu (S(:, p)). Он будет плохим порядком для chol (S(:, p)' * S(:, p)). Это самый важный параметр для ccolamd.

knobs(2)

если S имеет размер m×n, строки с более чем max (16, knobs(2) * sqrt (n)) элементами игнорируются.

knobs(3)

столбцы с более чем max (16, knobs(3) * sqrt (min (m, n))) элементами игнорируются и упорядочиваются последними в выходной перестановке (подчиняясь ограничениям cmember).

knobs(4)

если ненулевое, выполняется агрессивное поглощение.

knobs(5)

если ненулевое, печатаются статистика и параметры.

cmember — это необязательный вектор длиной n. Он определяет ограничения на порядок столбцов. Если cmember(j) = c, то столбец j находится в наборе ограничений c (c должно находиться в диапазоне от 1 до n). В выходной перестановке p все столбцы из набора 1 появляются первыми, за ними все столбцы из набора 2 и так далее. cmember = ones (1,n) если не указан или пуст. ccolamd (S, [], 1 : n) возвращает 1 : n

p = ccolamd (S) примерно то же, что и p = colamd (S). knobs и его значения по умолчанию отличаются. colamd всегда выполняет агрессивное поглощение и находит порядок, подходящий как для lu (S(:, p)), так и для chol (S(:, p)' * S(:, p)); он не может оптимизировать порядок для lu (S(:, p)) в той степени, в которой это может ccolamd (S, 1).

stats — это необязательный выходной вектор из 20 элементов, предоставляющий данные об упорядочивании и корректности входной матрицы S. Статистика упорядочения находится в stats(1 : 3). stats(1) и stats(2) — это количество плотных или пустых строк и столбцов, пропущенных CCOLAMD, а stats(3) — это количество сборок мусора, выполненных в внутренней структуре данных, используемой CCOLAMD (примерно размером 2.2 * nnz (S) + 4 * m + 7 * n целых чисел).

stats(4 : 7) предоставляют информацию о том, смог ли CCOLAMD продолжить работу. Матрица корректна, если stats(4) равно нулю, или 1, если некорректна. stats(5) — это индекс правого столбца, который не отсортирован или содержит повторяющиеся элементы, или ноль, если такого столбца нет. stats(6) — это последний обнаруженный повторяющийся или неотсортированный индекс строки в индексе столбца, заданном stats(5), или ноль, если такой индекс строки не существует. stats(7) — это количество повторяющихся или неотсортированных индексов строк. stats(8 : 20) всегда равен нулю в текущей версии CCOLAMD (зарезервировано для будущих применений).

Авторами кода являются С. Ларимор, Т. Дэвис и С. Раджаманичкам в сотрудничестве с Дж. Бильбертом и Э. Нгом. Поддержка Национального научного фонда (DMS-9504974, DMS-9803599, CCR-0203270) и грант от Национальной лаборатории Сандии. См. http://faculty.cse.tamu.edu/davis/suitesparse.html для ccolamd, csymamd, amd, colamd, symamd и других родственных упорядочений.

См. также: colamd, csymamd.

: p = colamd (S) ¶

: p = colamd (S, knobs) ¶

: [p, stats] = colamd (S) ¶

: [p, stats] = colamd (S, knobs) ¶

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

p = colamd (S) возвращает вектор перестановки столбцов приближённой минимальной степени для разреженной матрицы S. Для несимметричной матрицы S, S(:,p) имеет тенденцию к более разреженным факторам LU, чем S. Факторизация Холецкого S(:,p)' * S(:,p) также имеет тенденцию к большей разреженности, чем у S' * S.

knobs — это необязательный вектор от 1 до 3 элементов. Если S имеет размер m×n, то строки с более чем max(16,knobs(1)*sqrt(n)) элементами игнорируются. Столбцы с более чем max (16,knobs(2)*sqrt(min(m,n))) элементами удаляются перед упорядочиванием и упорядочиваются последними в выходной перестановке p. Если knobs(1) и knobs(2) < 0 соответственно, удаляются только полностью плотные строки или столбцы. Если knobs(3) ненулевое, stats и knobs выводятся на экран. По умолчанию — knobs = [10 10 0]. Обратите внимание, что knobs отличается от предыдущих версий colamd.

stats — это необязательный выходной вектор из 20 элементов, предоставляющий данные об упорядочивании и корректности входной матрицы S. Статистика упорядочения находится в stats(1:3). stats(1) и stats(2) — это количество плотных или пустых строк и столбцов, пропущенных COLAMD, а stats(3) — это количество сборок мусора, выполненных в внутренней структуре данных, используемой COLAMD (примерно размером 2.2 * nnz(S) + 4 * m + 7 * n целых чисел).

Встроенные функции Octave предназначены для генерации корректных разреженных матриц без дублирующихся элементов, с возрастающими индексами строк ненулевых элементов в каждом столбце, с неотрицательным числом элементов в каждом столбце (!) и так далее. Если матрица некорректна, COLAMD может или не может продолжить работу. Если есть повторяющиеся элементы (индекс строки появляется два или более раз в одном столбце) или если индексы строк в столбце не отсортированы, то COLAMD может исправить эти ошибки, игнорируя повторяющиеся элементы и сортируя каждый столбец своей внутренней копии матрицы S (входная матрица S не исправляется). Если матрица некорректна другими способами, COLAMD не может продолжить работу, выводится сообщение об ошибке, и выходные аргументы (p или stats) не возвращаются. COLAMD — это простой способ проверить разреженную матрицу на корректность.

stats(4:7) предоставляют информацию о том, смог ли COLAMD продолжить работу. Матрица корректна, если stats(4) равно нулю, или 1, если некорректна. stats(5) — это индекс правого столбца, который не отсортирован или содержит повторяющиеся элементы, или ноль, если такого столбца нет. stats(6) — это последний обнаруженный повторяющийся или неотсортированный индекс строки в индексе столбца, заданном stats(5), или ноль, если такой индекс строки не существует. stats(7) — это количество повторяющихся или неотсортированных индексов строк. stats(8:20) всегда равен нулю в текущей версии COLAMD (зарезервировано для будущих применений).

После упорядочения следует пост-упорядочение дерева исключения столбцов.

Авторами кода являются Стефан И. Ларимор и Тимоти А. Дэвис. Алгоритм был разработан в сотрудничестве с Джоном Гильбертом, Xerox PARC, и Эсмондом Нгом, Национальной лабораторией Ок-Риджа. (см. http://faculty.cse.tamu.edu/davis/suitesparse.html)

См. также: colperm, symamd, ccolamd.

: p = colperm (s) ¶

Возвращает перестановки столбцов, такие что столбцы s(:, p) упорядочены по возрастающему числу ненулевых элементов.

Если s симметрична, то p выбирается таким образом, что s(p, p) упорядочивает строки и столбцы по возрастающему числу ненулевых элементов.

: p = csymamd (S) ¶

: p = csymamd (S, knobs) ¶

: p = csymamd (S, knobs, cmember) ¶

: [p, stats] = csymamd (…) ¶

Для симметричной положительно определённой матрицы S, вернуть вектор перестановки p такой, что S(p,p) имеет более разреженный фактор Холецкого, чем S.

Иногда csymamd хорошо работает и для симметричных неопределённых матриц. Матрица S предполагается симметричной; используется только строго нижняя треугольная часть. S должна быть квадратной. Порядок определяется после упорядочивания по дереву исключения.

knobs — это необязательный вектор входных данных от 1 до 3 элементов, со значением по умолчанию [10 1 0]. Отсутствующие элементы устанавливаются по умолчанию.

knobs(1)

Если S имеет размер n×n, то строки и столбцы с более чем max(16,knobs(1)*sqrt(n)) элементами игнорируются и упорядочиваются в выходной перестановке в последнюю очередь (подчиняясь ограничениям cmember).

knobs(2)

Если ненулевое, выполняется агрессивное поглощение.

knobs(3)

Если ненулевое, выводятся статистика и параметры.

cmember — это необязательный вектор длины n. Он определяет ограничения на упорядочивание. Если cmember(j) = S, то строка/столбец j входит в набор ограничений c (c должен быть в диапазоне от 1 до n). В выходной перестановке p, строки/столбцы из набора 1 появляются первыми, за ними все строки/столбцы из набора 2 и так далее. cmember = ones (1,n) если не задан или пуст. csymamd (S,[],1:n) возвращает 1:n.

p = csymamd (S) аналогично p = symamd (S). knobs и его значения по умолчанию отличаются.

stats(4:7) сообщают, удалось ли CCOLAMD продолжить. Матрица корректна, если stats(4) равно нулю, или 1 если некорректна. stats(5) — это индекс правого столбца, который неупорядочен или содержит повторяющиеся элементы, или ноль, если такого столбца нет. stats(6) — это последний встреченный дублирующийся или не в порядке индекс строки в индексе столбца, заданном stats(5), или ноль, если такого индекса строки нет. stats(7) — это количество дублирующихся или не в порядке индексов строк. stats(8:20) всегда равен нулю в текущей версии CCOLAMD (зарезервировано для будущего использования).

Авторы кода — S. Larimore, T. Davis и S. Rajamanickam в сотрудничестве с J. Bilbert и E. Ng. Поддерживается Национальным научным фондом (DMS-9504974, DMS-9803599, CCR-0203270) и грантом от Sandia National Lab. См. http://faculty.cse.tamu.edu/davis/suitesparse.html для ccolamd, colamd, csymamd, amd, colamd, symamd и других связанных упорядочений.

См. также: symamd, ccolamd.

: p = dmperm (A) ¶

: [p, q, r, s, cc, rr] = dmperm (A) ¶

Выполнить перестановку по методу Дюльмежа-Мендельсона для разреженной матрицы A.

С единственным выходным аргументом dmperm, вернуть максимальное соответствие p такое, что p(j) = i если столбец j сопоставлен строке i, или 0, если столбец j не сопоставлен. Если A квадратная и имеет полный структурный ранг, p — это перестановка строк, и A(p,:) имеет ненулевую диагональ. Структурный ранг A равен sprank(A) = sum(p>0).

При вызове с двумя или более выходными аргументами, вернуть разложение Дюльмежа-Мендельсона матрицы A. p и q — векторы перестановок. cc и rr — векторы длины 5. c = A(p,q) разбивается на набор из 4×4 блоков:

A11 A12 A13 A14
    0  0   A23 A24
    0  0    0  A34
    0  0    0  A44

где A12, A23, и A34 — квадратные с ненулевой диагональю. Столбцы A11 — это несопоставленные столбцы, а строки A44 — это несопоставленные строки. Любой из этих блоков может быть пустым. В «грубом» разложении (i,j)-й блок равен C(rr(i):rr(i+1)-1,cc(j):cc(j+1)-1). С точки зрения линейной системы, [A11 A12] — это недоопределённая часть системы (всегда прямоугольная с большим количеством столбцов и строк, или 0×0), A23 — это хорошо определённая часть системы (всегда квадратная), а [A34 ; A44] — это переопределённая часть системы (всегда прямоугольная с большим количеством строк, чем столбцов, или 0×0).

Структурный ранг A равен sprank (A) = rr(4)-1, что является верхней границей числового ранга A. sprank(A) = rank(full(sprand(A))) с вероятностью 1 в точном арифметике.

A23 подматрица дополнительно разбивается на блочно-верхнетреугольную форму через «тонкое» разложение (сильно связные компоненты A23). Если A квадратная и структурно невырожденная, то A23 — это вся матрица.

C(r(i):r(i+1)-1,s(j):s(j+1)-1) — это (i,j)-й блок тонкого разложения. (1,1) блок — это прямоугольный блок [A11 A12], если этот блок не 0×0. Блок (b,b) — это прямоугольный блок [A34 ; A44], если этот блок не 0×0, где b = length(r)-1. Все остальные блоки вида C(r(i):r(i+1)-1,s(i):s(i+1)-1) — это диагональные блоки A23, и являются квадратными с ненулевой диагональю.

Используемый метод описан в: A. Pothen & C.-J. Fan. Computing the Block Triangular Form of a Sparse Matrix. ACM Trans. Math. Software, 16(4):303–324, 1990.

См. также: colamd, ccolamd.

: p = symamd (S) ¶

: p = symamd (S, knobs) ¶

: [p, stats] = symamd (S) ¶

: [p, stats] = symamd (S, knobs) ¶

Для симметричной положительно определённой матрицы S возвращает вектор перестановок p такой, что S(p, p) имеет более разреженный фактор Холецкого, чем S.

Иногда symamd хорошо работает и для симметричных неопределённых матриц. Предполагается, что матрица S симметрична; используется только строго нижняя треугольная часть. S должна быть квадратной.

knobs — это необязательный входной вектор из одного или двух элементов. Если S является n×n матрицей, то строки и столбцы с более чем max (16,knobs(1)*sqrt(n)) элементами удаляются перед упорядочиванием и упорядочиваются в конце выходной перестановки p. Если knobs(1) < 0, то строки/столбцы не удаляются. Если knobs(2) не равно нулю, выводятся stats и knobs. По умолчанию knobs = [10 0]. Обратите внимание, что knobs отличается от предыдущих версий symamd.

stats — это необязательный выходной вектор из 20 элементов, который предоставляет данные об упорядочивании и корректности входной матрицы S. Статистики упорядочивания содержатся в stats(1:3). stats(1) = stats(2) — это количество плотных или пустых строк и столбцов, проигнорированных SYMAMD, и stats(3) — это количество выполненных сборки мусора на внутренней структуре данных, используемой SYMAMD (приблизительно размером 8.4 * nnz (tril (S, -1)) + 9 * n целых чисел).

Встроенные функции Octave предназначены для генерации корректных разреженных матриц без дублирования элементов, с возрастающими индексами строк ненулевых элементов в каждом столбце, с неотрицательным числом элементов в каждом столбце (!) и т. д. Если матрица некорректна, SYMAMD может или не может продолжить работу. Если есть дублируемые элементы (индекс строки встречается два или более раз в одном столбце) или если индексы строк в столбце неупорядочены, SYMAMD может исправить эти ошибки, проигнорировав дублируемые элементы и отсортировав каждый столбец своей внутренней копии матрицы S (входная матрица S не исправляется). Если матрица некорректна другими способами, SYMAMD не может продолжить, выводится сообщение об ошибке, и никакие выходные аргументы (p или stats) не возвращаются. Таким образом, SYMAMD — простой способ проверить разреженную матрицу на корректность.

stats(4:7) предоставляют информацию о том, смог ли SYMAMD продолжить работу. Матрица корректна, если stats (4) равно нулю, или 1, если некорректна. stats(5) — это индекс правого столбца, который не отсортирован или содержит дублируемые элементы, или ноль, если такого столбца нет. stats(6) — это последний увиденный дублируемый или неупорядоченный индекс строки в столбце, указанном stats(5), или ноль, если такого индекса строки нет. stats(7) — это количество дублируемых или неупорядоченных индексов строк. stats(8:20) всегда равен нулю в текущей версии SYMAMD (зарезервировано для будущего использования).

За упорядочиванием следует пост-упорядочивание дерева устранения столбцов.

Авторами кода являются Стефан И. Ларимор и Тимоти А. Дэвис. Алгоритм был разработан совместно с Джоном Гильбертом (Xerox PARC) и Эсмондом Нгом (Oak Ridge National Laboratory). (см. http://faculty.cse.tamu.edu/davis/suitesparse.html)

См. также: colperm, colamd.

: p = symrcm (S) ¶

Возвращает симметричную перестановку обращённой схемы Катилла-Макки для S.

p — вектор перестановок, такой что S(p, p) имеет элементы диагонали, более близкие к диагонали, чем S. Это хорошее предварительное упорядочение для LU или факторизации Холецкого для матриц, которые происходят из «длинных и узких» задач. Работает как для симметричных, так и для несимметричных S.

Алгоритм представляет собой эвристический подход к задаче минимизации ширины, которая является NP-полной. Реализация основана на описаниях из

E. Cuthill, J. McKee. Reducing the Bandwidth of Sparse Symmetric Matrices. Proceedings of the 24th ACM National Conference, 157–172 1969, Brandon Press, New Jersey.

A. George, J.W.H. Liu. Computer Solution of Large Sparse Positive Definite Systems, Prentice Hall Series in Computational Mathematics, ISBN 0-13-165274-5, 1981.

См. также: colperm, colamd, symamd.