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

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

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

1. Если матрица диагональная, решить напрямую и перейти к 8
2. Если матрица является переставленной диагональной, решить напрямую с учетом перестановок. Перейти к 8
3. Если матрица квадратная, полосовая и если плотность полосы меньше, чем заданная spparms ("bandden"), продолжить, иначе перейти к 4.
   1. Если матрица треугольная и правая часть не разреженная, продолжить, иначе перейти к 3б.
      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 (Afcn, 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 используется функция Afcn (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, … являются аргументами Afcn (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.

Ссылка: N. J. Higham и F. Tisseur, A block algorithm for matrix 1-norm estimation, with and application to 1-norm pseudospectra, SIAM J. Matrix Anal. Appl., pp. 1185–1201, Vol 21, No. 4, 2000.

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

: cest = condest (A) ¶

: cest = condest (A, t) ¶

: cest = condest (A, Ainvfcn) ¶

: cest = condest (A, Ainvfcn, t) ¶

: cest = condest (A, Ainvfcn, t, p1, p2, …) ¶

: cest = condest (Afcn, Ainvfcn) ¶

: cest = condest (Afcn, Ainvfcn, t) ¶

: cest = condest (Afcn, Ainvfcn, t, p1, p2, …) ¶

: [cest, v] = condest (…) ¶

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

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

Входными данными может быть матрица A (алгоритм особенно подходит для больших разреженных матриц). В качестве альтернативы поведение матрицы может быть определено неявно с помощью функций. При использовании неявного определения condest требует следующих функций:

• Afcn (flag, x), которая должна возвращать
  • размерность n A, если flag равен "dim"
  • true, если A является вещественным оператором, если flag равен "real"
  • результат A * x, если flag равен "notransp"
  • результат A' * x, если flag равен "transp"
• Ainvfcn (flag, x), которая должна возвращать
  • размерность n inv (A), если flag равен "dim"
  • true, если inv (A) является вещественным оператором, если flag равен "real"
  • результат inv (A) * x, если flag равен "notransp"
  • результат inv (A)' * x, если flag равен "transp"

Любые параметры p1, p2, … являются дополнительными аргументами Afcn (flag, x, p1, p2, …) и Ainvfcn (flag, x, p1, p2, …).

Основным результатом является оценка числа обусловленности 1-нормы cest.

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

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

Ссылки:

• N.J. Higham и F. Tisseur, A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra. SIMAX vol 21, no 4, pp 1185–1201. https://dx.doi.org/10.1137/S0895479899356080
• N.J. Higham и F. Tisseur, A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra. https://citeseer.ist.psu.edu/223007.html

См. также: cond, rcond, 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, k) ¶

: d = eigs (Af, n, k, sigma) ¶

: d = eigs (Af, n, k, sigma, opts) ¶

: d = eigs (Af, n, B) ¶

: d = eigs (Af, n, B, k) ¶

: d = eigs (Af, n, B, k, sigma) ¶

: d = eigs (Af, n, B, k, sigma, opts) ¶

: [V, D] = eigs (…) ¶

: [V, D, flag] = eigs (…) ¶

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

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

Входная A — это квадратная матрица размера n-на-n. Как правило, A также является большой и разреженной.

Входная B для обобщённой задачи на собственные значения — это квадратная матрица того же размера, что и A (n-на-n). Обычно B также является большой и разреженной.

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

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

"lm"

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

"sm"

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

"la"

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

"sa"

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

"be"

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

"lr"

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

"sr"

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

"li"

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

"si"

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

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

issym

Если Af задан, этот флаг (true/false) определяет, определяет ли функция Af симметричную задачу. Он игнорируется, если задана матрица A. По умолчанию false.

isreal

Если Af задан, этот флаг (true/false) определяет, определяет ли функция Af вещественную задачу. Он игнорируется, если задана матрица A. По умолчанию true.

tol

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

maxit

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

p

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

v0

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

disp

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

cholB

Если вычисляется обобщённая задача на собственные значения, этот флаг (true/false) указывает, представляет ли вход B chol (B) или просто матрицу B. По умолчанию false.

permB

Вектор перестановок факторизации Холецкого для B, если cholB равно true. Он получен [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

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

(A - sigma * B) \ x

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

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

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

Замечания по программированию: для малых задач, n < 500, рассмотрите использование eig (full (A)).

Если ARPACK не сходится, рассмотрите увеличение количества векторов Ланцоша (opt.p), увеличение числа итераций (opt.maxiter) или уменьшение толерантности (opt.tol).

Ссылка: Эта функция основана на пакете 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.