Иерархическая матрица - Hierarchical matrix

В вычислительная математика, иерархические матрицы (H-матрицы)^[1][2][3]используются в качестве приближений нерезких матриц с разреженными данными. разреженная матрица измерения ${ displaystyle n}$ может быть эффективно представлена ​​в ${ Displaystyle О (п)}$ единиц хранения, сохраняя только свои ненулевые элементы, неразреженная матрица потребует ${ Displaystyle О (п ^ {2})}$ единиц хранения, поэтому использование этого типа матриц для больших задач было бы чрезмерно дорогим с точки зрения хранения и времени вычислений. Иерархические матрицы обеспечивают приближение, требующее только ${ Displaystyle О (пк , журнал (п))}$ единиц хранения, где ${ displaystyle k}$ - параметр, контролирующий точность аппроксимации. В типичных приложениях, например, при дискретизации интегральных уравнений^{[4][5][6][7][8]},^[9]предварительное кондиционирование полученных систем линейных уравнений,^[10]или решение эллиптических уравнений в частных производных^{[11][12][13][14]}, ранг, пропорциональный ${ displaystyle log (1 / epsilon) ^ { gamma}}$ с небольшой постоянной ${ displaystyle gamma}$ достаточно, чтобы обеспечить точность ${ displaystyle epsilon}$По сравнению со многими другими представлениями неразреженных матриц с разреженными данными, иерархические матрицы предлагают главное преимущество: результаты матричных арифметических операций, таких как матричное умножение, факторизация или инверсия, могут быть аппроксимированы в ${ Displaystyle О (пк ^ { альфа} , журнал (п) ^ { бета})}$ операции, где ${ Displaystyle альфа, бета в {1,2,3 }.}$^[2]

Основная идея

Иерархические матрицы полагаются на локальные аппроксимации низкого ранга: пусть ${ displaystyle I, J}$ быть индексными наборами, и пусть ${ displaystyle G in { mathbb {R}} ^ {I times J}}$ обозначают матрицу, которую мы должны аппроксимировать. Во многих приложениях (см. выше) мы можем найти подмножества ${ displaystyle t substeq I, s substeq J}$ такой, что ${ Displaystyle G | _ {т раз s}}$можно аппроксимировать рангом${ displaystyle k}$ Это приближение может быть представлено в факторизованном виде ${ displaystyle G | _ {t times s} приблизительно AB ^ {*}}$ с факторами${ displaystyle A in { mathbb {R}} ^ {t times k}, B in { mathbb {R}} ^ {s times k}}$.При стандартном представлении матрицы ${ Displaystyle G | _ {т раз s}}$ требует ${ Displaystyle О (( # т) ( # s))}$ единиц хранения, факторизованное представление требует только ${ Displaystyle О (к ( # т + # s))}$ ед., если ${ displaystyle k}$ не слишком большой, требования к хранилищу значительно снижаются.

Чтобы аппроксимировать всю матрицу ${ displaystyle G}$, он разбивается на семейство подматриц. Большие подматрицы хранятся в факторизованном представлении, а маленькие подматрицы хранятся в стандартном представлении для повышения эффективности.

Матрицы низкого ранга тесно связаны с вырожденными разложениями, используемыми в группировка панелей и быстрый мультипольный метод для аппроксимации интегральных операторов. В этом смысле иерархические матрицы можно рассматривать как алгебраические аналоги этих методов.

Приложение к интегральным операторам

Иерархические матрицы успешно используются для обработки интегральных уравнений, например, операторов потенциала одинарного и двойного слоя, появляющихся в метод граничных элементов Типичный оператор имеет вид

${ Displaystyle { mathcal {G}} [u] (x) = int _ { Omega} kappa (x, y) u (y) , dy.}$

В Метод Галеркина приводит к матричным элементам вида

${ displaystyle g_ {ij} = int _ { Omega} int _ { Omega} kappa (x, y) varphi _ {i} (x) psi _ {j} (y) , dy , dx,}$

куда ${ Displaystyle ( varphi _ {я}) _ {я in I}}$ и ${ Displaystyle ( psi _ {j}) _ {j in J}}$ являются семействами конечно-элементных базисных функций. Если функция ядра ${ displaystyle kappa}$ достаточно гладкий, мы можем аппроксимировать его полиномиальная интерполяция чтобы получить

${ Displaystyle { тильда { каппа}} (х, у) = сумма _ { ню = 1} ^ {к} каппа (х, хи _ { ню}) ell _ { nu} (y),}$

куда ${ Displaystyle ( хи _ { ню}) _ { ню = 1} ^ {к}}$ - семейство точек интерполяции и ${ displaystyle ( ell _ { nu}) _ { nu = 1} ^ {k}}$соответствующее семейство Полиномы Лагранжа.Замена ${ displaystyle kappa}$ к ${ displaystyle { tilde { kappa}}}$ дает приближение

${ displaystyle { tilde {g}} _ {ij} = int _ { Omega} int _ { Omega} { tilde { kappa}} (x, y) varphi _ {i} (x ) psi _ {j} (y) , dy , dx = sum _ { nu = 1} ^ {k} int _ { Omega} kappa (x, xi _ { nu}) varphi _ {i} (x) , dx int _ { Omega} ell _ { nu} (y) psi _ {j} (y) , dy = sum _ { nu = 1 } ^ {k} a_ {i nu} b_ {j nu}}$

с коэффициентами

${ Displaystyle a_ {я nu} = int _ { Omega} kappa (x, xi _ { nu}) varphi _ {i} (x) , dx,}$

${ displaystyle b_ {j nu} = int _ { Omega} ell _ { nu} (y) psi _ {j} (y) , dy.}$

Если мы выберем ${ displaystyle t substeq I, s substeq J}$ и использовать одни и те же точки интерполяции для всех ${ displaystyle i in t, j in s}$, мы получаем${ displaystyle G | _ {t times s} приблизительно AB ^ {*}}$.

Очевидно, что любое другое приближение, разделяющее переменные ${ displaystyle x}$ и ${ displaystyle y}$например, мультипольное разложение, также позволило бы нам разделить двойной интеграл на два одиночных интеграла и, таким образом, прийти к аналогичной факторизованной матрице низкого ранга.

Особый интерес представляют методы перекрестной аппроксимации.^[6][7][15]которые используют только элементы исходной матрицы ${ displaystyle G}$ построить приближение низкого ранга.

Приложение к эллиптическим уравнениям в частных производных

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

Функция Грина зависит от формы вычислительной области, поэтому она обычно не известна. Тем не менее, приближенные арифметические операции могут использоваться для вычисления приближенного обратного значения без явного знания функции.

Удивительно, но доказать можно^{[11][12][13][14]} что обратная величина может быть аппроксимирована, даже если дифференциальный оператор включает негладкие коэффициенты, и поэтому функция Грина негладкая.

Арифметические операции

Важнейшим нововведением метода иерархических матриц является разработка эффективных алгоритмов для выполнения (приближенных) матричных арифметических операций над не разреженными матрицами, например, для вычисления приближенных обратных функций, LU разложения и решения матричных уравнений.

Центральный алгоритм - эффективное умножение матрицы на матрицу, т. Е. Вычисление ${ Displaystyle Z = Z + альфа XY}$для иерархических матриц ${ displaystyle X, Y, Z}$ и скалярный множитель ${ displaystyle alpha}$Алгоритм требует, чтобы подматрицы иерархических матриц были организованы в блочную древовидную структуру, и использует свойства факторизованных матриц низкого ранга для вычисления обновленных ${ displaystyle Z}$ в${ Displaystyle О (пк ^ {2} , журнал (п) ^ {2})}$ операции.

Воспользовавшись блочной структурой, обратное можно вычислить с помощью рекурсии для вычисления обратных иШур дополняет диагональных блоков и комбинируя их с помощью умножения матрицы на матрицу. LU разложение^[16][17]могут быть построены с использованием только рекурсии и умножения. Обе операции также требуют ${ Displaystyle О (пк ^ {2} , журнал (п) ^ {2})}$ операции.

ЧАС²-матрицы

Для решения очень больших задач структура иерархических матриц может быть улучшена: H²-матрицы^[18][19]заменить общую низкоранговую структуру блоков иерархическим представлением, тесно связанным сбыстрый мультипольный метод чтобы уменьшить сложность хранения до ${ Displaystyle О (нк)}$.

В контексте граничных интегральных операторов замена фиксированного ранга ${ displaystyle k}$ блок-зависимыми рангами приводит к приближениям, которые сохраняют скорость сходимости лежащего в основе метода граничных элементов со сложностью ${ Displaystyle O (п).}$^[20][21]

Арифметические операции, такие как умножение, инверсия и факторизация Холецкого или LR H²-матрицы могут быть реализованы на основе двух основных операций: умножения матрицы на вектор на подматрицы и обновления подматриц с низким рангом. Хотя умножение матрицы на вектор является простым, реализация эффективных обновлений с низким рангом с адаптивно оптимизированными базами кластеров представляет собой значительную проблему.^[22]

Литература

Программного обеспечения

HLib - это программная библиотека на языке C, реализующая наиболее важные алгоритмы иерархических и ${ Displaystyle { mathcal {H}} ^ {2}}$-матрицы.

АХМЕД - это программная библиотека C ++, которую можно загрузить в образовательных целях.

HLIBpro представляет собой реализацию основных алгоритмов иерархической матрицы для коммерческих приложений.

H2Lib представляет собой реализацию алгоритмов иерархической матрицы с открытым исходным кодом, предназначенную для исследований и обучения.

удивительные иерархические матрицы - это репозиторий, содержащий список других реализаций H-матриц.

HierarchicalMatrices.jl представляет собой пакет Julia, реализующий иерархические матрицы.