22.2 Линейная алгебра на разреженных матрицах

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

Дерево выбора способа решения линейного уравнения:

1. Если матрица диагональная, решите напрямую и перейдите к 8
2. Если матрица — это переставленная диагональная, решите напрямую, учитывая перестановки. Перейдите к 8
3. Если матрица квадратная, полосового типа, и если плотность полосы меньше, чем задано в spparms ("bandden"), продолжайте, в противном случае перейдите к 4.
   1. Если матрица трёхдиагональная, а правая часть не разреженная, продолжайте, в противном случае перейдите к 3b.
      1. Если матрица эрмитова, с положительной действительной диагональю, попробуйте факторизацию Холецкого с использованием LAPACK xPTSV.
      2. Если это не удалось или матрица не эрмитова с положительной действительной диагональю, используйте метод Гаусса с выбором главного элемента, используя LAPACK xGTSV, и перейдите к 8.
   2. Если матрица эрмитова с положительной действительной диагональю, попробуйте факторизацию Холецкого, используя LAPACK xPBTRF.
   3. Если это не удалось или матрица не эрмитова с положительной действительной диагональю, используйте метод Гаусса с выбором главного элемента, используя LAPACK xGBTRF, и перейдите к 8.
4. Если матрица верхнетреугольная или нижнетреугольная, выполните разреженную подстановку вперёд или назад, и перейдите к 8
5. Если матрица — верхнетреугольная матрица с перестановкой столбцов или нижнетреугольная матрица с перестановкой строк, выполните разреженную подстановку вперёд или назад, и перейдите к 8
6. Если матрица квадратная, эрмитова с положительной действительной диагональю, попробуйте разреженную факторизацию Холецкого с использованием CHOLMOD.
7. Если разреженная факторизация Холецкого не удалась или матрица не эрмитова с положительной действительной диагональю, и матрица квадратная, выполните факторизацию, решение и одну итерацию уточнения с использованием UMFPACK.
8. Если матрица не квадратная или любой из предыдущих решателей отмечает сингулярную или почти сингулярную матрицу, найдите решение с минимальной нормой, используя CXSPARSE¹⁰.

Плотность полосы определяется как количество ненулевых значений в полосе, делённое на общее количество значений в полной полосе. Решатели для полосовых матриц могут быть полностью отключены, используя spparms для установки bandden в 1 (т.е., spparms ("bandden", 1)).

Решатель QR факторизует проблему с помощью разложения Дюльмаге-Мендельсона, чтобы разделить проблему на блоки, которые могут рассматриваться как переопределённые, несколько хорошо определённых блоков и окончательный переопределённый блок. Для матриц с блоками сильно связанных узлов это большое преимущество, так как разложение LU может быть использовано для многих блоков. Это также значительно повышает вероятность нахождения решения для переопределённых задач, а не просто возвращения вектора NaN.

Все вышеперечисленные решатели могут рассчитать оценку числа обусловленности. Это можно использовать для обнаружения проблем с числовой устойчивостью в решении и принудительно использовать решение с минимальной нормой. Однако для узких полосовых, треугольных или диагональных матриц стоимость вычисления числа обусловленности значительна и может фактически превышать стоимость факторизации матрицы. Поэтому число обусловленности не вычисляется в этих случаях, и Octave использует более простые методы для обнаружения сингулярных матриц или базового кода LAPACK в случае полосовых матриц.

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

nest = normest (A)

nest = normest (A, tol)

[nest, iter] = normest (…)

Оценить 2-норму матрицы A с помощью анализа степенного ряда.

Это обычно используется для больших матриц, где стоимость вычисления norm (A) является неприемлемой, и приближение к 2-норме является приемлемым.

tol — это точность вычисления 2-нормы. По умолчанию tol равен 1e-6.

Необязательный выход iter возвращает число итераций, необходимых для normest сходимости.

См. также: normest1, norm, cond, condest.

nest = normest1 (A)

nest = normest1 (A, t)

nest = normest1 (A, t, x0)

nest = normest1 (Afun, t, x0, p1, p2, …)

[nest, v] = normest1 (A, …)

[nest, v, w] = normest1 (A, …)

[nest, v, w, iter] = normest1 (A, …)

Оценить 1-норму матрицы A с использованием блочного алгоритма.

normest1 лучше всего подходит для больших разреженных матриц, когда требуется только оценка нормы. Для матриц малого и среднего размера используйте norm (A, 1). Кроме того, normest1 может использоваться для оценки 1-нормы линейного оператора A, когда матрично-векторные произведения A * x и A' * x могут быть вычислены недорого. В этом случае вместо матрицы A используется функция Afun (flag, x); она должна возвращать:

• размерность n A, если flag равно "dim"
• true, если A — вещественный оператор, если flag равно "real"
• результат A * x, если flag равно "notransp"
• результат A' * x, если flag равно "transp"

Типичный случай — A, определённая b ^ m, в котором результат A * x может быть вычислен, не формируя явно b ^ m, путём:

y = x;
for i = 1:m
  y = b * y;
endfor

Параметры p1, p2, … являются аргументами Afun (flag, x, p1, p2, …).

Значение по умолчанию для t — 2. Алгоритм требует матрично-матричных произведений размером n x n и n x t.

Начальная матрица x0 должна иметь столбцы с единичной 1-нормой. По умолчанию начальная матрица x0 имеет первый столбец ones (n, 1) / n и, если t > 1, остальные столбцы с случайными элементами -1 / n, 1 / n, делёнными на n.

На выходе, nest — искомая оценка, v и w — векторы, такие что w = A * v, с norm (w, 1) = c * norm (v, 1). iter содержит в iter(1) число итераций (максимум жёстко задан 5) и в iter(2) общее число произведений A * x или A' * x, выполненных алгоритмом.

Примечание к алгоритму: normest1 использует случайные числа во время вычислений. Поэтому, если требуются согласованные результаты, "state" генератора случайных чисел следует зафиксировать перед вызовом normest1.

Ссылка: Н. Дж. Хайгам и Ф. Тиссер, Блочный алгоритм для оценки 1-нормы матрицы, с применением к 1-норменным псевдоспектрам, SIAM J. Matrix Anal. Appl., стр. 1185–1201, том 21, № 4, 2000.

См. также: normest, norm, cond, condest.

cest = condest (A)

cest = condest (A, t)

cest = condest (A, solvefun, t, p1, p2, …)

cest = condest (Afcn, solvefun, t, p1, p2, …)

[cest, v] = condest (…)

Оценить 1-норменное число обусловленности квадратной матрицы A с использованием t тестовых векторов и случайного 1-норменного оценщика.

Необязательный вход t задаёт количество тестовых векторов (по умолчанию 5).

Если матрица не явная, например, при оценке числа обусловленности A по разложению LU, condest использует следующие функции:

• - Afcn, которая должна возвращать
  • размерность n a, если flag равно "dim"
  • true, если a — вещественный оператор, если flag равно "real"
  • результат a * x, если flag равно "notransp"
  • результат a' * x, если flag равно "transp"
• - solvefun, которая должна возвращать
  • размерность n a, если flag равно "dim"
  • true, если a — вещественный оператор, если flag равно "real"
  • результат a \ x, если flag равно "notransp"
  • результат a' \ x, если flag равно "transp"

Параметры p1, p2, … являются аргументами Afcn (flag, x, p1, p2, …) и solvefcn (flag, x, p1, p2, …).

Основной выход — оценка числа обусловленности 1-нормы cest.

Необязательный второй выход — приближённый нулевой вектор, когда cest велик; он удовлетворяет уравнению norm (A*v, 1) == norm (A, 1) * norm (v, 1) / est.

Примечание к алгоритму: condest использует случайный алгоритм для приближения 1-норм. Поэтому, если требуются согласованные результаты, "state" генератора случайных чисел следует зафиксировать перед вызовом condest.

Ссылки:

• Н.Дж. Хайгам и Ф. Тиссер, Блочный алгоритм для оценки 1-нормы матрицы, с применением к 1-норменным псевдоспектрам. SIMAX том 21, № 4, стр. 1185–1201. https://dx.doi.org/10.1137/S0895479899356080
• Н.Дж. Хайгам и Ф. Тиссер, Блочный алгоритм для оценки 1-нормы матрицы, с применением к 1-норменным псевдоспектрам. https://citeseer.ist.psu.edu/223007.html

См. также: cond, norm, normest1, normest.

spparms ()

vals = spparms ()

[keys, vals] = spparms ()

val = spparms (key)

spparms (vals)

spparms ("default")

spparms ("tight")

spparms (key, val)

Запрос или установка параметров, используемых разреженными решателями и функциями факторизации.

Четыре первых вызова выше получают информацию о текущих настройках, в то время как остальные изменяют текущие настройки. Параметры хранятся в виде пар «ключ-значение», где значения являются числами с плавающей точкой, а ключи — одной из следующих строк:

‘spumoni’

Уровень вывода отладочной информации решателей (по умолчанию 0)

‘ths_rel’

Включено для совместимости. Не используется. (по умолчанию 1)

‘ths_abs’

Включено для совместимости. Не используется. (по умолчанию 1)

‘exact_d’

Включено для совместимости. Не используется. (по умолчанию 0)

‘supernd’

Включено для совместимости. Не используется. (по умолчанию 3)

‘rreduce’

Включено для совместимости. Не используется. (по умолчанию 3)

‘wh_frac’

Включено для совместимости. Не используется. (по умолчанию 0.5)

‘autommd’

Флаг, указывающий, будут ли операторы LU/QR и ’\’ и ’/’ автоматически использовать функции mmd, сохраняющие разреженность (по умолчанию 1)

‘autoamd’

Флаг, указывающий, будут ли операторы LU и ’\’ и ’/’ автоматически использовать функции amd, сохраняющие разреженность (по умолчанию 1)

‘piv_tol’

Порог выбора опорного элемента решателей UMFPACK (по умолчанию 0.1)

‘sym_tol’

Порог выбора опорного элемента симметричных решателей UMFPACK (по умолчанию 0.001)

‘bandden’

Плотность ненулевых элементов в полосовой матрице перед обработкой полосовыми решателями LAPACK (по умолчанию 0.5)

‘umfpack’

Флаг, указывающий, используются ли решатели UMFPACK или mmd для операций LU, ’\’ и ’/’ (по умолчанию 1)

Значение отдельных ключей можно установить с помощью spparms (key, val). Значения по умолчанию можно восстановить с помощью специального ключевого слова "default". Специальное ключевое слово "tight" можно использовать для настройки решателей mmd на поиск более разреженного решения с потенциальной потерей времени выполнения.

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

p = sprank (S)

Вычисление структурного ранга разреженной матрицы S.

Обратите внимание, что в этом вычислении используется только структура матрицы на основе перестановки Далмага-Мендельсона для блочно-треугольной формы. Таким образом, числовой ранг матрицы S ограничен sprank (S) >= rank (S). Игнорируя ошибки с плавающей точкой sprank (S) == rank (S).

См. также: dmperm.

[count, h, parent, post, R] = symbfact (S)

[…] = symbfact (S, typ)

[…] = symbfact (S, typ, mode)

Выполнение символического анализа факторизации разреженной матрицы S.

Вводные переменные:

S

S — вещественная или комплексная разреженная матрица.

typ

Тип факторизации и может быть одним из:

"sym" (по умолчанию)

Факторизация S. Предполагает, что S симметрична и использует верхнюю треугольную часть матрицы.

"col"

Факторизация S' * S.

"row"

Факторизация S * S'.

"lo"

Факторизация S'. Предполагает, что S симметрична и использует нижнюю треугольную часть матрицы.

mode

Если mode не указано, возвращается факторизация Холецкого для R. Если mode равно "lower" или "L", возвращается сопряжённая транспонированная R', которая является нижнетреугольной матрицей. Версия сопряжённой транспонированной быстрее и использует меньше памяти, но всё ещё возвращает те же значения для всех других выходов: count, h, parent и post.

Выходные переменные:

count

Число строк факторизации Холецкого, определяемое typ. Вычислительная сложность выполнения истинной факторизации с использованием chol равна sum (count .^ 2).

h

Высота дерева устранения.

parent

Само дерево устранения.

post

Разреженная булева матрица, структура которой соответствует факторизации Холецкого, определяемой typ.

См. также: chol, etree, treelayout.

Для неквадратных матриц пользователь также может использовать функцию spaugment для поиска решения в методе наименьших квадратов для линейного уравнения.

s = spaugment (A, c)

Создание дополненной матрицы A.

Она задаётся следующим образом:

[c * eye(m, m), A;
            A', zeros(n, n)]

Это связано с решением методом наименьших квадратов A \ b, следующим образом:

s * [ r / c; x] = [ b, zeros(n, columns(b)) ]

где r — остаточная ошибка

r = b - A * x

Поскольку матрица s симметрична и не определена, её можно разложить с помощью lu, и поэтому можно найти решение с минимальной нормой без необходимости в факторизации qr. Поскольку остаточная ошибка будет zeros (m, m) для недоопределённых задач, примером является

m = 11; n = 10; mn = max (m, n);
A = spdiags ([ones(mn,1), 10*ones(mn,1), -ones(mn,1)],
             [-1, 0, 1], m, n);
x0 = A \ ones (m,1);
s = spaugment (A);
[L, U, P, Q] = lu (s);
x1 = Q * (U \ (L \ (P  * [ones(m,1); zeros(n,1)])));
x1 = x1(end - n + 1 : end);

Для нахождения решения переопределённой задачи требуется оценка остаточной ошибки r, поэтому более сложно сформулировать решение с минимальной нормой, используя функцию spaugment.

В общем случае, оператор левого деления более стабилен и быстрее, чем использование функции spaugment.

См. также: mldivide.

В заключение, функция eigs может использоваться для вычисления ограниченного числа собственных значений и собственных векторов на основе критерия отбора, и аналогично для svds, которая вычисляет ограниченное число сингулярных значений и векторов.

d = eigs (A)

d = eigs (A, k)

d = eigs (A, k, sigma)

d = eigs (A, k, sigma, opts)

d = eigs (A, B)

d = eigs (A, B, k)

d = eigs (A, B, k, sigma)

d = eigs (A, B, k, sigma, opts)

d = eigs (af, n)

d = eigs (af, n, B)

d = eigs (af, n, k)

d = eigs (af, n, B, k)

d = eigs (af, n, k, sigma)

d = eigs (af, n, B, k, sigma)

d = eigs (af, n, k, sigma, opts)

d = eigs (af, n, B, k, sigma, opts)

[V, d] = eigs (A, …)

[V, d] = eigs (af, n, …)

[V, d, flag] = eigs (A, …)

[V, d, flag] = eigs (af, n, …)

Вычисление ограниченного числа собственных значений и собственных векторов матрицы A на основе критериев выбора.

Количество вычисляемых собственных значений и векторов задаётся параметром k и по умолчанию равно 6.

По умолчанию, eigs решает уравнение, где — соответствующий собственный вектор. Если задана положительно определённая матрица B, то eigs решает общее уравнение для собственных значений.

Аргумент sigma определяет, какие собственные значения будут возвращены. sigma может быть скалярным значением или строкой. Если sigma является скалярным, возвращаются k собственных значений, наиболее близких к sigma. Если sigma является строкой, она должна иметь одно из следующих значений.

"lm"

Наибольший модуль (по умолчанию).

"sm"

Наименьший модуль.

"la"

Наибольшая алгебраическая часть (действительно только для вещественно симметричных задач).

"sa"

Наименьшая алгебраическая часть (действительно только для вещественно симметричных задач).

"be"

Оба конца, с одним дополнительным значением с высокой стороны, если k нечётное (действительно только для вещественно симметричных задач).

"lr"

Наибольшая вещественная часть (действительно только для комплексных или несимметричных задач).

"sr"

Наименьшая вещественная часть (действительно только для комплексных или несимметричных задач).

"li"

Наибольшая мнимая часть (действительно только для комплексных или несимметричных задач).

"si"

Наименьшая мнимая часть (действительно только для комплексных или несимметричных задач).

Если opts задан, это структура, определяющая возможные опции, которые eigs должен использовать. Поля структуры opts:

issym

Если af задан, флаг, указывающий, определяет ли функция af симметричную задачу. Игнорируется, если задана A. По умолчанию — ложь.

isreal

Если af задан, флаг, указывающий, определяет ли функция af вещественную задачу. Игнорируется, если задана A. По умолчанию — истина.

tol

Определяет требуемую толерантность сходимости, вычисляемую как tol * norm (A). По умолчанию — eps.

maxit

Максимальное количество итераций. По умолчанию — 300.

p

Количество векторов базиса Ланцоша для использования. Большее количество векторов приводит к более быстрой сходимости, но и к большему потреблению памяти. Оптимальное значение p зависит от задачи и должно находиться в диапазоне k + 1 до n. Значение по умолчанию — 2 * k.

v0

Начальный вектор для алгоритма. Начальный вектор, близкий к конечному вектору, ускорит сходимость. По умолчанию ARPACK генерирует начальный вектор случайным образом. Если указан, v0 должен быть n-строчным вектором, где n = rows (A)

disp

Уровень диагностических выводов (0|1|2). Если disp равно 0, диагностика отключена. Значение по умолчанию — 0.

cholB

Флаг, если chol (B) передаётся вместо B. По умолчанию — ложь.

permB

Вектор перестановки факторизации Холецкого B, если cholB равно истине. Получается с помощью [R, ~, permB] = chol (B, "vector"). По умолчанию — 1:n.

Также возможно представить A с помощью функции, обозначенной af. af должна быть сопровождаема скалярным аргументом n, определяющим длину вектора-аргумента, принятого af. af может быть функцией-обработчиком, встроенной функцией или строкой. Если af является строкой, она содержит имя функции для использования.

af является функцией вида y = af (x), где требуемое возвращаемое значение af определяется значением sigma. Четыре возможные формы:

A * x

если sigma не задан или является строкой, отличной от "sm".

A \ x

если sigma равен 0 или "sm".

(A - sigma * I) \ x

для стандартной задачи на собственные значения, где I — единичная матрица такого же размера, что и A.

(A - sigma * B) \ x

для общей задачи на собственные значения.

Возвращаемые значения eigs зависят от количества запрашиваемых возвращаемых значений. При одном возвращаемом значении возвращается вектор d длины k содержащий k найденных собственных значений. При двух возвращаемых значениях V — матрица n-на-k, столбцы которой — k собственных векторов, соответствующих возвращённым собственным значениям. Сами собственные значения возвращаются в d в виде n-на-k матрицы, где элементы на диагонали — собственные значения.

При третьем возвращаемом аргументе flag, eigs возвращает статус сходимости. Если flag равен 0, все собственные значения сошлись. Любое другое значение указывает на неудачную сходимость.

Эта функция основана на пакете ARPACK, написанном R. Lehoucq, K. Maschhoff, D. Sorensen, и C. Yang. Более подробная информация доступна по адресу http://www.caam.rice.edu/software/ARPACK/.

См. также: eig, svds.

s = svds (A)

s = svds (A, k)

s = svds (A, k, sigma)

s = svds (A, k, sigma, opts)

[u, s, v] = svds (…)

[u, s, v, flag] = svds (…)

Нахождение нескольких сингулярных значений матрицы A.

Сингулярные значения вычисляются с помощью

[m, n] = size (A);
s = eigs ([sparse(m, m), A;
                     A', sparse(n, n)])

Собственные значения, возвращаемые eigs, соответствуют сингулярным значениям A. Количество сингулярных значений для вычисления задаётся параметром k и по умолчанию равно 6.

Аргумент sigma определяет, какие сингулярные значения нужно найти. Когда sigma — строка 'L', по умолчанию, находятся наибольшие сингулярные значения A. В противном случае sigma должен быть действительным скаляром, и находятся сингулярные значения, наиболее близкие к sigma. Соответственно, sigma = 0 находит наименьшие сингулярные значения. Обратите внимание, что при сравнительно малых значениях sigma существует вероятность, что нужное количество сингулярных значений не будет найдено. В таком случае sigma следует увеличить.

opts — структура, определяющая опции, которые svds передаст eigs. Возможные поля этой структуры документированы в eigs. По умолчанию svds устанавливает следующие три поля:

tol

Требуемая толерантность сходимости для сингулярных значений. Значение по умолчанию — 1e-10. eigs передаётся tol / sqrt(2).

maxit

Максимальное количество итераций. По умолчанию — 300.

disp

Уровень диагностических выводов (0|1|2). Если disp равно 0, диагностика отключена. Значение по умолчанию — 0.

Если запрашивается более одного выходного значения, то svds вернёт приближение сингулярного разложения A

A_approx = u*s*v'

где A_approx — матрица такого же размера, что и A, но только ранга k.

flag возвращает 0, если алгоритм успешно сошёлся, и 1 в противном случае. Критерий сходимости:

norm (A*v - u*s, 1) <= tol * norm (A, 1)

svds лучше всего подходит для нахождения только нескольких сингулярных значений большой разреженной матрицы. В противном случае svd (full (A)) будет, скорее всего, более эффективным.

См. также: svd, eigs.