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"
• 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 (…)

Оцените числовую обусловленность квадратной матрицы 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.

Ссылки:

• 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, 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. По умолчанию false.

isreal

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

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. По умолчанию 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

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

(A - sigma * B) \ x

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

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

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

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