Using Gaussian process regression for efficient parameter reconstruction

TL;DR

This paper explores the use of Gaussian process regression within Bayesian optimization to efficiently reconstruct parameters in optical scatterometry, demonstrating accelerated convergence through pre-computed simulation data.

Contribution

It introduces a method combining Gaussian process regression with Bayesian optimization and pre-computed data to improve parameter reconstruction efficiency in optical scatterometry.

Findings

Bayesian optimization outperforms local minimization algorithms.

Pre-trained Gaussian processes accelerate the optimization process.

The approach reduces the number of simulations needed for accurate reconstruction.

Abstract

Optical scatterometry is a method to measure the size and shape of periodic micro- or nanostructures on surfaces. For this purpose the geometry parameters of the structures are obtained by reproducing experimental measurement results through numerical simulations. We compare the performance of Bayesian optimization to different local minimization algorithms for this numerical optimization problem. Bayesian optimization uses Gaussian-process regression to find promising parameter values. We examine how pre-computed simulation results can be used to train the Gaussian process and to accelerate the optimization.

Equations

$$\pi (\mathbf{x})\sim \exp \left(\sum _{i=1}^r\log g_i(\mathbf{x})\right).$$

$$\nabla \times {\mu }^{-1}\nabla \times \mathbf{E}-{\omega }^2\epsilon \mathbf{E}=0,$$

$$\pi ({\mathbf{y}}_M|\mathbf{x})\sim \exp \left(-\frac{1}{2}({\mathbf{y}}_M-\mathbf{y}(\mathbf{x}){)}^T{\mathbf{\Gamma }}_l^{-1}({\mathbf{y}}_M-\mathbf{y}(\mathbf{x}))\right).$$

$$\pi (\mathbf{x}|{\mathbf{y}}_M)=\pi ({\mathbf{y}}_M|\mathbf{x})\pi (\mathbf{x})$$

$${\mathbf{x}}_{\mathrm{MAP}}=\mathrm{missing}{argmax}_{\mathbf{x}\in X}\pi (\mathbf{x}|{\mathbf{y}}_M).$$

$$F(\mathbf{x})=\frac{1}{2}({\mathbf{y}}_M-\mathbf{y}(\mathbf{x}){)}^T{\mathbf{\Gamma }}_l^{-1}({\mathbf{y}}_M-\mathbf{y}(\mathbf{x}))-\sum _{i=1}^r\log g_i(\mathbf{x}).$$

$${\mathbf{\Gamma }}_p=F^{\prime \prime }({\mathbf{x}}_{\text{MAP}}{)}^{-1}$$

$${\mathbf{x}}_{min}=\underset{\mathbf{x}\in \mathcal{X}}{\mathrm{argmin}}f(\mathbf{x}).$$

$$P({\mathbf{Y}}^{\ast })=\frac{1}{(2\pi {)}^{N/2}|\mathbf{\Sigma }{|}^{1/2}}\exp \left[-\frac{1}{2}(\mathbf{Y}-\mathbf{m}{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{Y}-\mathbf{m})\right].$$

$$\mathbf{m}$$

$$\mathbf{\Sigma }$$

$$\overline{y}({\mathbf{x}}^{\ast })$$

$$\sigma ({\mathbf{x}}^{\ast }{)}^2$$

$$\begin{array}{rl}{\alpha }_{\mathrm{EI}}({\mathbf{x}}^{\ast },y_{\min })& =\mathbb{E}[\max (0,y_{\min }-f({\mathbf{x}}^{\ast }))]\\ & =\frac{1}{2}\left[1+\mathrm{erf}\left(\frac{{\mathrm{y}}_{\min }-\overline{\mathrm{y}}({\mathbf{x}}^{\ast })}{\sqrt{2}\sigma ({\mathbf{x}}^{\ast })}\right)\right](y_{\min }-\overline{y}({\mathbf{x}}^{\ast }))+\frac{\sigma ({\mathbf{x}}^{\ast })}{\sqrt{2\pi }}\exp \left(-\frac{(y_{\min }-\overline{y}({\mathbf{x}}^{\ast }){)}^2}{2\sigma ({\mathbf{x}}^{\ast }{)}^2}\right).\end{array}$$

$$k(\mathbf{x},{\mathbf{x}}^{\prime })=s^2\left(1+\sqrt{5}r(\mathbf{x},{\mathbf{x}}^{\prime })+\frac{5}{3}r(\mathbf{x},{\mathbf{x}}^{\prime }{)}^2\right)\exp \left(-\sqrt{5}r(\mathbf{x},{\mathbf{x}}^{\prime })\right)\text{with}r(\mathbf{x},{\mathbf{x}}^{\prime }{)}^2=\sum _i\frac{({\mathbf{x}}_i-{\mathbf{x}}_i^{\prime }{)}^2}{l_i^2}.$$

$$d=\underset{i=1,\cdots ,6}{\max }\frac{\left|(\mathbf{x}-{\mathbf{x}}_{\mathrm{MAP}}{)}_i\right|}{{\sigma }_i}$$

Full text

Copyright 2019 Society of Photo‑Optical Instrumentation Engineers (SPIE). One print or electronic copy may be made for personal use only. Systematic reproduction and distribution, duplication of any material in this publication for a fee or for commercial purposes, and modification of the contents of the publication are prohibited.

1 Introduction

Geometry reconstruction based on scatterometric data is a challenging numerical task. The sample structures and the measurement process are typically described using many parameters. This leads to high-dimensional optimization problems of finding shape parameters that are in agreement with the experimental data [1, 2, 3].

The forward-problem of computing the scattering behavior of the setup at a given point in the parameter space requires to rigorously solve Maxwell’s equations. We use our finite-element method (FEM) implementation JCMsuite[4, 5] to this aim. The optimization process can be often solved more efficiently by providing parameter derivatives. Therefore, we compute the gradient of the solution of Maxwell’s equations by automatic differentiation [1].

A large number of minimization algorithms can be used to solve the inverse problem of reconstructing the shape parameters. In the considered case the objective function has only a small number of local minima such that local minimization algorithms should be very efficient. Examples for gradient-based local optimization methods are the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm and its low-memory, bound-constrained extension L-BFGS-B  [6, and references therein] as well as the truncated Newton method [7]. An example for a gradient-free method is the downhill simplex algorithm (also known as Nelder-Mead method) [8].

Recently, techniques from the field of machine-learning have been employed for the optimization of photonic nanostructures. For example, deep neural networks trained with thousands of simulation results have been employed as accurate models for mapping a geometry to an optical response and vice versa almost instantaneously [9]. In this work, we consider consider Gaussian processes as a method to learn the behavior of the objective function [10]. A popular method that employs Gaussian processes is Bayesian optimization [11]. Bayesian optimization is regularly used in machine learning applications [11, 12, 13]. In the field of nano-optics it has been, e.g., employed to optimize ring resonator-based optical filters [14] and chiral scatterers [15].

Bayesian optimization derives promising parameter values by means of Bayesian inference based on all previous function evaluations. This is in contrast to local optimization methods, which only use few of the previous data points to determine new parameters. This statistical inference can often reduce the number of required simulations [16].

In the context of a parameter reconstruction it is possible to compute the system response for many parameter values in advance. Provided with scatterometry data of a specific structure, the deviation between numerical and experimental response for the pre-computed parameter values can serve as training data for the Gaussian process. We investigate to which extent this training can speed up the parameter reconstruction.

2 Scatterometric measurement configuration

In order to assess the optimization methods for solving the inverse problem, we use a critical dimension metrology setup studied already in previous publications [1, 17]. The scatterometric measurement was executed at Physikalisch-Technische Bundesanstalt (PTB) and a detailed investigation of the data and optical model can be found in a recently published paper [2].

For the sake of completeness, we briefly review the experimental setup and the optical model: A silicon grating (1D periodic lines) with nominal pitch of $p_{x}=50$ nm and nominal line-width of $CD=25$ nm was used as scattering target in a goniometric setup with an inspection wavelength of $\lambda=266\,$nm. A light beam with defined polarization and angle of incidence (inclination angle $\theta$, rotation angle $\phi$) illuminates the target. Due to the grating period and the wavelength only the zeroth diffraction order is present and the intensity of the reflected light in this diffraction order is recorded for S- and for P-polarized illumination at different inclination angles $\theta$. Two azimuthal rotations ($\phi=0$ and $\phi=90$°) are recorded.

A schematic of the measurement is shown in Figure 1 (left). The measured data set used in this study is plotted in Figure 2 (circles). See Ref. [1] for further explanations.

Optical model.

The aim of the optical model is to describe the measurements as a function of the microstructure’s parameters $\mathbf{x}\in\mathbb{R}^{n}$. The shape of the silicon line is parameterized with $n=6$ free parameters: the line height $h$, the line width (critical dimension) at $h/2$, $CD$, the oxide layer thickness, $t$, the sidewall angle, $SW\!\!A$, and the top and bottom corner roundings, $r_{\textrm{top}}$ and $r_{\textrm{bot}}$. Their definitions can be found in Figure 1 (center).

In the reconstruction we allow for parameter values within large intervals describing a wide range of line shapes. To avoid non-physical self-intersections, we demand these to stay in an admissible bounded region. The admissible region can in general be defined as $A=\{\mathbf{x}\in\mathbb{R}^{m}\mid g_{i}(\mathbf{x})\geq 0,i=1,\dots,r\}$ with smooth scalar valued functions $g_{i}$. For example, we demand that the corner rounding radius at the top of the line is smaller than half the width at the top. The admissible region is included into the objective function by means of a prior density as

[TABLE]

The scattering of monochromatic light off the nanoscopic line grating is described by the linear Maxwell equations in frequency domain. These lead to a single second order partial differential equation

[TABLE]

where $\mathbf{\varepsilon}$ and $\mathbf{\mu}$ are the permittivity and permeability tensors, and $\omega$ is the time-harmonic frequency. We employ the finite-element (FEM) electromagnetic field solver JCMsuite [4, 5], which has been successfully used in scatterometric investigations ranging from the optical [18] to the EUV and X-ray regimes [19, 20] on 2D (e.g., line masks) and 3D (e.g., FinFETs, contact holes) scattering targets.

Inverse problem.

To reconstruct likely parameters, we use the same objective function as in [1]. That is, we minimize the conditional probability $\pi(\mathbf{x}|\mathbf{y}_{M})$ of the parameter vector $\mathbf{x}$ given the measurement vector $\mathbf{y}_{M}$, also known as posterior probability. The posterior is given by Bayes’ theorem as the product of likelihood $\pi(\mathbf{y}_{M}|\mathbf{x})$ and prior probability $\pi(\mathbf{x})$. We assume independently normally distributed measurement errors with zero mean and diagonal covariance $\mathbf{\Gamma}_{l}$:

[TABLE]

The posterior density

[TABLE]

is maximized in order to find the most probable parameter values

[TABLE]

Equivalently, $\mathbf{x}_{\rm MAP}$ can be found by minimizing the logarithm of $\pi(\mathbf{x}|\mathbf{y}_{M})$,

[TABLE]

As the mapping $\mathbf{y}(\mathbf{x})$ from model parameters to scatterometry measurements is nonlinear, the posterior density $\pi(\mathbf{x}|\mathbf{y}_{M})$ is non-normal and can exhibit multiple local maxima. In the previous numerical study only two local maxima were found such that most of the runs of a gradient descent method converged efficiently to the global maximum $\mathbf{x}_{\textrm{MAP}}$ located at $CD=25.38$  nm, $h=48.08$ nm, $SW\!\!A=86.98$°, $t=4.94$ nm, $r_{\textrm{top}}=10.37$ nm, and $r_{\textrm{bot}}=4.79$ nm. The local uncertainties were quantified in terms of the covariance matrix

[TABLE]

yielding the standard deviations $\sigma_{CD}=0.395$  nm, $\sigma_{h}=2.484$ nm, $\sigma_{SW\!\!A}=0.999$°, $\sigma_{t}=0.162$ nm, $\sigma_{r_{\textrm{top}}}=4.289$ nm and $\sigma_{r_{\textrm{bot}}}=3.217$ nm [1].

In Figure 2 the experimental and simulated intensities are shown as function of inclination angle $\theta$. The four different angular spectra refer to the different polarizations and azimuthal orientations of the illumination. We observe an almost perfect alignment of the simulated data for $\mathbf{x}=\mathbf{x}_{\textrm{MAP}}$ and the PTB measurements.

3 Bayesian optimization method

The goal of every optimization algorithm is to identify the minimum (or maximum) of an unknown objective function $f$ in a certain design space $\mathcal{X}\subset\mathbb{R}^{d}$,

[TABLE]

The basic idea of Bayesian optimization is to treat the unknown objective as a random function, i.e. a stochastic model on a continuous domain $\mathcal{X}$. Based on the previous observation of the objective the algorithm identifies parameter values where it is expected to find a smaller function value [11, 21].

Gaussian processes.

Gaussian processes (GP) are frequently used as the stochastic model in Bayesian optimization. A stochastic process $(X_{x})_{x\in\mathcal{X}}$ is a Gaussian process if for any $N$ points $\mathbf{x}_{1}^{*},\cdots\mathbf{x}_{N}^{*}\in\mathcal{X}$ the probability of the objective to be equal to $\mathbf{Y}=(y_{1},\cdots,y_{N})$ follows a multivariate Gaussian distribution

[TABLE]

with a mean vector $\mathbf{m}$ and a covariance matrix $\mathbf{\Sigma}$ [22].

The Gaussian process is defined by a covariance function (or kernel) $k:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ and a mean function $\mu:\mathcal{X}\rightarrow\mathbb{R}$. Without prior information the mean vector $\mathbf{m}$ and a covariance matrix $\mathbf{\Sigma}$ evaluate to $\mathbf{m}=[\mu(\mathbf{x}_{1}^{*}),\cdots,\mu(\mathbf{x}_{N}^{*})]^{T}$ and $(\mathbf{\Sigma})_{ij}=k(\mathbf{x}_{i}^{*},\mathbf{x}_{j}^{*})$. After $M$ iterations, the function values $Y=(y_{1},\cdots,y_{M})=(f(x_{1}),\cdots,f(x_{M}))$ are known. Using Bayesian inference, the mean vector and the covariance matrix then reads

[TABLE]

with $\mathbf{m}_{1}=[\mu(\mathbf{x}_{1}),\cdots,\mu(\mathbf{x}_{M})]^{T}$, $\mathbf{m}_{2}=[\mu(\mathbf{x}_{1}^{*}),\cdots,\mu(\mathbf{x}_{N}^{*})]^{T}$, $(\mathbf{\Sigma}_{11})_{ij}=k(\mathbf{x}_{i},\mathbf{x}_{j})$, $(\mathbf{\Sigma}_{22})_{ij}=k(\mathbf{x}_{i}^{*},\mathbf{x}_{j}^{*})$, $(\mathbf{\Sigma}_{21})_{ij}=k(\mathbf{x}_{i}^{*},\mathbf{x}_{j})$, and $(\mathbf{\Sigma}_{12})_{ij}=k(\mathbf{x}_{i},\mathbf{x}_{j}^{*})$.

Gaussian processes can be easily extended to incorporate not only data on the objective function but also of its derivatives [23]. For example, a gradient observation $\nabla f(\mathbf{x}_{k})$ leads to additional entries in the covariance matrix formed by derivatives of the covariance function $\nabla_{\mathbf{x}_{k}}k(\mathbf{x}_{i},\mathbf{x}_{k})$ and $\nabla_{\mathbf{x}_{i}}\nabla_{\mathbf{x}_{k}}k(\mathbf{x}_{i},\mathbf{x}_{k})$.

From Eqs. (9–11) follows that at any position $\mathbf{x}^{*}\in\mathcal{X}$ the unknown function value is normally distributed, i.e. $f(\mathbf{x}^{*})\sim\mathcal{N}(\overline{y},\sigma^{2})$ with

[TABLE]

Based on the normal distribution the next parameters values are chosen to maximize a specific acquisition function. A popular choice is the expected improvement, i.e. the expectation value of the improvement $\max(0,y_{\rm min}-f(\mathbf{x}^{*}))$ with respect to the currently known minimal function value $y_{\rm min}$

[TABLE]

Other common acquisition functions are the probability of improvement or the lower confidence bound [11].

Hyperparameter choice.

Using different covariance functions $k(\mathbf{x},\mathbf{x}^{\prime})$, Gaussian processes allow to approximate a large class of random functions. A popular choice is the Matérn $5/2$ covariance function

[TABLE]

Moreover, we choose a constant mean function $\mu(\mathbf{x})=\mu_{0}$. The values of the hyperparameters $\mathbf{w}=(\mu 0,s^{2},$ $l_{1}^{2},\cdots,l_{d}^{2})$ are essential for the performance of the optimization procedure. The idea is to maximize the likelihood $P(\mathbf{Y})=P_{\mathbf{w}}(\mathbf{Y})$ of all known objective function values with respect to the hyperparameters. This hyper parameter optimization is computationally expensive and is only performed if the derivatives of $P_{\mathbf{w}}(\mathbf{Y})$ with respect to the length scales exceed a certain threshold [24].

Learning from offline calculations.

In a typical parameter reconstruction setup many specimens with the same type of geometry have to be probed. Therefore, we examine whether the reconstruction process can be accelerated by pre-comuting the optical model $\mathbf{y}(\mathbf{x})$ for many parameter values $\mathbf{x}\in X$. That is, provided with a measurement result $\mathbf{y}_{M}$ the parameters $X$ and the corresponding posterior probabilities $P_{0}=\{\pi(\mathbf{x}|\mathbf{y}_{M})\;|\;\mathbf{x}\in X\}$ and the parameter derivatives $P_{i}=\{\partial\pi(\mathbf{x}|\mathbf{y}_{M})/\partial\mathbf{x}_{i}\;|\;\mathbf{x}\in X\}$ for $i=1,\cdots,6$ are used to initialize the Gaussian process underlying Bayesian optimization. We draw $X$ from a pseudo-random Sobol sequence and use only parameters that meet the constraints, i.e. $X\subset A$ [25, 26].

4 Results

In order to assess the performance of the different optimization approaches we have conducted six independent optimization runs of 150 iterations with different initial conditions for each method. That is, the Downhill-Simplex method, L-BFGS-B, Newton Conjugate-Gradient, and Bayesian optimization were started from six different random initial points. Whenever a local minimization method converged to a local minimum, it was restarted at a different position. Moreover, six independent sets of points $(X^{(i)})_{i=1,\cdots,6}$ of 100 training samples and corresponding posterior probabilities and derivatives $(P_{0}^{(i)},\cdots,P_{6}^{(i)})_{i=1,\cdots,6}$ were prepared. The sets were used to initialize the Gaussian processes (”Bayesian optimization + training”). Apart from the downhill simplex methods, all methods make use of derivative information to determine the next sampling point.

Since the error function $F(\mathbf{x})$ [see Eq. (6)] varies between several orders of magnitude, we minimize its decadic logarithm $\lg[F(\mathbf{x})]$. Figure 3 compares the performance of the optimization approaches for minimizing $\lg[F(\mathbf{x})]$. On the left, the average objective value is shown as a function of the number of simulations for each of the different optimization methods. Because some optimization runs can fail by being trapped into a local minimum, the average objective value is not always meaningful. Therefore, the right plot shows also the median number of simulations needed to obtain objective values below a certain threshold. The median is less sensitive against failed runs, as long as at least four of the 6 runs are successful. Bayesian optimization needs on average significantly less simulations to find the minimum than the local optimization methods. Surprisingly, the performance of the derivative-free downhill-simplex method is similar to the gradient descent methods in minimizing $\lg[F(\mathbf{x})]$. Only close to the global minimum L-BFGS-B converges better to objective values below 1.1. Clearly, the 100 training samples lead to a significant speed-up of Bayesian optimization. Only after about 20 iterations the non-trained approach converges faster to the global minimum. We attribute this to the following behavior: Both Bayesian optimization methods converge into a region close to the global minimum. Because of small parameter derivatives, the expected improvement in this region becomes very small. Provided with many training samples the second approach now identifies other regions in the parameter space where some improvement can be expected and tends to explore these regions. If no training samples are available, the same happens at a later stage.

Figure 4 shows the maximum distance to the global minimum as a function of the number of simulations. The distance is measured in units of the measurement uncertainty for all six geometry parameters, i.e.,

[TABLE]

Bayesian optimization with training converges after a median number of 9 iterations to a region within the measurement uncertainty ($d=1$), i.e. where the derivatives of $F(\mathbf{x})$ become small. Without training, Bayesian optimization needs more than 25 iterations to converge to the same accuracy level. From the local minimization methods, only L-BFGS-B converges to the measurement uncertainty within 150 simulations.

5 Conclusion

We have compared the performance of Bayesian optimization to different local minimization algorithms for a specific example of a geometry parameter reconstruction based on scatterometry data. We find that Bayesian optimization finds the most probable parameters with a significantly smaller number of simulation results. A training with pre-computed simulations can further speed up the reconstruction up to the point where parameter values within the measurement uncertainty are identified. Thereafter, the convergence speed decreases since the expected improvement around other training samples is larger. Other acquisition functions such as the probability of improvement could prevent the exploration of other parts of the parameter space. In order to improve the convergence speed it might be beneficial to change to a different Bayesian optimization strategy or to a gradient descent method when the expected improvement becomes small.

Acknowledgments

This work is partially funded through the project 17FUN01 ”BeCOMe” within the Programme EMPIR. The EMPIR initiative is co-founded by the European Union’s Horizon 2020 research and innovation program and the EMPIR Participating Countries. Further, we acknowledge support by the Central Innovation Programme for SMEs of the German Federal Ministry for Economic Affairs and Energy on the basis of a decision by the German Bundestag (ZF4450901).