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 в квадрат не вызывает проблем. Однако возведение 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) и грант от Sandia National Lab. См. 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) и Эсмондом Нгом (Oak Ridge National Laboratory). (см. 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) — это индекс правого столбца, который не отсортирован или содержит дублирующиеся элементы, или 0, если такого столбца нет. stats(6) — это последний встреченный дублирующийся или неупорядоченный индекс строки в столбце, указанном stats(5), или 0, если такого индекса строки нет. stats(7) — это количество дублирующихся или неупорядоченных индексов строк. stats(8:20) всегда равен нулю в текущей версии CCOLAMD (зарезервировано для будущего использования).

Авторами кода являются С. Ларимор, Т. Дэвис и С. Раджаманикам в сотрудничестве с Дж. Бильбертом и Э. Нгом. Поддерживается Национальным научным фондом (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 (S)

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

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

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

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

Используемый метод описан в: А. Потен и К.-Дж. Фан. Вычисление блочно-треугольной формы разреженной матрицы. 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.