Лицензия MIT
© Алексей Александрович Щербаков, 2024
В настоящее время обратные и оптимизационные задачи в англоязчной литературе называются методами обратного проектирования. Решение таких задач происходит итеративно, а значит, должно включать в себя математическую модель и соответствующий численный метод для решения прямой задачи, а также оптимизационный алгоритм который позволяет прийти к нужным параметрам задачи в соответствии с заданной пользователем целевой производительностью рассматриваемого устройства.
Методы решения обратных задач можно разделить на методы, использующие информацию о градиенте целевой функции и не использующие. Градиентная оптимизация является наиболее широко используемым классом методов решения обратных задач фотоники. Она представляет собой итеративный алгоритм, в котором на каждой итерации рассчитывается вектор градиента целевой функции относительно всех параметров задачи, затем эти параметры изменяются в направлении градиента. Если параметров много, а прямой метод достаточно трудоёмкий, то вычисление частных производных по всем параметрам через конечные разности представляет требует чрезмерных вычислительных затрат. Для решения этой проблемы применяется так называемый метод сопряжения (adjoint method), требующий решения всего двух прямых задач для вычисления полного вектора градиента независимо от числа параметров.
В данной лекции дан краткий обзор градиентных методов и приведены формулировки метода сопряжения для линейных и нелинейных задач.
Рассмотрим линейную задачу. Пусть волновые уравнения решаются численным методом, сводящимся к системе линейных алгебраических уравнений: \begin{equation} M \boldsymbol{a} = \boldsymbol{a}_{inc} \end{equation} Обозначим вектор оптимизационных параметров как $\boldsymbol{p}$, а целевую функцию как $\mathcal{F}=\mathcal{F}(\boldsymbol{a}(\boldsymbol{p}))$. Будем считать целевую функцию аналитически известной функцией вектора $\boldsymbol{a}$, так что производные $\partial\mathcal{F}/\partial\boldsymbol{a}$ легко вычислимы. Рассмотрим частную производную целевой функций по одному из параметров: \begin{equation} \frac{\partial\mathcal{F}}{\partial p_j} = \frac{\partial\mathcal{F}}{\partial\boldsymbol{a}} \frac{\partial\boldsymbol{a}}{\partial p_j} \end{equation} Для вычисления последнего множителя рассмотрим исходную прямую задачу и продифференцируем её по рассматриваемому параметру: \begin{equation} \frac{\partial M}{\partial p_j} \boldsymbol{a} + M \frac{\partial\boldsymbol{a}}{\partial p_j} = \frac{\partial\boldsymbol{a}_{inc}}{\partial p_j} \; \Rightarrow \; \frac{\partial\boldsymbol{a}}{\partial p_j} = M^{-1} \left( \frac{\partial\boldsymbol{a}_{inc}}{\partial p_j} - \frac{\partial M}{\partial p_j} \boldsymbol{a} \right) \end{equation} Тода \begin{equation} \frac{\partial\mathcal{F}}{\partial p_j} = \left( \frac{\partial\mathcal{F}}{\partial\boldsymbol{a}} M^{-1} \right) \left( \frac{\partial\boldsymbol{a}_{inc}}{\partial p_j} - \frac{\partial M}{\partial p_j} \boldsymbol{a} \right) = \boldsymbol{a}^T_{adj} \left( \frac{\partial\boldsymbol{a}_{inc}}{\partial p_j} - \frac{\partial M}{\partial p_j} \boldsymbol{a} \right) \end{equation} Здесь $\boldsymbol{a}_{adj}$ - решение сопряженной задачи \begin{equation} M^T \boldsymbol{a}_{adj} = \left( \frac{\partial\mathcal{F}}{\partial\boldsymbol{a}} \right)^T \end{equation} Зачастую, наиболее вычислительно сложной частью оптмизационного процесса является решение прямой задачи, и метод сопряжения позволяет вычислять вектор градиента по заданному вектору параметров через решение всего двух прямых задач.
В более общем случае нелинейной задачи решение удовлетворяет некоторому уравнению \begin{equation} \boldsymbol{\varphi} \left( \boldsymbol{a}, \boldsymbol{a}_{inc}, \boldsymbol{p} \right) = 0 \end{equation} Его дифференцирование даёт \begin{equation} \frac{\partial\boldsymbol{\varphi}}{\partial\boldsymbol{a}} \frac{\partial\boldsymbol{a}}{\partial p_j} + \frac{\partial\boldsymbol{\varphi}}{\partial\boldsymbol{a}_{inc}} \frac{\partial\boldsymbol{a}_{inc}}{\partial p_j} + \frac{\partial\boldsymbol{\varphi}}{\partial p_j} = 0 \Rightarrow \frac{\partial\boldsymbol{a}}{\partial p_j} = -\left( \frac{\partial\boldsymbol{\varphi}}{\partial\boldsymbol{a}} \right)^{-1} \left( \frac{\partial\boldsymbol{\varphi}}{\partial\boldsymbol{a}_{inc}} \frac{\partial\boldsymbol{a}_{inc}}{\partial p_j} + \frac{\partial\boldsymbol{\varphi}}{\partial p_j} \right) \end{equation} Тогда для компоненты градиента по параметру $p_j$ имеем \begin{equation} \frac{\partial\mathcal{F}}{\partial p_j} = \frac{\partial\mathcal{F}}{\partial\boldsymbol{a}} \frac{\partial\boldsymbol{a}}{\partial p_j} = -\left[ \frac{\partial\mathcal{F}}{\partial\boldsymbol{a}} \left( \frac{\partial\boldsymbol{\varphi}}{\partial\boldsymbol{a}} \right)^{-1} \right] \left( \frac{\partial\boldsymbol{\varphi}}{\partial\boldsymbol{a}_{inc}} \frac{\partial\boldsymbol{a}_{inc}}{\partial p_j} + \frac{\partial\boldsymbol{\varphi}}{\partial p_j} \right) = - \boldsymbol{\chi} \left( \frac{\partial\boldsymbol{\varphi}}{\partial\boldsymbol{a}_{inc}} \frac{\partial\boldsymbol{a}_{inc}}{\partial p_j} + \frac{\partial\boldsymbol{\varphi}}{\partial p_j} \right) \end{equation} Основным отличием от линейного случая является то, что уравнение на $\boldsymbol{\chi}$ не является сопряженным исходному нелинейному уравнению, но сложность его решения, как правило, меньше сложность исходного, поскольку оно сводится к умножению обратной матрицы на заданный вектор.
В качестве примера рассмотрим задачу оптимизации фнуконала $\mathcal{F}(\boldsymbol{a},\lambda,\boldsymbol{p})$ для задачи на собственные значения \begin{equation} M \boldsymbol{a} = \lambda \boldsymbol{a} \end{equation} Для простоты будем считать матрицу системы действиетельной симметричной, а собственные числа невырожденными. Чтобы свести задачу к нелинейной рассмотрим вектор \begin{equation} \tilde{\boldsymbol{a}} = \left( \begin{array}{c} \boldsymbol{a} \\ \lambda \end{array} \right) \end{equation} Равенство нулю $M \boldsymbol{a} - \lambda \boldsymbol{a} = 0$ нужно дополнить еще одним условием, поскольку число параметров увеличилось на 1. Пусть это будет $\boldsymbol{a}^T\boldsymbol{a}=1$, тогда \begin{equation} \boldsymbol{\varphi} = \left( \begin{array}{c} M \boldsymbol{a} - \lambda \boldsymbol{a} \\ \boldsymbol{a}^T\boldsymbol{a}-1 \end{array} \right) \end{equation} Запишем искомое решение сопряжённой задачи как $(\boldsymbol{a}_{adj}; \alpha)^T$. Тогда сопряжённое уравнение примет вид \begin{align} (M - \lambda\mathbb{I}) \boldsymbol{a}_{adj} &= \left( \frac{\partial\mathcal{F}}{\partial\boldsymbol{a}} \right)^T - 2\alpha\boldsymbol{a} \\ -\boldsymbol{a}^T \boldsymbol{a}_{adj} &= \frac{\partial\mathcal{F}}{\partial\lambda} \end{align} Поскольку матрица $(M - \lambda\mathbb{I})$, выберем $\alpha$ таким, чтобы у первого уравнения существовало решение: $\boldsymbol{a}^T (\mathcal{F}_{\boldsymbol{a}}^T - 2\alpha\boldsymbol{a}) = 0$, откуда $\alpha = (1/2)\boldsymbol{a}^T \mathcal{F}_{\boldsymbol{a}}^T$. Тогда решения получившегося уравнения \begin{equation} (M - \lambda\mathbb{I}) \boldsymbol{a}_{adj} = (1-\boldsymbol{a}\boldsymbol{a}^T) \mathcal{F}_{\boldsymbol{a}}^T = P \mathcal{F}_{\boldsymbol{a}}^T \end{equation} запишем в виде $\boldsymbol{a}_{adj} = \boldsymbol{a}_{adj0} + \gamma \boldsymbol{a}$, где $\boldsymbol{a}^T\boldsymbol{a}_{adj0} = 0$. В итоге получаем градиент \begin{equation} \left. \frac{d\mathcal{F}}{d\boldsymbol{p}} \right|_{\boldsymbol{\varphi}=0} = \mathcal{F}_{\boldsymbol{p}} - \boldsymbol{a}_{adj0}^T M_p \boldsymbol{a} + \mathcal{F}_{\lambda} \boldsymbol{a}^T A_p \boldsymbol{a} \end{equation} Последняя член уравнения выражет известный в квантовой физике результат теоремы Гельмана-Фейнмана.
Теперь, имея в распоряжении алгоритм эффективного численного расчета градиентов, рассмотрим наиболее широко распространённые градиентные методы решения оптимизационных задач.
Как и выше, рассмотрим целевую функцию $\mathcal{F}=\mathcal{F}(\boldsymbol{a}(\boldsymbol{p}))$. Итерационная формула для алгоритма наискорейшего спуска имеет вид \begin{equation} \boldsymbol{p}_{k+1} = \boldsymbol{p}_{k} + \lambda_k \nabla\mathcal{F}(\boldsymbol{p}_{k})/||\nabla\mathcal{F}(\boldsymbol{p}_{k})|| \end{equation} при этом коэффициент выбирают исходя из решения задачи одномерной оптимизации \begin{equation} \lambda_k = \mathrm{arg}\min_{\lambda} \mathcal{F}(\boldsymbol{p}_{k} + \lambda \boldsymbol{S}_k) \end{equation} где $\boldsymbol{S}_k = \nabla\mathcal{F}(\boldsymbol{p}_{k})/||\nabla\mathcal{F}(\boldsymbol{p}_{k})||$ - нормированный градиент целевой функции - направление поиска нового приближения. В такой формулировке градиент в новой точке ортогонален направлению предыдущего шага спуска. Алгоритм сходится быстро вдали от экстремума и медленно вблизи такой точки.
В отличие от алгоритма наискорейшего спуска в методе сопряженных градиентов используется информация о производных функции на предыдущих шагах. Направление поиска на текущей итерации подбирается как линейная комбинация градиента на данном шаге и направлений поиска на предыдущих шагах, а коэффициенты в комбинации так, чтобы сделать направления сопряженными относительно Гессиана.
На первом шаге $\boldsymbol{p}_{1} = \boldsymbol{p}_{0} - \lambda_0 \nabla\mathcal{F}(\boldsymbol{p}_{0}) = \boldsymbol{p}_{0} + \lambda_0 \boldsymbol{S}_{0}$, где коэффициент выбирается исходя из условия $\lambda_0 = \mathrm{arg}\min_{\lambda}(\mathcal{F}(\boldsymbol{p}_{0})+\lambda\boldsymbol{S}_{0})$. Направление спуска на втором шаге выбирается как линейная комбинация $\boldsymbol{S}_{1} = -\nabla\mathcal{F}(\boldsymbol{p}_{1}) + \beta_1 \boldsymbol{S}_{0}$ так, чтобы выполнялось указанное выше условие сопряженности $\boldsymbol{S}_{0}^T H \boldsymbol{S}_{1} = 0$. Можно показать, что это условие приводит к следующему виду коэффициента: \begin{equation} \beta_k = -\frac{||\nabla\mathcal{F}(\boldsymbol{p}_{k})||}{||\nabla\mathcal{F}(\boldsymbol{p}_{k-1})||} \end{equation} В данной формулировке метод называется методом Флетчера-Ривса. В другой распространенной формулировке Полака-Рибьера коэффициент выглядит как \begin{equation} \beta_k = -\frac{\nabla\mathcal{F}(\boldsymbol{p}_{k}) \cdot \left[ \nabla\mathcal{F}(\boldsymbol{p}_{k}) - \nabla\mathcal{F}(\boldsymbol{p}_{k-1}) \right]}{\nabla\mathcal{F}(\boldsymbol{p}_{k-1}) \cdot \boldsymbol{p}_{k-1}} \end{equation}
Методы второго порядка основаны на разложении функции вида \begin{equation} \mathcal{F}(\boldsymbol{p}+\Delta\boldsymbol{p}) \approx \mathcal{F}(\boldsymbol{p}) + \nabla \mathcal{F}(\boldsymbol{p}) \cdot \Delta\boldsymbol{p} + \frac{1}{2} \Delta\boldsymbol{p} \cdot H(\boldsymbol{p}) \Delta\boldsymbol{p} \end{equation} Одним из наиболее распространенных алгоритмов явлется метод Бройдена-Флетчера-Гольдфарба-Шанно (BFGS). Трудоёмкое вычисление Гессиана здесь заменяется на вычисление приближенного значения соответствующей обратной матрицы на каждом шаге метода $\tilde{G}_k \approx H^{-1}$, так что минимум квадратичной задачи имеет вид $\Delta\boldsymbol{p}_k = -\tilde{G}_k\nabla\mathcal{F}(\boldsymbol{p}_k)$ При этом шаг обновления приближенного обратного Гессиана \begin{equation} \tilde{G}_{k+1} = \left[ \mathbb{I} - \rho_k\Delta\boldsymbol{p}_k \Delta(\nabla\mathcal{F})_k^T \right] \tilde{G}_{k} \left[ \mathbb{I} - \rho_k \Delta(\nabla\mathcal{F})_k \Delta\boldsymbol{p}_k^T \right] + \rho_k \Delta\boldsymbol{p}_k \Delta\boldsymbol{p}_k^T \end{equation} а в качестве начального приближения обратного Гессиана можно взять единичную матрицу.
Рассмотрим задачу на собственные значения в одномерном немагнитном фотонном кристалле с периодом $\Lambda$: \begin{equation} \left[ -\frac{d}{dz^2} + \omega^2\varepsilon(z)\mu_0 \right] E_x(z) =0 \end{equation} В прямой задаче это уравнение решается на неизвестные собственные частоты и соответствующие им собственные поля. Рассмотрим, например, обратную задачу, требуется найти такую функцию $\varepsilon(x)$, при которой поле первой моды представляет собой заданную функцию $E_{x1}^{(0)}(z) = e_{x1}^{(0)}(z)\exp(ik_Bz)$. Тогда целевую функцию можно задать как \begin{equation} \mathcal{F} = \int_{\Lambda} \left| E_{x1}(z)-E_{x1}^{(0)(z)} \right|^2 dz \end{equation} Если решать уравнение на сетке $z_n = n\Delta z$, $\Delta z=\Lambda/N_z$ конечно-разностным методом, уравнение сведется к матрично-векторной задаче на собственные значения с векторами собственной функции и значений диэлектрической проницаемости в узлах сетки, а оптимизационная функция запишется через норму вектора разности \begin{equation} \mathcal{F} \sim \left| \boldsymbol{E}_{x1}(z)-\boldsymbol{E}_{x1}^{(0)(z)} \right|^2 \end{equation} Для численного решения достаточно находить минимальное собственное значение и не решать полную задачу на собственные значения.