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) ¶

...

: 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 с использованием 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, Алгоритм блока для оценки 1-нормы матрицы, с приложением к 1-норме псевдоспектра. SIMAX том 21, № 4, стр. 1185–1201. https://dx.doi.org/10.1137/S0895479899356080
• N.J. Higham и F. Tisseur, Алгоритм блока для оценки 1-нормы матрицы, с приложением к 1-норме псевдоспектра. 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.