Accelerating Adjoint Variable Method Based Photonic Optimization with Schur Complement Domain Decomposition

TL;DR

This paper introduces a novel combination of the adjoint variable method with Schur complement domain decomposition to significantly improve computational efficiency in photonic device optimization by focusing calculations on the design region.

Contribution

The paper presents a new approach that integrates Schur complement domain decomposition with the adjoint variable method for faster photonic device optimization.

Findings

Reduces computational load by focusing on the design region

Enables faster convergence in photonic optimization tasks

Improves efficiency of inverse design processes in nanophotonics

Abstract

Adjoint variable method in combination with gradient descent optimization has been widely used for the inverse design of nanophotonic devices. In many of such optimizations, the design region is only a small fraction of the total computational domain. Here we show that the adjoint variable method can be combined with the Schur complement domain decomposition method. With this combination, in each optimization step, the simulation only involves the degrees of freedom that are inside the design region. Our approach should significantly improve the computational efficiency of adjoint variable method based optimization of photonic structures.

Equations

$$\bigg{(}\frac{1}{\mu_{0}}\nabla\times\nabla\times-\omega^{2}\epsilon_{0}\epsilon_{r}(\mathbf{r})\bigg{)}\mathbf{E(r)}=-i\omega\mathbf{J(r)}$$

$$A\mathbf{E}=-i\omega \mathbf{J}$$

$$\left[\begin{array}{cc}{A}_O& A_{OB}\\ A_{BO}& A_B\end{array}\right]\left[\begin{array}{c}{\mathbf{E}}_O\\ {\mathbf{E}}_B\end{array}\right]=-i\omega \left[\begin{array}{c}{\mathbf{J}}_O\\ {\mathbf{J}}_B\end{array}\right]$$

$$A_{BO}{\mathbf{E}}_O+A_B{\mathbf{E}}_B=-i\omega {\mathbf{J}}_B$$

$$S{\mathbf{E}}_{\mathbf{O}}=-i\omega {\mathbf{J}}_S$$

$$S\equiv A_O-A_{OB}{A}_B^{-1}{A}_{BO}$$

$${\mathbf{J}}_S\equiv {\mathbf{J}}_O-A_{OB}{A}_B^{-1}{\mathbf{J}}_B.$$

$${\mathbf{E}}_B=A_B^{-1}({\mathbf{J}}_O-A_{BO}{\mathbf{E}}_O).$$

$$\frac{\partial L}{\partial \varphi }=\frac{\partial L}{\partial \mathbf{E}}\frac{\partial \mathbf{E}}{\partial \varphi }$$

$$\frac{\partial \mathbf{E}}{\partial \varphi }=-A^{-1}\frac{\partial A}{\partial \varphi }\mathbf{E}$$

$$\frac{\partial L}{\partial\phi}=-\bigg{(}\frac{\partial L}{\partial\mathbf{E}}A^{-1}\bigg{)}\frac{\partial A}{\partial\phi}\mathbf{E}=-\hat{\mathbf{E}}^{T}\frac{\partial A}{\partial\phi}\mathbf{E}$$

$$A^{T}\hat{\mathbf{E}}=\bigg{(}\frac{\partial L}{\partial\mathbf{E}}\bigg{)}^{T}$$

$$\frac{\partial L}{\partial \varphi }=-{\hat{\mathbf{E}}}^T\frac{\partial A}{\partial \varphi }\mathbf{E}$$

$$\phi_{new}=\phi_{old}+\alpha\bigg{(}\frac{\partial L}{\partial\phi}\bigg{)}$$

$$\frac{\partial L}{\partial \varphi }=-{\hat{\mathbf{E}}}_O^T\frac{\partial S}{\partial \varphi }{\mathbf{E}}_O$$

$$L\propto {\mathbf{E}}_B^{†}(\hat{M}{\mathbf{E}}_B)$$

$${\mathbf{J}}_B=\frac{\partial L}{\partial {\mathbf{E}}_B}=\hat{M}{\mathbf{E}}_B$$

$$L=\bigg{|}\int_{\Omega_{L}}{E_{z}(x)H_{x,TE_{1}}^{*}(x)+E_{z,TE_{1}}^{*}(x)H_{x}(x)dx\bigg{|}^{2}}$$

$$\frac{\partial L}{\partial\epsilon_{r}}=2k_{0}^{2}\text{Re}\bigg{(}E_{z}(x)\hat{E}_{z}(x)L^{*}\bigg{)}$$

Full text

1 Introduction

The adjoint variable method has been well established for inverse design for passive and linear devices in nanophotonics (Piggott et al., 2015; Liu et al., 2011; Lalau-Keraly et al., 2013; Veronis et al., 2004; Georgieva et al., 2002; Lu and Vučković, 2012), and has been recently generalized to active and nonlinear systems (Wang et al., 2018; Hughes et al., 2018). In particular, the technique has the computational advantage of requiring exactly two full-field simulations in order to determine the gradient of an objective with respect to an arbitrarily large number of design parameters in the domain. Combined with a gradient descent search, the method significantly decreases the number of structures that need to be probed in a given design problem.

In many optical design problems, for a given structure, one typically optimizes only part of structure (referred to as the “design region"), while the rest of the structure (referred to as the “background") is kept unchanged. (Fig. 1). For many problems, this design region can often be a very small fraction of the entire computational domain (Estakhri et al., 2019). For the standard implementation of adjoint variable method based gradient-descent optimization, in each step of the optimization, one simulates the entire structure including the background, even though the background is unchanged between the optimization steps. In order to reduce the computational cost, it would certainly be advantageous if in each optimization step one instead only performs the simulation on the design region.

In this context, the technique of a Schur complement domain decomposition is of interest. The techniques of domain decomposition have been widely used to accelerate large-scale simulations by parallelization (Dolean et al., 2015) as well as to design nested dissection algorithms for direct solvers for sparse linear systems (George, 1973; Conroy, 1990). In photonic device design, Ref. (Verweij et al., 2014) has shown that the Schur-complement domain-decomposition technique can be fruitfully applied in a discrete optimization scheme. However, in photonics, there has not been any work in applying the Schur complement domain-decomposition technique in continuous optimization schemes such as those based on the adjoint variable methods.

In this paper, we detail how the Schur complement domain decomposition approach can be combined with adjoint variable optimization in photonics device design. Since the adjoint variable method is generally applicable to any linear systems, and the Schur complement domain decomposition method describes the photonic device of interest as linear systems, it should not come as a surprise that one can combine these two methods. However, there are subtleties in the combination. In particular, in many photonic device design problem, the excitation source, as well as the field components that define the figure of merits, are not in the design region. We show that the formalism of the Schur complement naturally deals with such situations. The proposed method here, by being able to selectively simulate only the subdomain of interest, points to an avenue for increasing the computational efficiency of adjoint-based optimizations in photonics.

The rest of the paper is organized as follows. In Section 2, we outline the mathematical formalism. In Section 3, we demonstrate the computational gains of our proposed method through the design example of a mode converter. We conclude in Section 4.

2 Mathematical description of the algorithm

In this section we discuss the approach for combining the method of adjoint variable with the method of Schur complement domain decomposition. First, we briefly outline the Schur complement domain decomposition as applied to a numerically discretized solution of Maxwell’s equations in the frequency domain. The Maxwell’s equation written for the $\mathbf{E}$-field is:

[TABLE]

where $\mathbf{E(r)}$ denotes the vector electric field, $\mathbf{J(r)}$ is the source, $\epsilon_{r}(\mathbf{r})$ is the dielectric permittivity and the magnetic permeability is assumed to be $\mu_{0}$ in the entire computational domain. In Eq. (1), $\mathbf{E(r)}$ and $\mathbf{J(r)}$ are vector complex functions.

By discretizing Eq. (1) via finite differences on a Yee grid, we can convert this into a system of linear equations:

[TABLE]

where $A$ represents discretized approximation of the left hand side of Eq. (1). The $\mathbf{E}$ in Eq. (2) is now a vector contains the field components, which are the unknowns, for each spatial grid point. $\mathbf{J}$ is a vector containing information about the sources on this discretized domain.

We now consider the example domain shown in Fig. 1a), which we will also use in the numerical example in the next section. The domain consists of the design region, as denoted by subscript $O$ since this is the region where the optimization takes place, and the background region, denoted by the subscript $B$, that contains the rest of the computational domain. The background in this case contains a perfectly matched layer (light gray), a lower and upper waveguide section (light blue). The input source is in the background region. By reordering the spatial degrees of freedom in $\mathbf{E}$ and appropriately permuting the rows of $A$ in Eq. (2), we can write Eq. (2) in a $2\times 2$ block form:

[TABLE]

The submatrices $A_{O}$ and $A_{B}$ are now the discretized representations of Eq. (1) restricted to the design region and the background, respectively. The off-diagonal blocks $A_{OB}$ and $A_{BO}$ represent the spatial couplings between grid points of these two regions. The coupling occurs along their borders.

Using the second row of Eq. (3)

[TABLE]

we can solve for $\mathbf{E}_{B}$ and substitute it into the linear equation in the first row to get the equation for $\mathbf{E}_{O}$ only:

[TABLE]

where:

[TABLE]

[TABLE]

Below, we will refer to Eq. (3) as the “original” system and Eq. (6) as its corresponding “reduced” system. From Eq. (7), we see that the sources outside the design region in the original system are now represented in the reduced system by sources on the surface of the reduced computational domain, which contains only the design region. The solution to the reduced system of Eq. (5) now contains only the degrees of freedom in $\mathbf{E}_{O}$. The most expensive step in forming Eq. (4) is to compute $A_{OB}A_{B}^{-1}A_{BO}$, and $A_{OB}A_{B}^{-1}\mathbf{J}_{B}$. This can be done by a direct factorization of $A_{B}$, or by solving for $A_{B}^{-1}A_{BO}$ and $A_{B}^{-1}\mathbf{J}_{B}$ using a linear system solver. Importantly, this step, while expensive, only needs to be done once for the entire optimization process. Moreover, if the background region has significant repeated structure or is highly uniform, this operation can be accelerated by breaking it into smaller subdomains and solving $A_{OB}A_{B}^{-1}A_{BO}$ in parallel (Dolean et al., 2015).

Finally, the solution field $\mathbf{E}_{B}$ in the background region can be recovered by simple substitution of the solution $\mathbf{E}_{O}$ into Eq. (4), yielding:

[TABLE]

In many optimization problems, the figure of merit is determined based on fields that lie outside the design region. Based on Eq. (8), such figure of merit can be straightforwardly evaluated using fields and sources inside the design region only. This observation is important when we discuss the combination with the adjoint variable method below.

We now briefly review the adjoint variable method. In adjoint method, the goal is to maximize a figure of merit as described by a scalar objective function $L$ with respect to all the design variables $\phi$ (Bendsøe and Sigmund, 2004). Thus, it is of critical importance to compute the sensitivity of $L$ with respect to each element of $\phi$, (i.e. to compute $\partial L/\partial\phi$). Here for notation simplicity we have suppressed all indices for the elements of $\phi$. $\partial L/\partial\phi$ for example is a vector with a dimension equal to the number of design variables. To proceed we note that:

[TABLE]

$\frac{\partial L}{\mathbf{\partial E}}$ can be determined analytically since the objective function $L$ is typically an analytic function of $\mathbf{E}$. To determine the second term $\partial\mathbf{E}/\partial\phi$, we take derivatives of Eq. (2) with respect to $\phi$ to obtain:

[TABLE]

With a simple substitution of Eq. (10) into Eq. (9), we get:

[TABLE]

The term $\hat{\mathbf{E}}$ is the adjoint field and can be determined by solving:

[TABLE]

where the term on the right hand side is the adjoint source. The combined solutions of $\mathbf{E}$ and $\hat{\mathbf{E}}$ are enough to compute the gradient of the objective function $L$ with respect to all design variables:

[TABLE]

In the adjoint-variable-method based optimization, at each step of the optimization, we update or modify the design variables $\phi$ according the sensitivity information $\partial L/\partial\phi$:

[TABLE]

This ensures that we modify the structure in a locally optimal way for the improvement of the objective function. $\alpha$ is the learning rate. Heuristic ways of choosing or adapting $\alpha$ can be found in Ref. (Ruder, 2016).

In combining the adjoint variable method with the Schur complement domain decomposition method, the steps in Eqs. (9) - (14) remain the same, but we now apply it to Eq. (5) instead of Eq. (2). In this case, the adjoint sensitivity can be determined by computing for both $\mathbf{E}_{O}$ as well as the adjoint field $\hat{\mathbf{E}}_{O}^{T}$:

[TABLE]

In that way, simulations are done only on the degrees of freedom inside the design region (i.e. directly solving the linear system in Eq. (5). During the optimization process, from Eq. (5), the only term in $S$ that changes between each optimization step is $A_{O}$. Consequently, the term $A_{OB}A_{B}^{-1}A_{BO}$ can be precomputed once and re-used throughout the optimization process.

In many optical design problems, the objective function is defined using the field inside the background region. As an example,

[TABLE]

where $\hat{M}$ is a linear operator. In order to evaluate such cost function using the Schur complement domain decomposition approach, where only the field inside the design region is specifically simulated, we note that the adjoint source corresponding to this is:

[TABLE]

where we use the subscript $B$ to signify that the relevant field, and the source distribution, are located in the background region and are outside the design region. The corresponding source, for the reduced system, can then be determined using Eq. (7). This adjoint source in the original system thus is mapped to a source in the reduced system. The computation required for this mapping only needs to be carried out once in the beginning of the optimization process. Thus, we have shown the adjoint variable method can be readily combined with the Schur complement domain decomposition method, including for the cases where the figure of merit involve fields that are outside the design region.

3 Application to Accelerating the Design of a Mode Converter

In this section, we demonstrate the performance gains of applying our formalism in the exemplary problem of designing a mode converter on the domain shown in Fig. 1a). For simplicity, we restrict our study to two-dimensional structures (i.e. structures that are infinite and translationally invariant in the out-of-plane $z$-direction) and study the TE polarization ($E_{z}$, $H_{x}$, $H_{y}$). In this example, the design variable $\phi$ is the spatial dielectric permittivity distribution $\epsilon_{r}(x,y)$.

The objective of the device design, as shown in Fig. 1a), will be to transmit the TE0 mode of an input waveguide to the TE1-mode of the output waveguide with minimal reflections. For this purpose we will optimize for the dielectric function permittivity distribution in the design region between the two waveguides, as shown in Fig. 1a). Typical mode converters in dielectric waveguides require relatively long structures (on the order of several wavelengths) (Chen et al., 2015). Designing a compact version of this device would be very useful.

To simulate this device, the computational domain, as shown in Fig. 1a), contains a square region with an edge length of $2\lambda$ with $\lambda=1.55$ $\mu$m. The square region contains the input and output waveguides, shown as sky blue. The waveguides have a guiding layer with relative permittivity $\epsilon$ = 6.25 and width of $w=1$ $\mu$m, surrounded by a background material with a permittivity $\epsilon_{B}=2.25$. Between the input and output waveguides is the design region, which is also a square. In the optimization process, the initial structure in the design region is assumed to be uniform with a relative dielectric permittivity of 2.25. The square region is surrounded on all sides by a perfectly matched layer (PML) (Berenger, 1994) with a thickness $L_{PML}$ of 0.75 $\mu$m. The entire computational domain is discretized on a square lattice grid with a grid spacing of 50 nm.

The figure of merit $L$ for this mode converter is proportional to the overlap integral involving the electric field along a line perpendicular to the output waveguide, and the spatial profile of the TE1 magnetic field mode of the output waveguide (Lalau-Keraly et al., 2013):

[TABLE]

where $\Omega_{L}$ denotes the line along which we integrate. The derivative of this figure of merit with respective to the dielectric permittivity is $\epsilon_{r}$ is:

[TABLE]

To search for the structure which optimizes the figure of merit, we perform gradient descent optimizations for 450 iterations using both Eq. (3) and Eq. (5) and compare the total time taken to optimize the device. To solve these equations, we use an LU factorization (Davis, 2004). To solve for the adjoint field, we then use back-substitution of the LU factorization, which is considerably faster. The number of iterations is chosen to be sufficiently large so that the figures of merit converge at the end of the iteration processes. In the gradient descent optimization, we gradually reduce the learning rate $\alpha$ in Eq. (14) as the optimization progresses. We denote the total time to optimize the device using the reduced system of Eq. (5) as $T_{S}$. $T_{S}$ does not includes the overhead cost $T_{0}$ of evaluating $A_{OB}A_{B}^{-1}A_{BO}$ and $A_{OB}A_{B}^{-1}\mathbf{J}_{B}$ in Eq. (5). We denote the total time to optimize the device using the original system of Eq. (3) as $T_{A}$.

We compare the performance of solving for either original or the reduced systems, in two different optimization approaches. In the first approach, which we will refer to as the "permittivity" approach, we vary the permittivity in each pixel continuously according the gradient of the cost function with respect to the permittivity in each pixel, while restricting the permittivity of each pixel to be within the range of [2.25, 6.25]. The resulting optimal structure then can have any real value between $[2.25,6.25]$ in the design region. The second approach, which will be referred to as the "level-set" approach, involves using a level set technique so that the final structure is binary with only two permittivity values (Piggott et al., 2015; Osher and Sethian, 1988; Sethian, 1999). In the level set optimization, we first use the permittivity approach to update the permittivity continuously in each pixel until the figure of merit exceeds a threshold value. Then, we binarize the structure by using $(\epsilon+\epsilon_{B})/2$ as the threshold. The permittivity of every pixel is then changed to be either $\epsilon$ or $\epsilon_{B}$, depending on whether the existing permittivity is above or below the threshold. The binarized structure thus obtained is then further optimized using the level set technique. In general, the figure of merit achieved by the permittivity approach will be larger than that for level set methods, which is sensible given that optimizing the permittivity over a continuous range of values probes a much larger parameter space.

In Fig. 2a) and b), we compare the speed-up factor, defined as $T_{A}/(T_{0}+T_{S})$ for varying design region sizes for the two optimization approaches. The line segment for the figure-of-merit calculation is always outside the design region. The sizes are measured as the length fraction, which is the ratio of the design region side length to the side length of the total computational domain. A speed-up factor larger than one (including the overhead) indicates that solving the reduced system is more efficient than solving the original system. For both types of optimization approaches, solving the reduced system results in the substantial speed up. Even for relatively large design regions where the length fraction is almost 0.7 (in this example, this means the design region encompasses effectively the whole domain except the PML), which corresponds to almost 50% of the grid points in the domain being optimized, a speedup factor of two is still achievable. However, we observe that in this design example, a design region with a length fraction of about 0.4 is sufficient to achieve roughly the same performance in the optimized device.

In our plot, for each size of the design region we perform the gradient descent optimization only once. In practice, the number of optimization iterations may be a lot larger. For example, typically one will need to run the entire gradient descent optimization multiple times either for different initial structures, or to determine the optimal values for parameters such as the learning rate $\alpha$ in Eq. (14). In the limit of a large number of optimization iterations, the overhead cost $T_{0}$ becomes negligible. It is instead more relevant to consider the speed-up factor per epoch, defined as $T_{A}/T_{S}$. In Fig. 2a) and b) we also plot the speed-up factor per epoch as a function of design region sizes (in cyan), which gives a measure of the speed-up by our algorithm in the limit of large number of iterations.

In Fig. 2a) and b), we see that the speed-up factor decreases as the design region increases. Certainly, the matrix $S$ for the reduced system should always have a dimension that is smaller than the matrix $A$ for the original system. The matrix $S$, however, is more dense. The speed of many linear system solvers are strongly influenced by both the dimension and the sparsity of the system matrix. As a simple estimate of the size of the design region beyond which point we can reduce our expectation of a significant speed up with our approach, we can calculate the size of the design region at which the number of non-zero elements of in the matrices $S$ and $A$ are the same.

Given that our design region and our computational domain are square in geometry, let $N$ denote the number of grid points per side of the computational domain and $n$ the same for the design region. The matrix $A$ is a sparse matrix with $5N^{2}$ nonzero elements. The matrix $S$ contains two contributions to its sparsity pattern. One is from the $4n-2$ grid points on its boundary, which are densely coupled in $S$ and contribute $(4n-2)^{2}$ nonzero elements. The remaining $(n-2)^{2}$ grid points in the interior of the design region retain the same sparsity pattern as $A$ and contributes $5(n-2)^{2}$ nonzero elements to $S$, giving a total of $(4n-2)^{2}+5(n-2)^{2}$ nonzero elements. Equating the number of nonzero elements for $S$ and $A$ gives a ratio $n/N\approx 0.5$ at which the nonzero in elements in $S$ and $A$ are approximately equal (though $S$ is still a factor of 4 smaller in dimension). In our simulation, at $n/N=0.5$, the speed-up factor for both optimization approach is about $3$, which is indeed relatively low compared to smaller design regions.

In Fig. 2c), for the optimal structure, we show the modal transmission coefficient from the TE0 mode in the input waveguide to the TE1 mode in the output waveguide as a function of design region size for both the permittivity approach (blue) and the level-set approach (green). The optimal structures we obtained by solving either the original or the reduced system are the same, which is expected since mathematically these two systems are equivalent to each other. For both optimization approaches, at small size of the design region (length fraction less than 0.4), the figure of merit of the optimal structures increases as the design region size increases. However, above this length fraction, the figure of merit saturates for the permittivity optimization, which suggests that the optimal design region for this problem lies in the regime where significant speed-ups can be achieved using our formalism.

We show, in Fig. 3a), the final optimized structure at the end of 450 iterations for a design region width that is 0.37 times the width of the total domain using adjoint optimization over a continuous range of permittivities and the same in Fig. 3c) with level set methods. As expected, the optimizing over a continuous range of permittivities results in a structure where the permittivity varies between 2.25 and 6.25 in the design region, whereas the level set method results in a binary structure. In Fig. 3b), we show the steady state Re$(E_{z})$ field corresponding to the structure in Fig. 3a) and in Fig. 3d), we show the steady state Re$(E_{z})$ field corresponding to the binary structure in Fig. 3c). In both cases, the optimal structure performs the expected functionality of mode conversion.

4 Conclusion

In conclusion, we have developed a formalism based on the Schur complement to significantly reduce the computational degrees of freedom being simulated in each step of adjoint method based optimization of photonic devices. We have demonstrated the use of this algorithm by considering the optimization of a mode converter. The approach should be applicable in nanophotonic design in general. In addition, since the dimension of the reduced system can be a lot smaller than the original system, our approach may enable efficient local searches of the optimization landscape.

This work was supported by the DOE “Photonics at Thermodynamic Limits" Energy Frontier Research Center (Grant No DE-SC0019140) and by an AFOSR MURI project (Grant No. FA9550-17-1-0002).