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 (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"
• истинно, если 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.

Ссылка: N. J. Higham и F. Tisseur, Блочный алгоритм для оценки 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, 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 (…) ¶

Оценивает числовую обусловленность квадратной матрицы A по норме 1, используя 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. Поэтому, если требуются согласованные результаты, состояние генератора случайных чисел должно быть зафиксировано перед вызовом condest.

Ссылки:

• N.J. Higham and 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 and 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, написанном Р. Лехуком, К. Машхоффом, Д. Соресеном и Ц. Янгом. Дополнительную информацию см. на странице 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.