Математика / Интерполяция

Алгоритм AAA: адаптивное рациональное приближение

Алгоритм AAA строит рациональное приближение табличных данных в барицентрической форме. Он постепенно выбирает опорные узлы с наибольшей текущей ошибкой, а затем находит веса по задаче наименьших квадратов.

Опубликовано: Обновлено:

Формула

$$r(z)=\frac{\displaystyle\sum_{j=1}^{m}\frac{w_j f_j}{z-z_j}}{\displaystyle\sum_{j=1}^{m}\frac{w_j}{z-z_j}}$$
Схема Барицентрическая рациональная модель
1234500.20.40.60.81r(z)zAAA-форма с заданными весами

Опоры: (−1;1) и (1;0), оба веса равны 1. На показанном участке модель r(x)=(x−1)/(2x).

Обозначения

$z$
точка вычисления, единица аргумента
$z_j$
j-й опорный узел, единица аргумента
$f_j$
значение данных в опорном узле, единица функции
$w_j$
ненулевой барицентрический вес, —
$m$
число выбранных опорных узлов, —
$r(z)$
значение рационального приближения, единица функции

Условия применения

Узлы z_j должны быть различны; значения f_j должны быть заданы в одной и той же шкале. Для полного алгоритма задают исходный набор пробных узлов, выбирают очередную опору с наибольшим модулем остатка и решают задачу наименьших квадратов для новых весов на оставшихся пробах при единичной евклидовой норме вектора весов. В точной опорной точке формулу понимают по непрерывному продолжению: r(z_j)=f_j, если соответствующий вес ненулевой. В остальных точках знаменатель не должен обращаться в ноль.

Ограничения

Рациональное приближение может иметь полюса, и нули знаменателя дают кандидатов в полюса, если одновременно не зануляется числитель; вычисление рядом с полюсом чувствительно. AAA стремится уменьшать ошибку на предоставленном наборе, но это не гарантия точности между пробами или вне их области. Шумные данные могут провоцировать лишние опоры и паразитные пары нуль–полюс, поэтому проверяют невязку на независимых точках и при необходимости применяют обрезку таких пар. Метод не является линейным интерполянтом через все исходные узлы: значения совпадают в выбранных опорах с ненулевыми весами.

Подробное объяснение

Барицентрическая форма рациональной функции записывает отдельно сумму вкладов значений и сумму вкладов весов. У каждого слагаемого один и тот же множитель 1/(z−z_j), поэтому в опорном узле предел отношения равен соответствующему f_j. Если опор две, отношение упрощается до рациональной функции низкой степени; с большим количеством опор доступна более гибкая форма. В AAA множество опор не фиксируют заранее: после первого узла выбирают ещё одну пробу, где текущая модель ошибается сильнее всего.

При добавлении опоры веса подбирают из линейной задачи наименьших квадратов, сформированной по непринятым пока пробам. Процесс повторяют до достижения допуска или заданного числа опор. Поэтому алгоритм совмещает две идеи: рациональное барицентрическое представление и жадный выбор опорных узлов по максимальному остатку. В опубликованной статье авторы называют метод adaptive Antoulas–Anderson в честь использованных ими идей Антоуласа и Андерсона.

Для ручной проверки уже найденной модели вычислите числитель и знаменатель отдельно, сохраняя дроби до финального деления. Например, при опорах z₀=−1 и z₁=1, значениях 1 и 0, весах 1 и 1 формула равна [1/(z+1)]/[1/(z+1)+1/(z−1)]=(z−1)/(2z). При z=2 получаем 1/4. Этот расчёт иллюстрирует оценку из фиксированной AAA-формы; он не воспроизводит шаг выбора опор и решения SVD, который нужен для построения весов по большому набору данных.

В отличие от формулы Флоатера—Хорманна, где локальный параметр задаёт степень и веса для интерполяции всей сетки, AAA выбирает небольшое число опор глобально и обычно приближает оставшиеся узлы, минимизируя дискретную невязку. Такой подход может быть компактнее, но результат зависит от набора проб и возможных полюсов. Практическая проверка включает график, независимые точки и анализ нулей и полюсов знаменателя.

На пробе u вне опор умножают остаток на барицентрический знаменатель. Получается линейная комбинация весов с коэффициентами (f(u)−fⱼ)/(u−zⱼ). Из этих строк строят матрицу Лёвнера, а единичный вектор весов берут из правых сингулярных векторов для наименьшего сингулярного числа.

Как пользоваться формулой

  1. Сформируйте набор пробных узлов и значений; проверьте, что масштаб переменной и данных удобен для вычислений.
  2. Выберите начальную опору, затем вычислите остаток модели на всех ещё не выбранных пробах.
  3. Добавляйте узел с наибольшим модулем остатка и находите веса из задачи наименьших квадратов для оставшихся точек.
  4. Остановитесь по заранее выбранному допуску или числу опор и сохраните узлы, значения и веса.
  5. Вычисляйте r(z) барицентрической формулой; отдельно проверьте знаменатель, опорные точки и независимые контрольные значения.

Историческая справка

AAA был предложен Юдзи Накацукасой, Оливье Сете и Ллойдом Трефетеном в статье *The AAA Algorithm for Rational Approximation*. Рукопись поступила в SIAM в декабре 2016 года, была принята в декабре 2017 года и опубликована онлайн в мае 2018 года в *SIAM Journal on Scientific Computing*. Авторы объединили барицентрическое представление рациональной функции с жадным выбором опор, что позволило строить приближения по вещественным и комплексным наборам проб без заранее заданных параметров типа.

Название расшифровывается как adaptive Antoulas–Anderson: авторы прямо указали, что оно отдаёт должное Димитрису Антоуласа и Брайану Андерсону за идеи, на которых основан формат. Это название конкретного алгоритма, а не утверждение о создании всей рациональной интерполяции. В статье 2018 года описаны вычислительный шаг, MATLAB-реализация и примеры применения, включая редукцию моделей и идентификацию систем. Поэтому для исторической атрибуции достаточно указывать авторов оригинальной статьи, дату и названную ими связь с предшествующей работой.

Историческая линия формулы

Алгоритм AAA опубликован Юдзи Накацукасой, Оливье Сете и Ллойдом Трефетеном в 2018 году. Название adaptive Antoulas–Anderson отсылает к Димитрису Антоулаcу и Брайану Андерсону; статья объединяет барицентрическую рациональную форму с жадным выбором опор.

Пример

Пусть после адаптивного выбора остались две опоры z₀=−1, z₁=1, значения f₀=1, f₁=0 и веса w₀=w₁=1. Найдём значение рациональной модели при z=2. Числитель равен 1·1/(2−(−1))+1·0/(2−1)=1/3. Знаменатель равен 1/(2−(−1))+1/(2−1)=1/3+1=4/3. Поэтому r(2)=(1/3)/(4/3)=1/4. Проверка упрощением: общая формула сокращается до (z−1)/(2z), что при z=2 также даёт 1/4. Значения в опорах сохраняются: r(−1)=1 и r(1)=0 по пределу исходного отношения. Это пример вычисления по готовым весам, а не ручного выполнения жадного выбора опор: веса в реальной задаче получают из набора проб и остатков.

Частая ошибка

Не подставляйте z=z_j напрямую в числитель и знаменатель: там возникает 0/0, а значение берут по пределу. Не считайте все исходные узлы интерполяционными опорами: AAA обычно совпадает с данными только в выбранном подмножестве и приближает остальные. Наконец, малый остаток на обучающих пробах не исключает полюс или большой выброс между ними; проверяйте знаменатель и отдельные точки, которых не было при подборе весов.

Практика

Задачи с решением

Выбрать первую опору

Условие. На узлах 0, 1, 2 даны значения 0, 3, 9. Начальное приближение AAA — среднее значение. Какую пробу берут первой?

Решение. Начальное среднее (0+3+9)/3=4. Остатки равны −4, −1 и 5, поэтому максимальный модуль 5 достигается при z=2. После добавления единственной опоры барицентрическое отношение равно её значению 9 при ненулевом весе.

Ответ. Первая опора z=2, одноопорная модель r=9.

Отыскать полюс по весам

Условие. Даны опоры 0 и 2, значения 1 и 4, веса 2 и −1. Где у рациональной модели полюс?

Решение. После умножения числителя и знаменателя на z(z−2) получаем r(z)=[2(z−2)−4z]/[2(z−2)−z]=(−2z−4)/(z−4). При z=4 числитель равен −12 и сокращения нет.

Ответ. Простой полюс z=4.

Дополнительные источники

  • https://epubs.siam.org/doi/abs/10.1137/16M1106122
  • https://people.maths.ox.ac.uk/trefethen/AAAfinal.pdf

Связанные формулы

Математика

Рациональная интерполяция Флоатера — Хорманна при d = 1

$\begin{aligned}r(x)&=\frac{\sum_{j=0}^n w_jy_j/(x-x_j)}{\sum_{j=0}^n w_j/(x-x_j)},\\ w_j&=(-1)^{j-1}\left(\frac{\mathbf1_{j>0}}{h_{j-1}}+\frac{\mathbf1_{j<n}}{h_j}\right).\end{aligned}$

При d = 1 рациональный интерполянт Флоатера — Хорманна сочетает соседние линейные полиномы. На равномерной сетке внутренние веса вдвое больше крайних по модулю.

Математика

Барицентрическая интерполяция: веса и вычисление полинома

$P(x)=\frac{\sum_{k=0}^{n}\dfrac{w_ky_k}{x-x_k}}{\sum_{k=0}^{n}\dfrac{w_k}{x-x_k}},\qquad w_k=\frac{1}{\prod_{j\ne k}(x_k-x_j)},\quad x\ne x_k$

Барицентрическая форма вычисляет полином Лагранжа через заранее рассчитанные веса узлов. Она удобна для множества запросов по одной таблице; при точном совпадении аргумента с узлом нужно сразу вернуть узловое значение.

Математика

Интерполяция Лагранжа: полином через заданные точки

$P_n(x)=\sum_{k=0}^{n}y_kL_k(x),\qquad L_k(x)=\prod_{\substack{0\leq j\leq n\\j\ne k}}\frac{x-x_j}{x_k-x_j}$

Формула Лагранжа строит полином степени не выше n по n+1 различным узлам. Каждый базисный множитель выбирает значение своего узла и обращает вклады остальных узлов в ноль.

Математика

Константа Лебега: усиление ошибки интерполяции

$L_n(x)=\sum_{j=0}^n|\ell_j(x)|,\qquad\Lambda_n=\max_{x\in[a,b]}L_n(x)$

Константа Лебега измеряет наибольшее усиление ошибок значений при полиномиальной интерполяции на фиксированной сетке. Она зависит от расположения узлов и выбранного отрезка, но не от конкретных значений функции.