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 для примера структуры простой положительно определенной матрицы.

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

Авторы кода — 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 — это необязательный вектор от 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 (зарезервировано для будущего использования).

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

Авторы кода — Stefan I. Larimore и Timothy A. Davis. Алгоритм был разработан в сотрудничестве с John Gilbert, Xerox PARC, и Esmond Ng, 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) — это индекс правого столбца, который не отсортирован или содержит дубликаты, или ноль, если такого столбца нет. 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, 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. Вычисление блочно-треугольной формы разреженной матрицы. 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 (зарезервировано для будущего использования).

Порядок следуется посторядированием дерева элиминации столбцов.

Авторами кода являются Stefan I. Larimore и Timothy A. Davis. Алгоритм был разработан в сотрудничестве с Джоном Гилбертом, 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 по вычислительной математике, ISBN 0-13-165274-5, 1981.

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