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

Рисунок 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).

Рисунок 22.4: Структура непереставленного разложения Холецкого вышеуказанной матрицы.

Рисунок 22.5: Структура переставленного разложения Холецкого вышеуказанной матрицы.

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

Авторами кода являются S. Larimore, T. Davis и S. Rajamanickam в сотрудничестве с J. Bilbert и E. Ng. Поддержка Национального научного фонда (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 — это необязательный вектор входных данных от одного до трех элементов. Если 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)

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

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 (зарезервировано для использования в будущем).

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

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

: p = dmperm (S)

: [p, q, r, S] = dmperm (S)

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

При использовании одного выходного аргумента dmperm выполняет перестановку строк p таким образом, что S(p,:) не имеет нулевых элементов на диагонали.

При вызове с двумя или более выходными аргументами возвращает перестановки строк и столбцов, такие что S(p, q) имеет блочно-треугольную форму. Значения r и S определяют границы блоков. Если S квадратная, то r == S.

Используемый метод описан в: 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, и Эсмондом Нгом, Национальной лабораторией Ок-Ридж. (см. 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. Сведение ширины разреженных симметричных матриц к минимуму. Труды 24-й национальной конференции ACM, 157–172 1969, Brandon Press, Нью-Джерси.

A. George, J.W.H. Liu. Компьютерное решение больших разреженных систем с положительно определёнными матрицами, Prentice Hall Series in Computational Mathematics, ISBN 0-13-165274-5, 1981.

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