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)

Если значение ненулевое, выводится статистика и значения knobs.

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) и грантом Национальной лаборатории Сандия. Описание ccolamd, csymamd, amd, colamd, symamd и других родственных алгоритмов упорядочения см. на странице http://faculty.cse.tamu.edu/davis/suitesparse.html.

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

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

Авторы кода — Stefan I. Larimore и Timothy A. Davis. Алгоритм разработан совместно с John Gilbert из Xerox PARC и Esmond Ng из Окриджской национальной лаборатории (см. 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) и грантом Национальной лаборатории Сандия. Описание ccolamd, colamd, csymamd, amd, colamd, symamd и других родственных алгоритмов упорядочения см. на странице http://faculty.cse.tamu.edu/davis/suitesparse.html.

См. также: 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. Алгоритм разработан совместно с John Gilbert из Xerox PARC и Esmond Ng из Окриджской национальной лаборатории (см. 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.