# Effect of Choice of Boundary Condition on the Numerical Efficiency of Direct
Solutions of Finite Difference frequency Domain Systems with Perfectly Matched
Layers

Nathan Zhao
Department of Applied Physics
Stanford University
Stanford, CA 94305
<EMAIL_ADDRESS>
&Shanhui Fan
Department of Electrical Engineering
Stanford University
Stanford, CA 94305
<EMAIL_ADDRESS>

###### Abstract

Direct solvers are a common method for solving finite difference frequency
domain (FDFD) systems that arise in numerical solutions of Maxwell’s
equations. In a direct solver, one factorizes the system matrix. Since the
system matrix is typically very sparse, the fill-in of these factors is the
single most important computational consideration in terms of time complexity
and memory considerations. As a result, it is of great interest to determine
ways in which this fill-in can be systematically reduced. In this paper, we
show that in the context of commonly used perfectly matched boundary layer
methods, the choice of boundary condition behind the perfectly matched
boundary layer can be exploited to reduce fill-in incurred during the
factorization, leading to significant gains of up to 40% in the efficiency of
the factorization procedure. We illustrate our findings by solving linear
systems and eigenvalue problems associated with the FDFD method.

##  1 Introduction

In any finite difference simulation, the domain must be truncated and a
boundary condition is imposed. Common boundaries include the periodic boundary
condition or the Dirichlet boundary condition. However, in many applications,
one wants to simulate an "infinite" domain in the sense that any wave incident
on the boundary should have no reflection. In one dimension, it is possible to
eliminate all reflections via the Mur absorbing boundary condition [1].
However in $n>2$ dimensions, there is no absorbing boundary condition that can
eliminate reflections for waves incident from all angles. Instead, it is
common to surround the computational domain by perfectly matched layers (PML),
a specially designed anisotropic absorber [1, 2]. The PML region is typically
truncated using a boundary condition such as the periodic or the Dirichlet
boundary condition.

The PML has been extensively used in computational electromagnetics, for both
frequency and time-domain solutions of the Maxwell’s equations [1, 3, 4, 5,
2]. In the time domain, significant efforts have been invested in reducing the
reflection for a given number of layers in the PML region [6, 7, 8]. Various
formalisms of the PML have also been developed to ensure strong suppression of
reflection over broad bandwidth [9, 2]. In the frequency domain, there are
additional considerations in the applications of PML. For example, Ref. 10 has
shown that the condition number of the system matrix, and hence the
performance of iterative solutions of Maxwell’s equations in the frequency
domain, can vary drastically for different formalisms of PML, even though
these formalisms offer similar wall-clock performance in time-domain
algorithms.

In this paper we provide a discussion of PML in the context of direct solver
for Maxwell’s equations in the frequency domain. In such direct solvers, one
typically performs an LU decomposition of the system matrix. We show that the
fill-in in the LU decomposition strongly depends on the boundary condition
that is used to truncate the PML. In particular, the use of the Dirichlet
boundary condition leads to far less fill-in and hence significant efficiency
gain, as compared with the use of the periodic boundary condition. This result
is in contrast with both time-domain algorithms, as well as with iterative
algorithms in the frequency domain for the solutions of Maxwell’s equations,
where the use of these two boundary conditions do not significantly affect the
performance of the algorithms. Our result is important for the development of
direct solvers for Maxwell’s equations in the frequency domain.

##  2 A Brief Review of the Finite Difference Frequency Domain Method

###  2.1 Linear system formulation

The Maxwell’s equation written for the $\mathbf{H}$-field is:

$\nabla\times\frac{1}{\overleftrightarrow{\epsilon_{r}}(\mathbf{r})}\nabla\times\mathbf{H(r)}-\omega^{2}\mu_{0}\mathbf{H(r)}=-i\omega\mathbf{M(r)}$
(1)

where $\overleftrightarrow{\epsilon_{r}}$ denotes a three-dimensional
dielectric tensor, $\omega$ is the frequency, $\mathbf{H(r)}$ denotes the
vector magnetic field, $\mathbf{M(r)}$ is the vector magnetic source, and the
magnetic permeability is assumed to be $\mu_{0}$ in the entire computational
domain. In the finite-difference frequency-domain (FDFD) method, one
discretizes Eq. (1) via finite differences on a Yee grid [11, 12, 2, 13] to
obtain a system of linear equations:

$\hat{A}\mathbf{x}=\mathbf{b}$ (2)

The matrix $\hat{A}$ corresponds to the left hand side of Eq. (1) and
$\mathbf{x}$ is a vector corresponding to the discretized version of
$\mathbf{H}$. $\mathbf{b}$ corresponds to the source term on the right hand
side.

In order to obtain a solution to Eq. (2), a number of algorithms exist, such
as Gaussian elimination [14]. However, a common method is to factorize the
matrix $\hat{A}$ into $\hat{A}=\hat{L}\hat{U}$, where $\hat{L}$ is a lower
triangular matrix and $\hat{U}$ is an upper triangular matrix. The solution
can then be obtained:

$\mathbf{x}=\hat{U}^{-1}(\hat{L}^{-1}\mathbf{b})$ (3)

Eq. (3) is typically called back-substitution. The computation of $L^{-1}b$
and $U^{-1}(L^{-1}b)$ is cheaper as compared with the full Gaussian
elimination of the system matrix $\hat{A}$ [15, 16]. However, the process of
determining the factorization requires a modified version of Gaussian
elimination, so the overall time complexity of using the LU method is the same
as the Gaussian elimination [14]. On the other hand, if one has multiple right
hand sides to solve, then the LU factorization becomes efficient as the $L$
and $U$ factors can be re-used.

###  2.2 Linear Eigenvalue System Formulation

The LU decomposition is also useful for solving eigenvalue problems in
electromagnetics, an example of which is:

$\nabla\times\frac{1}{\overleftrightarrow{\epsilon_{r}}(\mathbf{r})}\nabla\times\mathbf{H(r)}=\omega^{2}\mu_{0}\mathbf{H(r)}$
(4)

Following the same discretization procedure as that leading to Eq. (2), we
obtain the eigenvalue problem:

$\hat{A}\mathbf{x}=\lambda\mathbf{x}$ (5)

where $\hat{A}$ now corresponds to the left-hand side of Eq. (5) and
$\lambda=\mu_{0}\omega^{2}$. Eigenvalues are typically determined iteratively
by using the Arnoldi iteration and its variants where one repeatedly applies
the matrix $\hat{A}$ to a test vector [17, 18]. However, direct application of
these methods can only produce the largest $k$ eigenvalues of a system. In
computational electromagnetics, one typically is interested instead in
eigenvalues that are near a target value $\sigma$. For such a problem, a
reformulation of the eigenvalue problem is required [17, 19]. Instead of
solving Eq. (5) we solve:

$(\hat{A}-\sigma\hat{I})^{-1}\mathbf{x}=\beta\mathbf{x}$ (6)

where $\hat{I}$ represents the identity matrix. The eigenvalues $\beta$ can be
related to eigenvalues of $A$ via: $\beta=1/(\lambda-\sigma)$. Thus the
eigenvalues of $\hat{A}$ that are closest to $\sigma$ correspond to the
largest eigenvalues of $(\hat{A}-\sigma I)^{-1}$. In applying iterative
methods to obtain the largest eigenvalues of $(\hat{A}-\sigma I)^{-1}$, we
need to repeatedly apply $(\hat{A}-\sigma\hat{I})^{-1}$. In the FDFD method,
$\hat{A}$ is sparse, but $(\hat{A}-\sigma\hat{I})^{-1}$ is fully dense. For
large problems it is too costly memory-wise to form
$(\hat{A}-\sigma\hat{I})^{-1}$ explicitly. The efficient way to do this is to
compute an LU factorization of $\hat{A}-\sigma\hat{I}$ [20]. Thus, the LU
factorization is important for the eigenvalue problems that arise in the FDFD
method as well.

###  2.3 Brief Review of Boundary Conditions

In the FDFD method the computational domain is finite and a boundary condition
needs to be applied on the surface of the computational domain. This is
illustrated in Fig. 1, where we show an undirected graph representation of a
simple square grid consisting of $N=N_{x}\times N_{y}$ nodes with
$N_{x}=N_{y}=10$. The points at which the fields are evaluated are termed
nodes. We denote the fields evaluated at the nodes as $\mathbf{x}(i,j)$, where
$i$ is the row index and $j$ is the column index. For a 2D problem, the
Maxwell’s equations can be written in terms of only $E_{z}$ or $H_{z}$ fields,
where $z$ is the out-of-plane direction, hence $\mathbf{x}$ can correspond
either to the $E_{z}$ or the $H_{z}$ fields. The outermost nodes are the
boundary nodes where the boundary condition is enforced. The edges represent
the node-to-node coupling as derived from finite-difference approximation of
Eq. (1). In Fig. 1a, we show the periodic boundary, which can be described as:
$\mathbf{x}(N_{x},j)=\mathbf{x}(0,j)$ and
$\mathbf{x}(i,N_{y})=\mathbf{x}(i,0)$. In Fig. 1b, we show a Dirichlet
boundary condition:
$\mathbf{x}(N_{x},j)=\mathbf{x}(0,j)=\mathbf{x}(i,N_{y})=\mathbf{x}(i,0)=0$.
When $\mathbf{x}$ corresponds to the $\mathbf{H}$($\mathbf{E}$) fields, the
Dirichlet boundary condition corresponds to the perfectly magnetic (electric)
conductor boundary condition.

Figure 1: Two dimensional finite difference grid with a five point stencil. a)
A grid truncated with the periodic boundary condition. The boundary
connections which connect left boundary to right and top to bottom are shown
as trailing black edges. b) Same grid with the Dirichlet boundary condition.

In typical FDFD simulations one often surrounds a computation domain of
interest with perfectly match layers (PML) [1, 10]. These layers are used to
ensure that any wave exiting the computational domain is nearly perfectly
absorbed, which is important for simulating an open electromagnetic system.
When the PML are used, the entire computational domain consists of the domain
of interest and the PML regions. Outside the PML regions, the computational
domain still needs to be truncated. Either the periodic or the Dirichlet
boundary conditions can be used for this purpose.

The choice of either the periodic or the Dirichlet boundary conditions outside
the PML regions does not affect their absorption properties. Moreover, in the
finite difference time domain method, or when iterative solvers are used in
the finite difference frequency domain method, these two boundary conditions
result in very similar computational cost. Therefore, there has been very
little attention paid to the choice between these two boundary conditions. As
we will show in this paper, however, this choice in fact makes a significant
difference when direct solvers are used in the finite difference frequency
domain method, since the choice significantly influences the fill-in behaviors
in the LU decomposition of the system matrix.

##  3 A Brief Review of the Fill-in of Sparse Linear Systems

In general, $\hat{A}$ for Eq. (2) or Eq. (5) is sparse. A key consideration in
the LU factorization is the issue of fill-in, which is the accumulation of
additional non-zero elements in $\hat{L}$ and $\hat{U}$ outside those of
$\hat{A}$. In general, fill-in during the factorization process for a sparse
linear system depends sensitively on the detailed sparsity pattern of the
system matrix. Consider, for example, the following factorization of a matrix:

$\begin{bmatrix}p&p&p&p&p\\\ p&p&&&\\\ p&&p&&\\\ p&&&p&\\\ p&&&&p\\\
\end{bmatrix}=\begin{bmatrix}p&&&&\\\ p&p&&&\\\ p&p&p&&\\\ p&p&p&p&\\\
p&p&p&p&p\\\ \end{bmatrix}\begin{bmatrix}1&p&p&p&p\\\ &1&p&p&p\\\ &&1&p&p\\\
&&&1&p\\\ &&&&1\\\ \end{bmatrix}$ (7)

where $p$ denotes a matrix element which is nonzero (element values may be
different). Because the first row has no zeros, the factorization process
results in $\hat{L}$ and $\hat{U}$ factors that are dense [14]. In contrast,
consider the factorization of an equivalent system as obtained from the left
hand size of Eq. (7) by permuting the row index:

$\begin{bmatrix}p&&&&p\\\ &p&&&p\\\ &&p&&p\\\ &&&p&p\\\ p&p&p&p&p\\\
\end{bmatrix}=\begin{bmatrix}p&&&&\\\ &p&&&\\\ &&p&&\\\ &&&p&\\\ p&p&p&p&p\\\
\end{bmatrix}\begin{bmatrix}1&&&&p\\\ &1&&&p\\\ &&1&&p\\\ &&&1&p\\\ &&&&1\\\
\end{bmatrix}$ (8)

The $\hat{L}$ and $\hat{U}$ factors resulting from Gaussian elimination do not
have any fill-in. In general, the order of the equations as well as the order
of the variables can drastically change the fill-in. The operations which
generate these reorderings can be expressed as preconditioners on $\hat{A}$
via a set of permutation matrices $\hat{P}$ and $\hat{Q}$:
$\hat{A}\rightarrow\hat{P}\hat{A}\hat{Q}$, where $\hat{P}$ performs row-wise
and $\hat{Q}$ performs column-wise permutations. Unfortunately, the problem of
finding the optimal permutation to produce the least amount of fill-in is NP-
complete [21], which means there are no known, at-worst polynomial time
algorithms to determine the optimal permutation. Instead, there is a wide
array of different heuristic or approximate methods to reduce fill-in.

The state-of-the-art re-orderings are based on the exact minimum degree
algorithm and their variants, such as the methods of approximate minimum
degree (AMD) and nested dissection [22, 23, 24, 25, 26]. All of these methods
effectively rely on the minimum degree heuristic, which is the observation
that first eliminating variables in the factorization process with fewer
couplings (or nonzeros in the row/column) tends to yield less fill-in overall.

##  4 Reducing Fill-in by Reducing Connectivity within the PML

As discussed in Section 2.3, the choice of boundary condition backing a PML
can vary. As a result, we have the freedom to choose the boundary condition
behind the PML. In this section, we now show how the choice of boundary
condition behind the PML can affect the efficiency of sparse direct solvers,
specifically impacting the fill-in during factorization.

From Figure 1, we can apply the minimum degree heuristic by counting the
neighbors of the individual nodes on the grid. The number of neighbors
corresponds to the number of nonzero elements in each row of $\hat{A}$. For
the periodic boundary case shown in Fig. 1a, all nodes including the boundary
have 4 neighbors. By the minimum degree heuristic, there is no preference in
picking which rows to eliminate first. Meanwhile, for the Dirichlet boundary
in Fig. 1b, the four corner nodes have 2 neighbors, the edge nodes have 3
neighbors and then all other nodes have 4 neighbors. By the minimum degree
heuristic in this case, we can see that the rows corresponding to the boundary
nodes should be eliminated first.

Finally, we can take advantage of the fact that the nodes adjacent to the
Dirichlet boundary have fields that are generally close to zero, since the
fields have been strongly by the PML. As a result, we can also further modify
the connectivity to minimize the number of neighbors per node, as we show in
Fig. 2:

Figure 2: a) Optimized grid with reduced connectivity of the boundary. b) Same as a) with the nodes color coded by their degree (deg) in the graph representation. Boundary Condition of Grid | Total Grid Couplings for ($N_{x}=N_{y}$) | $N_{x}=100$
---|---|---
Periodic | $4N_{x}^{2}$ | 40000
Dirichlet | $4(N_{x}-2)^{2}+4(N_{x}-2)+8$ | 38816
Modified Dirichlet | $4(N_{x}-2)^{2}+N_{x}-2$ | 38514
Table 1: Changes in the number of couplings in a grid of $N$ nodes for
different boundary conditions. The 2nd column shows the analytic expression
for the number of total couplings in the grid. The third column shows the
values of these expressions for the specific case of $N=10000$.

In Fig. 3, we show the influence of the boundary conditions on the fill-in
behavior for the small 10$\times$10 grids shown in Figs. 1 and 2. The use of
such a small grid allows easier visualization of various quantities involved.
We solve a 2D problem corresponding to setting all terms involving
$\partial_{z}$ to zero in Eq. (1). We simulate a magnetic source in vacuum for
the magnetic field perpendicular $H_{z}$ to the 2D plane:

$\bigg{(}\frac{\partial}{\partial
x}\frac{1}{\epsilon_{r}(x,y)}\frac{\partial}{\partial
x}+\frac{\partial}{\partial
y}\frac{1}{\epsilon_{r}(x,y)}\frac{\partial}{\partial
y}+\mu_{0}\omega^{2}\bigg{)}H_{z}=i\omega M_{z}$ (9)

where we assume $\epsilon_{r}$ is isotropic, $M_{z}$ is the $z$-component of
the magnetic source current. We consider the periodic boundary conditions, the
Dirichlet boundary conditions, and the modified Dirichlet boundary conditions
as discussed above where the connectivity at the boundary is reduced. The
first column of Figure 3 shows the sparsity pattern of the system matrices for
the three boundary conditions. The sparsity patterns of the three matrices are
similar since they differ only on the boundary. Compared with the system
matrix from the Dirichlet boundary condition (Fig. 3d), the matrix from the
periodic boundary condition (Fig. 3a) has extra non-zero elements at the
lower-left and upper-half corners in the figure due to the long-range coupling
induced by the periodic boundary condition. The matrix from the modified
Dirichlet (Fig. 3g) boundary condition has less non-zero elements near the
diagonal as results from the modifications.

The second column in Figure 3 shows the $\hat{L}$ factor as obtained by
directly factoring the corresponding matrices in the first column as is
without any reordering. For the case with periodic boundary condition (Fig.
3b), the $\hat{L}$ factor fills in near the diagonal and the lower edge in the
figure. The fill-in at the lower edge arises from the long-range coupling as
induced by the periodic boundary condition. By switching to the Dirichlet
boundary condition (Fig. 3e), the fill-in near the lower edge disappears and
the number of non-zero elements in the $\hat{L}$ factor significantly
decreases. Modifying the Dirichlet boundary condition further reduces the
fill-in and the number of non-zero elements (Fig. 3h).

The third column in Figure 3 shows the $\hat{L}$ factor as obtained by
factoring the corresponding matrices in the first column using a minimum
degree heuristics to reorder the matrix [22]. Compared with the second column,
the use of the minimum degree heuristics significantly reduces the number of
fill-in for all three cases. Also, the number of non-zero elements reduces as
we go from the case with the periodic boundary condition (Fig. 3c), to the
case with the Dirichlet boundary condition (Fig. 3f), to the case with the
modified Dirichlet boundary condition (Fig. 3i). The $\hat{L}$ factor shown in
Fig. 3i, which was obtained for the modified Dirichlet boundary condition
using the minimum degree heuristics, has 456 non-zero element, as compared
with the number of non-zero elements of 500 in the system matrix with the
periodic boundary condition. The result in Figure 3 provides a direct
demonstration that the use the modified Dirichlet boundary condition, in
combination with the minimum degree heuristics, can significantly reduce the
fill-ins.

Figure 3: In the first column of Fig. 3, we show the sparsity pattern of
$\hat{A}$ corresponding to the finite difference discretization of the
$10\times 10$ grid with the a) periodic, d) Dirichlet boundary, and g)
modified Dirichlet boundary with reduced connectivity at the boundary
respectively. In the second column, we show the corresponding $\hat{L}$
factors for factorizing these three matrices as is. In the third column, we
show corresponding $\hat{L}$ factors for factorizing these three matrices
using a minimum degree heuristic to reorder the matrices. The number of
nonzero elements (nnz) for each matrix is also provided in each figure.

##  5 Numerical Experiments

In this section, we show that the concept as described above can be applied to
typical simulations of electromagnetic waves. We consider both the driven case
with a source, and eigenmode calculations without a source. For both cases, we
show that choosing either the periodic or the Dirichlet boundary conditions
behind the PML does not affect the accuracy of the results, but the choice of
the Dirichlet boundary condition leads to better numerical performance.

###  5.1 Driven Problem

We solve Eq. (9) as in the previous section but now the interior domain
consists of $N_{x}\times N_{y}$ with $N_{x}=N_{y}=301$ grid cells surrounded
with a PML with a thickness of 30 cells to ensure minimal reflection [27]. The
domain has a length $a=4\lambda$ where $\lambda$ is the wavelength of the
source.

First, we confirm that the use of either the periodic boundary or the modified
Dirichlet boundary in conjunction with a PML indeed has no effects on the
solution within the interior domain. In Fig. 4 we show the result where the
interior domain consits of vacuum and we put a point source at the center of
the interior domain. For the two boundary conditions, the field distributions
inside the computational domain are almost identical, as can be seen by
comparing visually the field distributions shown in Fig. 4a and b, or by
exampling the field distributions along a line in the computational domain as
shown in Figure 4c.

Figure 4: a) Field solution for $H_{z}$ for Eq. (9) with the periodic
boundary. The periodic boundary is denoted by a red line. b) Field solution on
the same grid but with the Dirichlet boundary condition (shown in orange). c)
Comparison of the field distribution along the cyan lines in a) and b). The
transparent gray regions in (a) and (b) denote where the PML lies.

For the matrix $\hat{A}$ corresponding to the driven problem, we show in Fig.
5 the difference in factorizing the matrix with the periodic boundary vs the
Dirichlet boundary. In Fig. 5a and b, we show the fill-in of the $\hat{L}$
factor with the modified Dirichlet and periodic boundary respectively for the
system shown in Figure 4. The $\hat{L}$ factor from the modified Dirichlet
boundary is noticeably sparser than the $\hat{L}$ corresponding to the
periodic boundary.

In Fig. 5c, we vary the size the domain by keeping the domain to be a square
shape and varying $N_{x}=N_{y}$. We plot the percentage reduction in the
number of non-zero elements in the $\hat{L}$ factor as we change from the
periodic to the modified Dirichlet boundary condition. The percentage
reduction fluctuates significantly as $N_{x}$ varies, since the reordering
strategies are heuristic. But in general, we see significant reduction for all
$N_{x}$’s. On average the percentage reduction decreases as $N_{x}$ increases.
This is expected since the boundary, being a smaller and smaller fraction of
the domain as $N_{x}$ increases, would lead to a smaller effect on the
relative fill-in. What is remarkable is that the boundary couplings clearly
have a disproportionately large effect on the fill-in of $\hat{L}$ factor,
even though the boundary nodes are often just 1% or less of the total nodes in
the entire grid. Even with $N_{x}\approx 1000$, there are cases where the use
of Dirichlet boundary condition results in 33% fewer nonzeros in the $\hat{L}$
as compared with the use of the periodic boundary condition.

Figure 5: a) Fill-in pattern of the $\hat{L}$ factor for the system shown in
Fig. 4 with the Dirichlet boundary. b) Fill-in pattern for the $\hat{L}$
factor for the system with the periodic boundary. c) The percentage in the
reduction of the number of nonzeros in $\hat{L}$ factor in going from the
periodic to the modified Dirichlet boundary condition for varying values of
$N_{x}=N_{y}$.

###  5.2 Eigenvalue Problem

As a second example, we consider an eigenvalue problem that commonly arises in
the study of nanophotonic waveguides [28, 29]. In this problem, we consider a
structure that is periodic along the $x$-direction. For a given operating
frequency $\omega$, we solve for the propagating wavevector $k_{x}$. To
develop the formalism, we substitute via the Bloch theorem
$H_{z}(x,y)\rightarrow H_{z}(x,y)e^{ik_{x}x}$ in Eq. (9) to obtain:

$\frac{1}{\mu_{0}}\bigg{(}\partial_{x}\bigg{(}\frac{1}{\epsilon_{r}}\bigg{)}\partial_{x}+\partial_{y}\bigg{(}\frac{1}{\epsilon_{r}}\bigg{)}\partial_{y}+\omega^{2}\mu_{0}+ik_{x}\bigg{(}\partial_{x}\epsilon_{r}^{-1}+\epsilon_{r}^{-1}\partial_{x}\bigg{)}\bigg{)}H_{z}e^{ikx}=\frac{1}{\mu_{0}\epsilon_{r}}k_{x}^{2}H_{z}e^{ikx}$
(10)

which is a quadratic eigenvalue problem of the form:

$\hat{M}\mathbf{v}+\lambda^{2}\hat{K}\mathbf{v}+\lambda\hat{C}\mathbf{v}=0$
(11)

where $\lambda=k_{x}$ and:

$\hat{K}=\frac{1}{\mu_{0}}\epsilon_{r}$ (12)
$\hat{M}=\frac{1}{\mu_{0}}\bigg{(}\partial_{x}\bigg{(}\frac{1}{\epsilon_{r}}\bigg{)}\partial_{x}+\partial_{y}\bigg{(}\frac{1}{\epsilon_{y}}\bigg{)}\partial_{y}+\omega^{2}\bigg{)}$
(13)
$\hat{C}=\frac{i}{\mu_{0}}\bigg{(}\partial_{x}\epsilon_{r}^{-1}+\epsilon_{r}^{-1}\partial_{x}\bigg{)}$
(14)

Eq. (11) can be reformulated into an equivalent linear, generalized eigenvalue
problem of the form:

$\hat{D}\lambda\mathbf{w}=\hat{B}\mathbf{w}$ (15)

where $\mathbf{w}=[\mathbf{x},\lambda\mathbf{x}]$ with $\mathbf{x}$ denoting
the vector formed by $H_{z}(x,y)$ and where:

$\hat{D}=\begin{bmatrix}\hat{K}&\hat{0}\\\ \hat{0}&\hat{I}\\\ \end{bmatrix}$
(16)

and

$\hat{B}=\begin{bmatrix}\hat{C}&\hat{M}\\\ -\hat{I}&\hat{0}\\\ \end{bmatrix}$
(17)

After discretization using a computational domain consisting of a single unit
cell of the waveguide structure, $\hat{0}$ denotes a matrix of zeros of
dimension $N_{x}N_{y}$. In the waveguide eigenvalue problem, similar to the
discussion in Section (2.2), we are typically interested in modes with
wavevector near a particular value $\sigma$. Therefore, following the same
derivation as in Section (2.2), the matrix that needs to factored is
$\hat{B}-\sigma\hat{D}$ [17].

As an illustration we consider the waveguide structure schematically shown in
Figure 6a [30, 31]. The permittivity of the materials used is provided in the
caption of Figure 6a. The waveguide consists of two walls. Each wall consists
of a mixture of dielectric and a Drude metal. The thickness of each wall is
0.3 $\mu$m. The separation between the walls is 1.6 $\mu$m. The structure is
periodic along the $x$-direction. The periodicity along the $x$-direction is
0.2 $\mu$m. The metal region in the wall has a width along the $x$-direction
of 0.04 $\mu$m. Two unit cells are shown in Figure 6a. By operating this
system at $\lambda=2$ $\mu$m, the walls become near-completely reflective,
creating a hollow-core waveguide with light guided in the vacuum region
between the walls.

To simulate this structure we use a computational domain consisting of a
single unit cell shown in Figure 6a. We discretize the computational domain
with $N_{x}\times N_{y}$ grid cells where $N_{x}=100$ and $N_{y}=150$. Along
the $x$-direction a periodic boundary condition is imposed at the edges of the
computational domain. Along the $y$-direction we place PML regions at both
edges. The PML regions consist of 10 grid cells. Outside the PML region we
truncate the grid with either the modified Dirichlet or the periodic boundary
conditions.

First, we validate that the use of different boundary conditions backing the
PML does not affect the guided modes of the system. In Fig. 6b and c, we
compare the real part of the field profiles of the structure for a sample mode
with the same solved $k_{x}$ eigenvalue but with either the Dirichlet or
periodic boundary conditions. We see an excellent agreement between the modal
profiles. In Fig. 7, we compare the bandstructure of the waveguide with the
periodic and Dirichlet boundary conditions backing the PML in the
$y$-direction. To generate these bandstructures, we scan a discrete set of
$\omega$ values, and look for the $k_{x}$ eigenvalues closest to 0.
Additionally, due to spurious solutions that come from discretizing Eq. (10),
a mode-filter must be implemented to produce the band-structure. Fig. 7a is
the band structure with the Dirichlet boundary and Fig. 7b is the
corresponding one with the periodic boundary. Again, we see that the band
structures in both case are essentially identical, as expected.

Figure 6: a) Schematic of the waveguide structure and the computational
domain. The violet regions consist of a metal described by the Drude model
with $\epsilon(\omega)=1-\omega^{2}/(\omega_{p}^{2}-i\gamma\omega)$, with
$\omega_{p}=0.72\pi\times 10^{15}$ rad/s, $\gamma=5.5\times 10^{12}$ $s^{-1}$.
The yellow regions consist of a dielectric material with a relative
permittivity of $16$. The blue region is vacuum. The gray regions are the PML
regions. The structure is periodic along the $x$-direction. Plotted here are
two unit cells of the structure. b) $Re(H_{z})$ modal profile of the structure
with a Dirichlet boundary at the top and bottom edges of the computational
domain at an operating frequency corresponding to a free space wavelength of 2
$\mu$m. b) The corresponding result with the periodic boundary at the top and
bottom edges of the computational domain. Figure 7: Band structures obtained
with a) the modified Dirichlet boundary, b) with the periodic boundary
condition. The blue and red dots correspond to the real and imaginary parts of
the $k_{x}$, respectively. Figure 8: a) Sparsity pattern of
$\hat{B}-\sigma\hat{D}$ with $\sigma=0$ in Eq. (11) of the eigenvalue problem.
The gray lines partition the matrix into $N_{x}N_{y}\times N_{x}N_{y}$ sub-
blocks corresponding to those in Eqs. (17) and (16). The corresponding
sparsity pattern of the $\hat{L}$ factor with the periodic boundary condition.
c) The corresponding sparsity pattern of the $\hat{L}$ factor for the modified
Dirichlet boundary condition. In the second row, we show the analogous result
with $\sigma=1$ (in units of $\lambda_{p}$, the free-space wavelength at the
plasma frequency).

Now we compare the fill-in of factorizing $(\hat{B}-\sigma\hat{D})$ with the
periodic versus the modified Dirichlet boundary conditions. In Fig. 8a-c we
consider the case where $\sigma=0$. In Fig. 8a, we show the sparsity pattern
of the matrix $\hat{B}$ in Eq. (17). In Fig. 8b and c, we show the fill-in of
the $\hat{L}$ factor of $\hat{B}$ with the domain truncated by the periodic
boundary and the modified Dirichlet boundary conditions, respectively. The
fill-in patterns vary substantially between the two cases. Moreover, even
though we have PML regions only along the $y$ boundaries, we still achieve a
$40\%$ reduction in the number of nonzeros. In the second row of Fig. 8, we
show a case with $\sigma\neq 0$. Since $\hat{B}$ and $\hat{D}$ have different
sparsity patterns, the sparsity pattern shown in Fig. 8d differ from that in
Fig. 8a. In comparing Fig. 8e with Fig. 8f, we again see that the use of the
modified Dirichlet boundary condition leads to less fill-in.

##  6 Conclusion

In summary, we consider direct solvers for the FDFD methods applied to
Maxwell’s equations. We demonstrate that the choice of boundary condition used
behind the PML can have significant effects on the numerical efficiency of the
LU factorization of the system matrix. With the PML the physical properties of
the solutions do not depend on the choices of the boundary conditions. But we
observe that the use of a modified Dirichlet boundary condition can lead to a
reduction of fill-in by up to 40%, as compared to the commonly used periodic
boundary conditions. Such a reduction translates to significant increase in
computational speed and decrease in memory requirements for the direct
solvers.

While in our paper we focus on the PML, our results should be applicable for
other numerical absorbers, such as the adiabatic absorbers [32]. Also, while
we illustrate our observations with two-dimensional simulations, the same
observations should be applicable for three-dimensional simulations as well.
Our observation should also be relevant for finite-element frequency-domain
solvers.

Our result is important for enhancing the performance of direct solvers in the
FDFD methods. In addition, this result should have implications for several
numerical tasks in computational electromagnetics that can benefit from
efficient LU decomposition of the system matrix, such as adjoint variable
method[12, 33], efficient line searches for gradient-based photonic
optimization [34, 35, 36], and domain decomposition methods utilizing the
Schur complement approach [37, 38].

## References

  * [1] Jean-Pierre Berenger.  A perfectly matched layer for the absorption of electromagnetic waves.  Journal of Computational Physics, 114(2):185–200, oct 1994.
  * [2] Allen Taflove.  Computational Electrodynamics.  Artech House, Inc., USA, 3rd edition, 2005.
  * [3] J.-P. Berenger.  Perfectly matched layer for the fdtd solution of wave-structure interaction problems.  IEEE Transactions on Antennas and Propagation, 44(1):110–117, 1996\.
  * [4] Z.S. Sacks, D.M. Kingsland, R. Lee, and Jin-Fa Lee.  A perfectly matched anisotropic absorber for use as an absorbing boundary condition.  IEEE Transactions on Antennas and Propagation, 43(12):1460–1463, 1995.
  * [5] Steven G. Johnson.  Notes on perfectly matched layers (pmls), 2021.
  * [6] F Collino and P B Monk.  Optimizing the perfectly matched layer.  Computer Methods in Applied Mechanics and Engineering, 164(1):157–171, 1998.
  * [7] Arti Agrawal and Anurag Sharma.  Perfectly matched layer in numerical wave propagation: factors that affect its performance.  Appl. Opt., 43(21):4225–4231, Jul 2004.
  * [8] Ardavan Oskooi and Steven Johnson.  Distinguishing correct from incorrect pml proposals and a corrected unsplit pml for anisotropic, dispersive media.  Journal of Computational Physics, 230:2369–2377, 04 2011.
  * [9] J.-P. Berenger.  Perfectly matched layer for the fdtd solution of wave-structure interaction problems.  IEEE Transactions on Antennas and Propagation, 44(1):110–117, 1996\.
  * [10] Wonseok Shin and Shanhui Fan.  Choice of the perfectly matched layer boundary condition for frequency-domain Maxwell’s equations solvers.  Journal of Computational Physics, 231(8):3406–3431, apr 2012.
  * [11] Kane Yee.  Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media.  IEEE Transactions on Antennas and Propagation, 14(3):302–307, 1966\.
  * [12] Georgios Veronis, Robert W. Dutton, and Shanhui Fan.  Method for sensitivity analysis of photonic crystal devices.  Optics Letters, 29(19):2288, oct 2004.
  * [13] Randall LeVeque.  Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems (Classics in Applied Mathematics Classics in Applied Mathemat).  Society for Industrial and Applied Mathematics, USA, 2007.
  * [14] G. Strang.  Introduction to Linear Algebra.  Wellesley-Cambridge Press, 2016.
  * [15] Xiaoye S. Li.  An overview of SuperLU: Algorithms, implementation, and user interface, sep 2005.
  * [16] Timothy A. Davis.  Algorithm 832: UMFPACK V4.3 - An unsymmetric-pattern multifrontal method.  ACM Transactions on Mathematical Software, 30(2):196–199, jun 2004\.
  * [17] R. B. Lehoucq, D. C. Sorensen, and C. Yang.  Arpack users guide: Solution of large scale eigenvalue problems by implicitly restarted arnoldi methods., 1997.
  * [18] Gene H. Golub and Henk A. van der Vorst.  Eigenvalue computation in the 20th century.  Journal of Computational and Applied Mathematics, 123(1):35–65, 2000.  Numerical Analysis 2000. Vol. III: Linear Algebra.
  * [19] Carmen Campos and Jose E Roman.  Strategies for spectrum slicing based on restarted Lanczos methods.  Numerical Algorithms, 60(2):279–295, 2012.
  * [20] Chang Peng and Daniel Boley.  Methods for large sparse eigenvalue problems from waveguide analysis.  In Proceedings of the 1996 12th Annual Review of Progress in Applied Computational Electromagnetics, pages 645–652, January 1996.
  * [21] Mihalis Yannakakis.  Computing the minimum fill-in is np-complete.  SIAM Journal on Algebraic Discrete Methods, 2(1):77–79, 1981.
  * [22] Patrick R. Amestoy, Timothy A. Davis, and Iain S. Duff.  An approximate minimum degree ordering algorithm.  SIAM Journal on Matrix Analysis and Applications, 17(4):886–905, 1996.
  * [23] Alan George.  Nested dissection of a regular finite element mesh.  SIAM Journal on Numerical Analysis, 10(2):345–363, 1973.
  * [24] John R Gilbert and Robert Endre Tarjan.  The analysis of a nested dissection algorithm.  Numerische Mathematik, 50(4):377–404, 1986.
  * [25] Pratik Agarwal and Edwin Olson.  Variable reordering strategies for slam.  In 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 3844–3850, 2012.
  * [26] George Karypis and Vipin Kumar.  MeTis: Unstructured Graph Partitioning and Sparse Matrix Ordering System, Version 4.0.  http://www.cs.umn.edu/~metis, 2009.
  * [27] Allen Taflove, A. Oskooi, and S. Johnson.  Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology.  Artech House, 01 2013.
  * [28] Georgios Veronis and Shanhui Fan.  Guided subwavelength plasmonic mode supported by a slot in a thin metal film.  Opt. Lett., 30(24):3359–3361, Dec 2005.
  * [29] Andre Nicolet and Christophe Geuzaine.  Waveguide propagation modes and quadratic eigenvalue problems.  In 6th International Conference on Computational Electromagnetics, pages 1–3, 2006.
  * [30] Nathan Z. Zhao, Ian A. D. Williamson, Zhexin Zhao, Salim Boutami, and Shanhui Fan.  Penetration depth reduction with plasmonic metafilms.  ACS Photonics, 6(8):2049–2055, 2019.
  * [31] Nathan Zhao, Ian A. D. Williamson, Zhexin Zhao, Salim Boutami, and Shanhui Fan.  Penetration depth engineering in plasmonic metafilms for enhanced reflection and confinement.  In Conference on Lasers and Electro-Optics, page FF3E.7. Optica Publishing Group, 2020.
  * [32] Ardavan F. Oskooi, Lei Zhang, Yehuda Avniel, and Steven G. Johnson.  The failure of perfectly matched layers, and towards their redemption by adiabatic absorbers.  Opt. Express, 16(15):11376–11392, Jul 2008.
  * [33] Christopher M. Lalau-Keraly, Samarth Bhargava, Owen D. Miller, and Eli Yablonovitch.  Adjoint shape optimization applied to electromagnetic design.  Optics Express, 21(18):21693, sep 2013.
  * [34] Salim Boutami, Nathan Zhao, and Shanhui Fan.  Large permittivity increments for efficient predictive photonic devices optimization.  In Connie J. Chang-Hasnain, Andrei Faraon, and Weimin Zhou, editors, High Contrast Metastructures IX, volume 11290, pages 55 – 67. International Society for Optics and Photonics, SPIE, 2020.
  * [35] Salim Boutami, Nathan Zhao, and Shanhui Fan.  Determining the optimal learning rate in gradient-based electromagnetic optimization using the shanks transformation in the lippmann–schwinger formalism.  Opt. Lett., 45(3):595–598, Feb 2020.
  * [36] Nathan Z. Zhao, Salim Boutami, and Shanhui Fan.  Efficient method for accelerating line searches in adjoint optimization of photonic devices by combining schur complement domain decomposition and born series expansions.  Opt. Express, 30(4):6413–6424, Feb 2022.
  * [37] Nathan Z. Zhao, Salim Boutami, and Shanhui Fan.  Accelerating adjoint variable method based photonic optimization with Schur complement domain decomposition.  Optics Express, 27(15):20711, jul 2019.
  * [38] Nathan Zhao, Sacha Verweij, Wonseok Shin, and Shanhui Fan.  Accelerating convergence of an iterative solution of finite difference frequency domain problems via schur complement domain decomposition.  Opt. Express, 26(13):16925–16939, Jun 2018.


# Investigating the Relationship Between Dropout
Regularization and Model Complexity in Neural Networks

Jai Sharma, Milind Maiti, Christopher Sun

(August 2021)

## Abstract

Dropout Regularization, serving to reduce variance, is nearly ubiquitous in
Deep Learning models. We explore the relationship between the dropout rate and
model complexity by training 2,000 neural networks configured with random
combinations of the dropout rate and the number of hidden units in each dense
layer, on each of the three data sets we selected. The generated figures, with
binary cross entropy loss and binary accuracy on the z-axis, question the
common assumption that adding depth to a dense layer while increasing the
dropout rate will certainly enhance performance. We also discover a complex
correlation between the two hyperparameters that we proceed to quantify by
building additional machine learning and Deep Learning models which predict
the optimal dropout rate given some hidden units in each dense layer. Linear
regression and polynomial logistic regression require the use of arbitrary
thresholds to select the cost data points included in the regression and to
assign the cost data points a binary classification, respectively. These
machine learning models have mediocre performance because their naive nature
prevented the modeling of complex decision boundaries. Turning to Deep
Learning models, we build neural networks that predict the optimal dropout
rate given the number of hidden units in each dense layer, the desired cost,
and the desired accuracy of the model. Though, this attempt encounters a
mathematical error that can be attributed to the failure of the vertical line
test. The ultimate Deep Learning model is a neural network whose decision
boundary represents the 2,000 previously generated data points. This final
model leads us to devise a promising method for tuning hyperparameters to
minimize computational expense yet maximize performance. The strategy can be
applied to any model hyperparameters, with the prospect of more efficient
tuning in industrial models.

## Introduction

Dropout Regularization is a technique used in Deep Learning and other branches
of Machine Learning to combat overfitting on training data with respect to the
validation set [2]. Another hyperparameter, the architecture of a neural
network, determines the complexity of the functions the model will be able to
learn, with a deeper model that contains more hidden units able to reduce bias
in the network [6].

Superficially, the implementation of dropout regularization seemingly
counteracts the effect of building a model with more hidden units. While
increased model complexity mitigates bias, higher dropout rates serve to
restore bias to prevent overfitting. Thus, it is crucial for the Machine
Learning practitioner to equilibrate these two hyperparameters to find the
tuning that induces the finest performance. The common approach is to build a
complex model accompanied by substantial regularization. However, the problem
with this attempted one-size-fits-all approach is that training the model can
require extensive computational resources and prolonged train times, not to
mention the greater issue of suboptimal convergence.

Thus, we seek to identify the relationship between the dropout rate and the
number of hidden units in a neural network. We aim to demonstrate brief
hyperparameter tuning with respect to cost and accuracy, specifically
intending to illustrate that adding complexity in a neural network is not
always the steadfast solution to unlock a better-performing model. Rather, a
simpler model, which can be trained in much less time, can offer similar
performance. We experiment with the relationship between variance and bias by
studying model complexity and dropout regularization, showing that multiple
model configurations can result in similar performance. In addition, we
explore a method to find the best dropout rate to use in a Deep Learning model
given the model complexity of each dense hidden layer. Finally, we share a
generalized version of this insightful procedure that Machine Learning
practitioners can utilize to more quickly perform hyperparameter tuning to
balance bias and variance while optimizing validation cost and accuracy.

### Data Availability Statement

For experimentation, we used three publicly available cardiovascular disease
data sets. Each of the training examples in these three data sets represents
individual patients. Medical information about the patients and a binary
diagnosis label for cardiovascular disease are included for each data example.
We selected these data sets because they are all relatively small, allowing
for faster hyperparameter tuning. The following describes each of the data
sets used (data sets will henceforth be referred to by the names below in
bold):

##### High-Level Data Set

The High-Level Data Set is available on the Kaggle [9] database and contains
70,000 examples and 11 high-level features, most binary and ternary. Prior to
experimentation, this data set underwent slight preprocessing to remove
statistical outliers using the 1.5 x IQR Rule with a modified coefficient of
2.5, which we believed better represented the spread of the data set. Training
examples with features beyond the range permitted by the IQR Rule were
excluded from the data set.

##### Cleveland Data Set

The Cleveland Data Set is a famous source of cardiovascular disease data
containing 303 examples and 75 features available on the UCI Machine Learning
Repository [3]. Since only 14 of the 75 features are used by Machine Learning
researchers when studying cardiovascular disease, we also experimented with
this subset of features.

##### Combined Data Set

We assembled a data set that is a fusion of cardiovascular disease data
sourced from hospitals in Cleveland, Hungary, Switzerland, and Long Beach,
California. The individual data sets are available on the UCI Machine Learning
Repository [3]. The Combined Data Set contains 1,025 examples with the same 14
features included in the Cleveland Data Set.

## Methods

### Measuring Performance of Deep Learning Models Configured With Unique
Combinations of Dropout Rate and Layer Size

To investigate the relationship between dropout regularization and model
complexity, we explored hyperparameter tuning. We randomly chose combinations
of model complexity and regularization, then trained 2,000 different neural
networks 111All Deep Learning models constructed during our research were
trained using the highly optimized TensorFlow Keras library [4]. using these
randomly chosen hyperparameters to predict the cardiovascular disease
diagnosis label given the medical information of the patients. Finally, we
determined the binary cross entropy cost and binary accuracy of each model
after 150 iterations of training. The architecture of the networks comprised
six dense hidden layers each with the same number of hidden units and the same
dropout rate. However, each model trained used a unique number of hidden units
and a unique dropout rate.

Let x be the number of hidden units in each of the six hidden layers of a
constructed Deep Learning model.
---

Let y be the dropout rate for this Deep Learning model.
---

The value of x ranged from $2^{3}=8$ to $2^{10}=1024$. To sample as many
models with $2^{3}=8$ to $2^{4}=16$ hidden units as models with $2^{9}=512$ to
$2^{10}=1024$ hidden units, we chose x based on a logarithmic scale.
Specifically, $\textit{{x}}=\lfloor 2^{c}\rfloor$, where $c$ was a uniformly
random number chosen from the interval $(3,10)$. (The number x must be an
integer, so the floor function was performed.) The value of y was a uniformly
random number chosen from the interval $(0,1)$.

All other hyperparameters remained constant during experimentation: each
model’s weights were initialized with Xavier initialization using the same
random seed, and each model was optimized with Adam Optimization [5] using a
mini-batch [7] size of 128 examples, undergoing 150 epochs of training.

Figure 1: Generalized structure of each of the 2,000 neural networks built,
where X is the number of hidden units in each hidden layer and Y is the
dropout rate.

After training 2,000 such models on each of the three data sets and recording
the performance metrics of each model, we proceeded to plot all generated
points on three-dimensional graphs: the number of hidden units on the $x$
axis, the dropout rate on the $y$ axis, and the attained cost/accuracy on the
$z$ axis.

Figure 2: Cost plots generated by training 2,000 models on each of the three
data sets, where the color map is determined by the cost, the z-axis
parameter. Notice that all plots have distinct contours that are able to be
modeled. Figure 3: Accuracy plots generated by training 2,000 models on each
of the three data sets, where the color map is determined by the binary
accuracy, the z-axis parameter.

### Modeling a Generalized Relation Between Dropout Rate and Layer Size

Analyzing the figures above, we identified respective regions of lower cost
and higher accuracy. For example, we observed a trough spanning the length of
the cost plot of the High-Level data set that looked linear in nature.
Considered together, our observations compelled us to answer the question:
What is the correlation between dropout rate and model complexity that unlocks
superior performance? We hypothesized the optimal dropout rate to be a
function of the number of hidden units in a hidden layer (assuming the number
of hidden layers upon other variables stayed constant), and thus built a
series of predictive models to fit the collected 2,000 points of data and
quantify the hypothesized relationship.

##### Linear Regression

Linear Regression was the first technique utilized to fit points in the
“trough” regions of all graphs. Initially, the boundaries of these troughs
needed to be specified, requiring the manual setting of numerical and
percentile thresholds of the 2,000 cost data points. Combinations of dropout
rate and model complexity that achieved a cost lower than this numerical
threshold or within this percentile threshold were included in the regression,
while other inferior combinations were excluded. We experimented with
numerical thresholds incrementing in 0.001 from 0.54 to 0.55 and percentile
thresholds incrementing in 5% from 5% to 25%. Clearly, the selection of
thresholds was inherently subjective, not to mention that abundant amounts of
usable data were being wasted because of the procedure’s selective nature.

##### Polynomial Logistic Regression

To utilize all data points collected, we turned to Logistic Regression to
assign models a binary classification based on their cost. Because cost is on
a continuous scale, we transformed cost into a binary feature by labeling
costs under a hand-chosen percentile threshold as 1 (superior) and all other
costs as 0 (inferior). We tested 10% and 25% as these percentile thresholds,
after which we built Logistic Regression models to classify combinations of
dropout rate and number of hidden units depending on their achieved cost.
Feeding second-order and third-order polynomial features into the models
allowed for the learning of nonlinear decision boundaries.

Though the Logistic Regression models employed all 2,000 generated data
points, the arbitrary selection of a percentile for the classification of
costs was still required. In addition, the selection of the degree of
polynomial features was not methodical. Furthermore, the Logistic Regression
models were quite naive in the sense that they assigned the trained Deep
Learning models just two elementary labels, which meant they were incapable of
differentiating between a very poor model and a below average model or between
sufficient model and an exceptional model.

#### Neural Networks

##### 1\. Neural Networks for Optimal Dropout Rate Prediction

Deficiencies of previous machine learning models led us to utilize the complex
hypothesis functions of neural networks [6]. With these models, we wished to
rigorously quantify the performance of model using the exact cost rather than
using discrete labels. Accordingly, we constructed three neural networks
222These neural networks are not to be confused with the networks trained on
the cardiovascular disease data or the network that outputs a continuous
dropout rate. that received the number of hidden units, the desired cost, and
the desired accuracy as inputs. At test time, the desired cost was interpreted
as the minimum cost attained in our generated data set, and the desired
accuracy was interpreted as the maximum accuracy attained in our generated
data set. Thus, we were able to feed an array of hidden units ranging from
$2^{3}$ to $2^{10}$ into these neural networks to obtain the optimal dropout
rate for each value. Instead of yielding a binary label, the network yielded a
continuous output between 0 and 1, representing the optimal dropout rate given
the inputs.

Figure 4: Generalized architecture of the neural network trained on the 2,000
generated points to predict optimal dropout rate given a number of hidden
units in each dense layer, the desired cost, and the desired accuracy.

The suboptimal performance of these neural networks gave us insight on how to
implement the same type of model more effectively.

##### 2\. Neural Networks for Surface Plots

We turned to a more rational approach, training neural networks to predict the
cost of a Deep Learning model given the number of hidden units in the hidden
layers and the dropout rate of the Deep Learning model.

Two neural networks were trained for each of the three data sets to model cost
and accuracy, these networks containing six hidden layers with 16 hidden units
each and regularized with a dropout rate of 0.1. We created three-dimensional
surface plots which were representative of the continuous cost predictions
given hyperparameter combinations. These decision boundaries were juxtaposed
with the generated data points to visualize the models’ performance. (These
plots are shown in the Results section.)

## Results

### Performance of Deep Learning Models With Unique Hyperparameter
Combinations

Analyzing the generated data led to various findings about the nature of the
interactions between dropout rate, number of hidden units in a hidden layer,
and the cost and accuracy achieved. Some figures, such as the accuracy plot of
the Cleveland Data Set, depicted succinct regions in which the cost or
accuracy was determined entirely by only a single hyperparameter of the two we
analyzed. Other figures, such as the cost plot of the High-Level Data Set,
showed that using a more hidden units with a higher dropout rate yielded the
same performance as using less hidden units with a lower dropout rate. In
contrast, still other figures, such as the accuracy plot of the Combined Data
Set, favored the use of a lower dropout rate without regard for the number of
hidden units.

Such findings compelled us to rethink common presupposed relationships between
the these two hyperparameters, challenging the frequent assumption that a
larger number of hidden units in a dense layer simply requires a larger
dropout rate to enhance performance. In addition, the plots discredited the
notion that increasing the number of hidden units results in better
performance. In fact, many spatial regions of the highest cost were achieved
by models that had a relatively large amount of hidden units. We urge network
designers to consider that an increase in model complexity, though
theoretically reducing bias, can actually lower the performance of a model.
The involuted balance created by regularization parameters such as the dropout
rate cannot be readily oversimplified into a questionable cliche that adding
hidden nodes unlocks superior performance.

### Generalized Relation Between Dropout Rate and Layer Size

##### Linear Regression

When examining the performance of the linear models which were fitted on a
subset of the generated data to predict the optimal dropout rate, we noticed
that while effective for the High-Level data set, linear models were unable to
perform well on the Combined data set and the Cleveland data set, particularly
due to the non-linearity of the problem and the disarray of the selected
points. In addition, with the threshold percentile being set at the relatively
low level of 25%, 75% of the data or 1500 points in our case were unused,
amounting to waste in computational resources. Arbitrary selection of the
threshold percentile was another undesirable byproduct of this model, as it
was unclear which threshold percentile would lead to better performance.

Mean Absolute Error of Linear Regression
---
|  High-Level Data Set |  Cleveland Data Set |  Combined Data Set
10% Threshold |  0.0609 |  0.0976 |  0.0768
25% Threshold |  0.0942 |  0.1084 |  0.1137
Figure 5: Visualizations of Linear Regression on all three data sets, with
different percentile thresholds being used to select points included in the
regression. Notice that this approach had different levels of effectiveness
across data sets and performed especially poorly on the Combined Data Set.

##### Polynomial Logistic Regression

The Logistic Regression models greatly improved upon the performance of the
Linear Regression due to the extra data employed. When using Logistic
Regression, the non-linearity of separating optimal models from non-optimal
models became apparent, and thus we employed higher-order features. When
utilizing these higher-order features, we found that including third-order
features greatly improved performance without overfitting the data. Logistic
Regression with higher order features yielded promising results for
distinguishing “superior” models from “inferior” models. Nonetheless, Logistic
Regression still fell short in that it offered a binary classification of the
data. Due to the absence of clear definitions of optimal and non-optimal
models, arbitrary thresholds had to be selected. The more desirable solution
would be to utilize the continuous cost values.

Binary Accuracy of Polynomial Logistic Regression (25% Threshold)
---
|  High-Level Data Set |  Cleveland Data Set |  Combined Data Set
Second Degree Features |  0.8490 |  0.8695 |  0.9825
Third Degree Features |  0.9005 |  0.8870 |  0.9820
Figure 6: Visualizations of the Polynomial Logistic Regression decision
boundaries on all three data sets, using a percentile threshold of 25%. Red
areas are regions that the models classify as having a superior cost (“good”,
1), while the background color represents the models’ classification of
inferior cost (“bad”, 0). Blue points symbolize combinations of
hyperparameters that achieve costs within the 25% threshold, while green
points achieve costs that exceed the threshold.

##### Neural Networks for Optimal Dropout Rate Prediction

The relationship between the number of hidden units in each dense layer and
the optimal dropout rate, as modeled by the first type of neural network
introduced in the Methods section, is plotted below.

Mean Absolute Error of Neural Network Used to Predict Dropout Rate
---
High-Level Data Set |  Cleveland Data Set |  Combined Data Set
0.1384 |  0.0995 |  0.1640
Figure 7: These figures plot the number of hidden units in each dense layer on
the $x$ axis against the predicted optimal dropout rate on the $y$ axis.

When we plotted the decision boundaries of the the neural networks in three
dimensions, we saw that the models had disreputable performances, which caused
us to speculate that this was due to the fundamental nature of the idea
itself. We realized that theoretically, more than one dropout rate can be
configured with the same number of hidden units to achieve the same cost and
accuracy. Thus, we ventured to discover why the dropout rate is not always a
function of the number of hidden units, cost, and accuracy. To explain this
mathematically:

Let $n$ be the number of hidden units in each hidden layer for some Deep
Learning model.
---

Let $p$ be the optimal dropout rate given $n$.
---

Let $\epsilon_{1}$ and $\epsilon_{2}$ both be some positive numbers.
---

By definition, the configuration of $p$ dropout rate achieves the lowest
possible cost for all possible models with $n$ hidden units in each hidden
layer. (Side Note: This configuration is not guaranteed to achieve the highest
possible accuracy, only the lowest possible cost.)

Using $p+\epsilon_{1}$ or $p-\epsilon_{2}$ as the dropout rate will result in
a higher cost than using $p$ as the dropout rate. The cost function graph,
though not appearing as completely continuous, strongly suggests that there
exists some values for $\epsilon_{1}$ and $\epsilon_{2}$ such that using
$p+\epsilon_{1}$ or $p+\epsilon_{2}$ as the dropout rate will result in the
same cost. (If not exactly the same, then our data indicates that
$\epsilon_{1}$ and $\epsilon_{2}$ can be found such that the resulting costs
have a negligible difference.) This observation explains why it is impossible
to perfectly predict the optimal dropout rate given the number of hidden units
in each of the hidden layers, the cost of the model, and the accuracy of the
model.

##### Network Networks for Surface Plots

These neural networks were significant because they enabled us to test
thousands of configurations of hyperparameters through forward propagation in
a matter of seconds, demonstrating an immense reduction in time of the
hyperparameter tuning process.

Mean Absolute Error of Neural Networks for Surface Plots
---
|  High-Level Data Set |  Cleveland Data Set |  Combined Data Set
Cost Points |  0.0115 |  0.0897 |  0.1156
Accuracy Points |  0.0100 |  0.0459 |  0.0479

Figure 8: Visualization of the decision boundaries of the neural networks used
to model the 2,000 generated points, creating surface plots.

## Discussion

Allowing more flexibility in designing network architectures opens up many new
avenues for the network designer. Our findings regarding the effects of
regularization and hyperparameter tuning should motivate network designers to
prioritize the conservation of valuable computational resources by designing
more efficient architectures.

Rapid prototyping in industrial applications can also be made more accessible
when designers are not as inclined to tediously tune bigger models, but
instead implement a larger volume of smaller models. While each neural network
used in our research was small and was trained in roughly 10 minutes on GPU
runtime, most industrial Deep Learning models utilize massive data sets with
numerous features to tackle difficult problems (e.g. image recognition,
sentiment analysis, self-supervised learning tasks). In these cases, training
many models in a short period of time and iterating through the development
process is essential to creating deployable products. Hence, the results of
our research call for a more insightful approach to tune hyperparameters in
neural networks.

### A Rapid and Effective Hyperparameter Tuning Method

The following is an outline of a new method to accelerate the meticulous
process of hyperparameter adjustment in Deep Learning models.

Let n be the initial number of models trained on random hyperparameter
settings chosen according to their respective range of values (in our case n =
2,000; the number of hidden units in each dense layer was chosen on a
logarithmic scale from $2^{3}$ to $2^{10}$ and the dropout rate was chosen on
a linear scale from 0 to 1).
---

The n models with can be trained and the random hyperparameter settings, the
cost, and the final accuracy after training can be recorded.
---

Using this generated collection of points, a network designer can then train a
separate neural network (or some other Machine Learning model as they see fit)
to fit either the cost or the accuracy, in essence to predict the
cost/accuracy of the model given the value of each hyperparameter.
---

The network designer can use this previous model to sample thousands of
hyperparameter configurations and generate a surface plot.
---

Using this surface plot, the designer can “zoom” into areas which are
predicted by the model to have relatively lower costs. Note that using this
surface plot is theoretically equivalent to testing out thousands of
configurations using forward propagation in a matter of seconds.
---

The network designer could then repeat this process by training another n
models initialized with random hyperparameter settings sampled from this
smaller value range, iterating through the process until a model with the
least cost is reached.
---

### Misconceptions

##### Misconception 1:

Some may suppose that even with using the proposed method, the designer still
needs to sample many configurations of hyperparameters. However, the method
ensures that the designer samples fewer points than would be sampled with
traditional hyperparameter tuning, yet still has the same result.

##### Misconception 2:

Some may suppose that training the neural network to generate the surface plot
will take more time than simply trying new configurations. However, even while
n grows very large, the neural network can be easily trained because there are
only two features and one label.

### Generalizations and Implementation Details

As one might imagine, this hyperparameter tuning process is not restricted to
only the hyperparameters we discussed in this example. Rather, the process can
be extrapolated to a larger scale containing different hyperparameters or more
hyperparameters, allowing designers to more quickly determine the
configurations of their models.

We remind developers while using this technique to consider the trade off
between time and accuracy that is intrinsic to this process, in that the
generation of more points expends more computational resources but ultimately
results in a better network, while training less points, though timely, may
result in a flawed hyperparameter selection. Symbolically, if n is too small,
then the surface plot will not generalize to untrained points and the neural
network will fail to its job. On the other hand, if n is too large, then this
defeats the purpose of using the proposed method altogether. One potential
solution to calibrate the number of models trained is to incrementally
decrease the number of configurations tested. For example, a designer could
test 100 configurations, zoom in, test another 10 configurations, zoom in,
test another 5 configurations, and so on.

## Need For Further Research

Further experiments can be conducted to determine if the phenomena we observed
exist in larger data sets and diverse data sets not limited to cardiovascular
disease. The Deep Learning models we trained with 150 epochs can instead be
trained for a larger number of iterations to verify the usability of our
results before extrapolation to industrial models. Different variations of
hyperparameters, such as varying the number of hidden units and dropout rate
between each dense layer of a model, can also be tested. This is crucial
because in more sophisticated applications of Deep Learning (e.g.
convolutional neural networks), different layers of the network learn either
“specific” or “general” features, making it practically impossible for each
layer to have the same architecture [10]. Thus, more experimentation needs to
be done to gauge the usefulness of the procedures when applied to different
types of Deep Learning models.

In addition, research can be done to gauge the accuracy and reliability of the
surface plot neural networks by training many models with random
hyperparameters and comparing the results to those predicted by the surface
plot neural networks. The significant trade off between surface plot
generalizability and computational expense should also be explored to validate
or nullify our proposed method of hyperparameter tuning, and to ultimately
find the answer to the optimal number of models to train before generating a
surface plot.

## References

  * [1] Glorot, X. and Bengio, Y. “Understanding the difficulty of training deep feedforward neural networks.”
In International Conference on Artificial Intelligence and Statistics, 2010.

  * [2] Hinton, G.E., Srivastava, N., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R.R. “Improving neural
networks by preventing co-adaptation of feature detectors.” arXiv preprint
arXiv:1207.0580, 2012.

  * [3] Janosi, A., Steinbrunn, W., Pfisterer, M., Detrano, R., and Aha, D.W. “Heart Disease Data Set.”
UCI Machine Learning Repository, archive.ics.uci.edu/ml/datasets/Heart+Disease

  * [4] Keras Team.“Keras Documentation: The Sequential Model.” Keras,
keras.io/guides/sequential_model/.

  * [5] Kingma, D. P. and Ba, J. L. “Adam: A Method for Stochastic Optimization.”
arXiv preprint arXiv:1412.6980, 2015.

  * [6] Larochelle, H., Bengio, Y., Louradour, J., and Lamblin, P. “Exploring strategies for training deep neural
networks.” In The Journal of Machine Learning Research, pages 1–40, 2009.

  * [7] Li, M., Zhang, T., Chen, Y., and Smola, A.J. “Efficient mini-batch training for stochastic optimization.”
In Proc. of the 20th ACM SIGKDD international conference on Knowledge
discovery and data mining,
pages 661–670, 2014.

  * [8] Ng, A. and deeplearning.ai. “Deep Learning Specialization.”
coursera.org/specializations/deep-learning?utm_source=gg.

  * [9] Ulianova, S. “Cardiovascular Disease Dataset.” 2019,
kaggle.com/sulianova/cardiovascular-disease-dataset.

  * [10] Yosinski, J., Clune, J., Bengio, Y., and Lipson, H. “How transferable are features in deep neural networks?”
In Advances in Neural Information Processing Systems, pages 3320–3328, 2014.

### Acknowledgements

We would like to thank Mr. Govind Tatachari for his tremendous help and
mentorship over the course of this research. We very much appreciate his
valuable inputs and feedback.

## Appendix

##### Cross-Entropy Cost Function

When training the Deep Learning models, we assess model performance using the
Cross-Entropy Cost Function, defined as the following:

$J(\theta)=\frac{1}{m}*\sum^{m}_{i=1}\left[y^{(i)}\log{(h(x^{(i)}))}+(1-y^{(i)})\log{(1-h(x^{(i)}))}\right],$

where $\theta$ contains the parameters for the neural network, $m$ is the
number of examples in the validation data set, $y^{(i)}$ is the ground truth
label for example $i$, $x^{(i)}$ is the input for example $i$, $h(x^{(i)})$ is
the prediction for example $i$, and $\log$ is the natural logarithm with base
$e\approx 2.71828\ldots$

##### Binary Accuracy

Let $m$ be the number of examples in the validation data set. Let $x$ be the
number of patients who are correctly classified by the Machine Learning model.
The Binary Accuracy is then defined as:

$\boxed{\text{Binary Accuracy}=\left(100\cdot\frac{x}{m}\right)\%}.$

Binary Accuracy ranges from 0% to 100%, where 100% corresponds to perfect
model performance on validation data.