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 (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(1) аргумента iter содержится число итераций (максимальное значение 5 задано в коде), а в iter(2) — общее число произведений A * x или A' * x, выполненных алгоритмом.

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

Ссылка: 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 следует зафиксировать "state" генератора случайных чисел.

Литература:

• 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’

Флаг, указывающий, будут ли для операций LU, ’\’ и ’/’ использоваться решатели UMFPACK или mmd (по умолчанию 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.