Heuristic Methods and Performance Bounds for Photonic Design

In the photonic design problem, a scientist or engineer chooses the physical parameters of a device to best match some desired device behavior. Many instances of the photonic design problem can be naturally stated as a mathematical optimization problem that is computationally difficult to solve globally. Because of this, several heuristic methods have been developed to approximately solve such problems. These methods often produce very good designs, and, in many practical applications, easily outperform `traditional' designs that rely on human intuition. Yet, because these heuristic methods do not guarantee that the approximate solution found is globally optimal, the question remains of just how much better a designer might hope to do. This question is addressed by performance bounds or impossibility results, which determine a performance level that no design can achieve. We focus on algorithmic performance bounds, which involve substantial computation to determine. We illustrate a variety of both heuristic methods and performance bounds on two examples. In these examples (and many others not reported here) the performance bounds show that the heuristic designs are nearly optimal, and can considered globally optimal in practice.


Introduction
An important part of photonics, and many other scientific and engineering fields, is the design and construction of physical devices.A 'physical device', as used in this paper, includes anything as simple as a spherical lens, where a scientist can easily find an optimal lens for a given application by using basic algebra and ray optics, all the way to potentially very complicated devices and applications such as range detection and mapping using Li-DAR [YSC + 20], where building an 'optimal' device, in nearly any practically useful sense, is an open research problem.Though the two examples we give are in the field of photonics (and this is, indeed, the focus of our review), the term refers to any device that can, at least theoretically, be built, and whose desired behavior can be mathematically specified.
Traditionally, the design of physical devices was done by an engineer or scientist, whom we will generally call a designer, for a specific application.The designer would have a library, either physically or through experience, of well-understood components or materials, each performing a specific function.These components would then be carefully pieced together, often with a good amount of ingenuity, in order to perform the desired task.In many cases, the resulting designs could then be modified in part or in whole, or combined with other designs, in order to perform even more complicated functions.This procedure, while effective in practice, is time consuming and sometimes even tangential to the final application of the design itself.
A second approach to constructing physical devices was initially explored in the early 1960s within the field of electrical engineering, originally for the purpose of recognizing printed letters [KL63] and has since been extended to many other fields [LAJ + 98, GCL + 03, EK10, ASKA18, MLP + 18], with sometimes very surprising results [Tho97].In this approach, the designer specifies a mathematical objective function which, given a design, outputs a number representing how well the input design matches the desired specifications; the lower this number is, the better the design.This function is then fed into an optimization algorithm, which attempts to minimize this objective function by finding a design that is good in the sense specified by the designer.In almost all practical cases, the algorithm will fail to find the best possible design, and, except in very specific scenarios, may never do so even when left to run for a very long time.But, in many applications, the resulting designs found often have much better performance than any design found by humans.This approach can be seen as a declarative approach to design: the user specifies what they want, while ceding control of how it should be done to the optimization algorithm.This idea has many names in different fields, with field-specific connotations and denotations; these include 'automated design', 'computational design', 'inverse design' (in photonics and aerospace engineering), 'shape design' and 'generative design' (in mechanical engineering), 'topology design' (in several fields), or 'synthesis' (in hardware design), among many others.We will simply call this optimization problem the 'physical design problem', with the understanding that many, if not all, of the previously mentioned problems are instances of the physical design problem.

The physical design problem
The usual physical design problem can be formally stated in many ways.In this review, we focus on a simple but general formulation, which, as we show in this section, includes many important problems in photonic design.

Physics
The design problem starts with a physical theory that describes the behavior of the field (which we will write as z) under some excitation (which we will write as b).The field z and excitation b are vectors in some (typically infinite dimensional) vector space.We focus here on the case when the physical theory is linear, in which case we can write the physics equation as where A is a linear operator.Almost all problems in physical design are governed by linear physics equations, e.g., Maxwell's equations, Helmholtz's equation, the heat equation, the Schrödinger equation, among many others.
Electromagnetic wave equation.For example, a common way of writing the electromagnetic (EM) wave equation in terms of the electric field E and the currents J, for a monochromatic wave with angular frequency ω is [JJWM08, §2], where ε denotes the permittivities at each point in space and µ 0 is the magnetic permeability, which we assume to be constant throughout space, in this example.We can then make the following correspondences: which naturally leads to an equation of the form of (1).
Discretization.We will work with an appropriate discretization of the field, excitation, and physics equation (1).We will overload notation to use the same symbols for their discretized versions.In the sequel, the field z will be a vector in R n , the excitation b will be a vector in R m , and the linear operator A will be a matrix in R m×n .The physics equation (1) is then a set of m linear equations in n scalar variables.Complex fields and excitations can be reduced to the real case by separating them into their real and imaginary parts.
Solutions and simulations.For a fixed A and excitation b, we will call any z which satisfies (1) a solution of the physics equation.In general that there can be a unique solution, many solutions, or no solution.We will focus on the case when there is a unique solution, i.e., m = n and A is invertible, so z = A −1 b.We refer to computing the field z = A −1 b, given A and b, as a simulation.There are many simulation methods, including generic methods for solving linear equations such as sparse-direct methods [ABLM19,CDHR08] or iterative methods [HS52,SS86], and custom methods crafted specifically for the particular physics equations [Gil11].We note that, in practical photonic design, the resulting linear systems can be very large, with the number of variables, n, often in the millions or tens of millions.Solving systems in the upper end of this range often requires the use of large-scale linear solvers [SVS + 20].
Approximate solutions and physics residual.Some of the methods we will see work with approximate solutions of the physics equation.For any field vector z and excitation b, we define the physics residual as r phys = Az − b.A reasonable numerical measure of the size of the residual is Ax − b / b , where • is a norm, typically the Euclidean norm • 2 .Simulations, especially those that use iterative methods, produce fields with small physics residuals.
Modes.A simple trick can be used to represent the modes of a system as a solution to (1).Suppose the original (discretized) physics equation is Hz = λz, where λ is the eigenvalue and z is an associated mode.Directly expressing this as Az = b with A = H − λI and b = 0 yield a physics equation that is not invertible and has multiple solutions (including, of course, z = 0).To fix a unique solution, we use a linear normalization and insist that c T z = 1, where c ∈ R n is some nonzero vector.We then represent the mode equation and normalization as Az = b with This has a unique solution, provided c is not an eigenvector of H with eigenvalue λ.Note that the linear normalization in (3) differs from the usual choice of normalization, z 2 = 1, where • 2 is the Euclidean norm.

Design parameters
In physical design, the designer is able to change the system physics equation (1), by choosing some parameters that affect the physics, i.e., A and b.Thus A and b depend on some design parameters θ ∈ R d .For example, in photonic design, θ is generally a variable that controls the permittivities inside of the device.We can then write the physics equation (1) with explicit dependence on the design parameters as When A(θ) is invertible we have z(θ) = A(θ) −1 b(θ); i.e., the field also depends (implicitly) on the design parameters θ.
Affine physics design.In many practical cases (and in all of the examples we will show), A and b are affine functions of the design parameters; i.e., where A i ∈ R m×n and b i ∈ R m .These A i are usually sparse matrices and vectors; that is, each design parameter θ i affects just a few entries of A and b.Because equation (4) is affine in θ when holding z fixed, and affine in z when holding θ fixed, this type of equation is sometimes called 'bi-affine' or 'multilinear' in θ and z.

Diagonal physics design.
A common and useful special case of (4) is when b(θ) is a constant, b(θ) = b 0 , and A(θ) can be written as where diag(θ) is a matrix whose diagonal entries contain the elements of the vector θ and is zero elsewhere.In other words, A i = E ii , where E ii is the matrix with only one nonzero entry (which is one), in the i, i entry.We will call this special case of (4) the diagonal physics equation.
EM wave equation in diagonal form.A specific example of a diagonal physics equation are the EM wave equation (2), whenever a designer is allowed to vary the permittivities.In this case, we can define θ to be proportional to the permittivities (with proportionality constant µ 0 ω 2 ) such that the following correspondences can be made: This correspondence results in an equation of the form of (5).(We can also similarly write the more general case where the designer is allowed to vary both the permittivities and permeabilities, in this form.See, e.g., the appendix in [AVB19].) Low-rank updates.Whenever A(θ) is sparse (which is very common in practice), it is often possible to compute an explicit factorization of A(θ) (e.g., the sparse Cholesky factorization, when A(θ) is positive semidefinite) which makes evaluating A(θ) −1 b(θ) inexpensive, after the factorization.Additionally, whenever the matrices A i are also low rank (as in the diagonal case (5), for example), updates to the factorization of A(θ) can be efficiently computed [CDHR08].This implies that we can also efficiently evaluate A(θ ) −1 b(θ ), when only a small number of entries of θ differ from those of θ, given the factorization for A(θ).
Parameter constraints.In general, a designer has constraints on the parameters that can be chosen.Because of this, we will define the feasible parameter set, Θ ⊆ R d , such that only design parameters satisfying θ ∈ Θ are feasible or valid.In many applications, Θ is a hyperrectangle (or box) indicating that each component of θ must lie in some interval given by where the inequalities are elementwise.
In photonic design, it is often not possible to vary the permittivities along an interval but are instead allowed to be one of two possible values.This leads to another common set of parameter constraints, where each component of θ is constrained to be exactly one of two elements (we will call this class of constraints Boolean constraints): In this class of constraints, the total number of parameters that are feasible is |Θ| = 2 d .The feasible parameter set Θ can also include fabrication constraints such as minimum possible feature sizes, among other possibilities [PPSV17], but we will focus on the common cases of box or Boolean constraints.
Normalization.We can re-parametrize the design parameters to lie between −1 and 1 (or any other limits).So, without loss of generality, we can always consider the upper and lower bounds to be θ min = −1 and θ max = 1, where 1 is the vector with all entries equal to one.To do this, we introduce a new parameter δ ∈ R d and define where θ = (θ max + θ min )/2 is the parameter midpoint, while ρ = (θ max − θ min )/2 is the parameter radius, and • denotes the elementwise (Hadamard) product.The constraint θ ∈ Θ becomes −1 ≤ δ ≤ 1, where the inequalities are elementwise, in the box-constrained case, or δ ∈ {−1, 1} d in the Boolean case.We then have A(θ) = Ã(δ), b(θ) = b(δ), with Ã and b affine, and

Optimization problem
The design objective is often written as a function of the fields, specifying how well the resulting field matches the desired objective.This objective function could specify the power in a given direction, the field overlap (i.e., the inner product between the current field and a desired one), or the total energy, all of which can be written as functions depending only on the field, that the designer may wish to optimize.For example, the designer may wish to maximize the power transmitted through a specific port of a device at a given frequency [PLL + 15], or to maximize the focusing efficiency of a lens within a specific region [CM20, BPC + 20].Finding a device that best matches this objective can be directly phrased as a mathematical optimization problem.
Objective function.In other words, we seek to find a design whose field optimizes an objective function f : R n → R ∪ {+∞}.The function's input is a field z (generated by some design θ ∈ Θ) and its output is a number that specifies how good or bad this field is, or how well the field matches the designer's specification.Without loss of generality, we will assume that a higher number is worse (i.e., a designer wishes to minimize f ), but we can just as well maximize f by, equivalently, minimizing its negative, −f .We allow the objective function f to take on infinite values to denote hard constraints on the desired field: if f (z) = +∞ for some field z, then z is not a feasible field.
Problem statement.We can then compactly write the problem that a designer wishes to solve (or approximately solve), which we will call the physical design problem: where the problem variables are the fields z ∈ R n and design parameters θ ∈ R d , while the problem data include the matrix A(θ) ∈ R m×n , the excitation b(θ) ∈ R m , and the parameter constraint set Θ ⊆ R d .We will call the special case of (6) where the physics equation is the diagonal physics equation (5) the diagonal physical design problem.(As a reminder, in photonic design, θ is usually proportional to the permittivities, while A(θ) is the operator corresponding to the electromagnetic wave equation (2).) Problem attributes.Problem (6) has several important properties.In many practical cases, the function f is a convex function and the set Θ is a convex set, which implies that problem (6) is a convex problem in the variable z, when holding θ fixed, and is a convex problem in θ, when holding z fixed.This property leads to some useful heuristics for approximately solving problem (6); cf., [LV10,LV13].Additionally, problem (6) often has a smooth (differentiable) objective function f , while Θ can almost always be represented as a number of smooth equality and inequality constraints, which happens in many practical applications and all examples presented in this review.In this case, we can apply general nonlinear optimization solvers such as IPOPT [WB06] directly to problem (6).
Computational hardness.On the other hand, it is not difficult to show that even finding a feasible design and field for problem (6) is, in general, a computationally difficult problem (i.e., it is NP-hard) even when the design parameters θ are unconstrained; that is, even if Θ = R n .To do this, we will reduce the subset sum problem [Kar72], a problem known to be NP-hard, to an instance of (4).This would imply that, if we could efficiently solve problem (4), then we could efficiently solve the subset sum problem, which is widely believed to be computationally hard to solve.(See, e.g., [Aar11] for a good overview of P vs. NP and its implications.) The subset sum problem asks: given c ∈ R n , is there a nonzero binary vector x ∈ {0, 1} n such that c T x = 0? We will show that, given c, we can answer this question by finding a field z and a design θ that satisfy the constraints of (6), which will imply that, in general, problem (6) is computationally difficult.
First, note that we can write the following conditions on z and θ, as an instance of (4), by appropriately stacking the conditions into a matrix form.Now, the last two conditions are true if, and only if, z i (z i − 1) = 0, which also happens only when z i ∈ {0, 1} for i = 1, . . ., n.The second condition implies that 1 T z = 0, and, when combined with the first condition, this statement is true if, and only if, there exists a nonzero solution to the subset sum problem, with the provided vector c.This, in turn, shows that finding a feasible design θ and field z such that the constraints of (6) are satisfied must, in general, be a problem that is at least as hard as the subset sum problem.Because problem ( 6) is likely to be computationally difficult to solve exactly when the number of parameters is large, we will focus on heuristics which approximately solve the problem for the remainder of this paper.
Multi-scenario design.A common design task is to find a single device that has good performance across many different scenarios.For example, the device might need to be robust against temperature variations, or the device might be required to filter out a number of specific wavelengths, while allowing others through.In this case, we will assume that the device satisfies S instances of the physics equation ( 1), where the design parameters θ are held fixed across the S instances, but the physics equation or the excitation, is allowed to vary: Here A s (θ) ∈ R ms×ns , b s (θ) ∈ R ms and z s ∈ R ns for s = 1, . . ., S. This leads to the multi-scenario physical design problem: where the variables are the fields z s ∈ R ns in each of the s = 1, . . ., S scenarios and the design parameters θ ∈ R d , while the objective function f : R on the fields at any of the S scenarios.
In fact, it turns out that we can write any instance of the multi-scenario physical design problem (7) as an instance of (6).To do this, we collect all of the individual physical equations into a single constraint by placing all of the A s (θ) along the diagonal of a lager matrix A(θ), and stacking the excitations b s (θ) and fields z s .More specifically, define So we can write the S distinct physical equations, A s (θ)z s = b s (θ) for scenarios s = 1, . . ., S, as the single equation, Overloading notation slightly, such that f (z) = f (z 1 , . . ., z S ), then reduces problem (7) to one of the form of (6).In other words, it suffices to only consider a problem of the form of (6).
Eliminating the field variables.A simple (and relatively common) equivalent formulation of problem ( 6) is to note that, because we have assumed A(θ) is invertible, we can write it as a problem that depends only on the design variables; i.e., since z = A(θ) −1 b(θ), we can write problem (6) as the following optimization problem over θ ∈ R d : When the function f is differentiable, we can easily compute its derivatives with respect to each component of θ: where z = A(θ) −1 b(θ) is the solution to (4) for design parameters θ.We can write this in a slightly more compact form by defining y = A(θ) −T (∇f (z)), such that Note that, to find z and y, we only need to solve two systems of linear equations, one over A(θ), and one over A(θ) T .This observation can be used to efficiently compute the gradient of the objective with respect to the design variables and is called the adjoint method [CLPS03].
Other methods of computing the derivative when, for example, the matrix A(θ) is not invertible include automatic differentiation through the simulation (as in [MWA + 20]) among others.
Eliminating the design variables.In the case of diagonal design, and for some choices of the parameter constraint set Θ, it is also possible to eliminate the corresponding field variables [AVB20].For example, in the case where the parameter constraint set Θ is a box, −1 ≤ θ i ≤ 1, then, given a field z, there exist design parameters θ satisfying (i.e., z satisfies the diagonal physics equation ( 5) with parameters θ) if, and only if the field z satisfies where the absolute value | • | is taken elementwise.
To see this, note that, from (5), we can write and, since we are free to choose any −1 ≤ θ ≤ 1, we have that such a θ exists, if, and only if, This lets us write the diagonal physical design problem in the following equivalent way: where the only problem variable is the field z ∈ R n .We will call this formulation of the diagonal physical design problem, as in [AVB20], the absolute-upper-bound formulation.This rewriting also shows an interesting property of the physical design problem, whenever f is convex: if you know the signs of any optimal field, then the design problem becomes convex and therefore easy to solve globally.That is, if we know s = sign(z ), where sign is the signum function and where z is optimal for (10), then any solution to the convex optimization problem with variable z ∈ R n is a solution to (10) with the same optimal value as z .This follows from two basic facts.First, any z that is feasible for (11) is feasible for (10) with the same objective value because And, second, that z is feasible for (11) since and, by definition, z is feasible for (10) and so satisfies |A 0 z − b 0 | ≤ |z |.This observation also leads to an optimization algorithm which iteratively updates the signs and yields a sequence of feasible fields with decreasing objective value [AVB20, §3], called sign-flip descent or SFD for short.We compare its performance against other basic solvers in §3.
A similar analysis holds when Θ = {−1, 1} n , i.e., θ is constrained to be Boolean.In this case, we have that, given some field z, there exists a Boolean design θ ∈ {−1, 1} n such that θ and z satisfy the diagonal physics equation if, and only if, z satisfies We note that, in photonics, the box-constrained formulation sometimes leads to designs that cannot be practically implemented because the resulting designs may have permittivities that lie along an interval (i.e., θ i may lie anywhere in [−1, 1]), while most fabrication methods only allow the use of two possible permittivities (i.e., we must have θ i ∈ {−1, 1}).Despite this, the box-constrained formulation is often used as a good initialization for current inverse design algorithms, which then approximately solve the Boolean case [SVS + 20].
Approximate solution methods.There are many practical methods for approximately solving (6).For example, many of the earliest solution methods approximately solve the physical design problem by applying zeroth order (or 'derivative free') optimization algorithms after eliminating the field variable, as shown in problem (8) [Hau95, JB00, CBT01, KBDQ04, SHL + 04, HSS05, GFEV07].Such methods are easy to implement in practice, as they only require the ability to perform a basic simulation; i.e., to solve the physics equation (4) for a given θ ∈ Θ. Zeroth order optimization methods include hill-climbing, genetic algorithms [HGLL06], simulated annealing [KGV83], Nelder-Mead [NM65], and adaptive coordinate descent [LSS11], among many others [RS13].While effective at finding designs with moderate to good performance, zeroth order optimization methods scale poorly and suffer from slow convergence when compared to higher-order methods.(See [CSV09] and [AH17] for more information on zeroth order optimization methods.) A second important family of optimization algorithms, which include the algorithms most used in practice, are the first order optimization algorithms, which are also almost always applied to problem (8).In these cases, such methods additionally make use of gradient information (9), leading to better computational performance and faster convergence times, at the expense of higher implementation complexity, as these methods require additional information from the simulator.Examples of first order optimization algorithms used in practive include L-BFGS-B [BLNZ95], proximal gradient methods [Par14], and the method of moving asymptotes [Sva87], among many others.(See [Ber16, ES19, KW19] for a comprehensive overview.) We compare a few different methods in §3.1 in terms of computational performance and resulting design performance.

Performance limits
Any approximate optimization method for the physical design problem (6) can be used to generate approximately optimal designs.In other words, if we let p be the optimal value for (6), these procedures generate a design and field that satisfy the physics equation (1) whose objective value, say p, satisfies p ≥ p .In general, because problem (6) is hard to solve, it is hard to know how far away our designs are from the true optimal value p .For example, once we have approximately optimized (6) and received some design with objective value p, it is not clear if this design is close to optimal (and no design can do significantly better) or if there are designs that have much better performance than the one we've found.
It is an old tradition in physics to then ask: what is the best possible value that we can hope to achieve?More specifically: is there some lower bound d such that we can guarantee that the optimal value of problem (6), p , is never smaller than this bound; i.e., p ≥ d?Such a bound can be interpreted in many ways.For example, it can be interpreted as an 'impossibility result', which states that no device that satisfies the physics equation (4) and the parameter constraints, θ ∈ Θ can have objective value smaller than d.We can also interpret is as a 'certificate of optimality': given some design with objective value p, if p is close to d, then p must also be close to the optimal objective value, p , since p ≥ p ≥ d.
Of course the best performance bound is p , but computing this is intractable, and we seek bounds that can be computed at reasonable cost.
Additionally, performance bounds can be very important in speeding up the design process.For example, it is often not clear how large a design needs to be in order to achieve reasonable performance.This often results in designers having to experiment with the total device size in order to find a design which has at least the desired performance.Lower bounds on these values, if they are efficiently computable, would give an indication of how large a design needs to be in order to achieve the designer's goals without additional (potentially very computationally expensive) experimentation.
Methods.Roughly speaking, there are two main approaches to the problem of finding lower bounds to the optimal objective value of problem (6).The first, and likely the earliest of the methods, is to make basic physical assumptions about the system (for example, that the system size is substantially smaller than the wavelength of the excitation [Whe47,Chu48]) and derive bounds on the corresponding quantities [MBB82].While these methods are historically important and yield good rule-of-thumb heuristics for design, many of the bounds derived in this way require assumptions that are not satisfied by the devices found by inverse design, or result in weak bounds.The second approach, which has become relatively popular recently (see, e.g., [MPHR + 16, AVB19, MJVR19, TAS + 20, MCJR20, MCR20b, MCR20a, MVJR20, MBY + 20, KM20, GSJC20, CM20]), essentially uses basic properties of the constraints and objective function of problem (6) to derive bounds on the best possible performance of the problem.Such approaches include algebraic manipulations of the physics equation (1) combined with the parameter constraints θ ∈ Θ, and applications of Lagrange duality to problem (6).The resulting bounds often do not have analytical forms, but can be numerically evaluated by an efficient algorithm and are therefore called computational bounds.We will discuss such bounds in this section.

Lagrange duality
The basic tool in a number of these bounds is the use of Lagrange duality.The idea is as follows.Given the optimization problem for some objective function f : R n → R, constraint function h : R n → R m , and optimization variable x ∈ R n , we form the Lagrangian: where λ ∈ R m , with λ ≥ 0 is a Lagrange multiplier or dual variable.This lets us define the Lagrange dual function g : R m → R, given by See, e.g., [BV04, §5] for more information on Lagrange dual functions.
Lower bound property.The function g has a few interesting properties.First, for any λ ≥ 0, g(λ) is always smaller than p , the optimal value of (13), and is therefore a performance bound.More specifically, we have The first inequality follows from the fact that the set of x which satisfy h(x) ≤ 0 is no larger than the set of all x ∈ R n , and the second follows from the fact that, because h(x) ≤ 0 and λ ≥ 0, then λ T h(x) ≤ 0.
In general, evaluating g(λ) at some λ ≥ 0 is at least as hard as solving the original problem (13).On the other hand, when it is possible to efficiently evaluate g(λ), it is almost always possible to efficiently find the optimal value of the dual problem (14) because the function g is always a concave function, even when the objective function f and constraints h in the original problem are not convex [BV04, §5.1.2].
Initializations.A solution to the dual problem (14) often suggests a good initialization for heuristics which attempt to minimize problem (13).Given some dual variable λ that is optimal for (14), there exists some x 0 which minimizes the Lagrangian at this choice of dual variable under some basic assumptions on the objective function f and constraints h.In practice, x 0 is generally close to a reasonable design (see, e.g., [AVB19]) even if it is not feasible; i.e., it need not satisfy h(x 0 ) ≤ 0, except in some special scenarios such as when the functions f and h are both convex, in which case x 0 is globally optimal [BV04, §5.2].Because it is often true that x 0 is easy to evaluate whenever g(λ ) is easy to evaluate, this initialization can be seen as a by-product of finding a solution to (14).

Local power conservation
The first approach to constructing bounds for (6), presented originally in some generality in [MPHR + 16] and then [MVJR20] and later extended and fully clarified in [GSJC20, KZM20, MCJR20] and subsequently fully generalized in [MCR20b] and [KM20] in the case where the parameters are Boolean (i.e., Θ = {−1, 1} n ).We will present a further generalization to the case where the parameters are box-constrained (Θ = [−1, 1] n ), by a slightly different proof that considers a relaxation of the absolute-upper-bound formulation (10).
(The Boolean case follows from an identical argument by considering (12), instead.)Bounds of this form are essentially power conservation laws over a given subdomain, which, in photonics, are often included under the name 'optical theorem' [New76].
Power inequalities.Starting with the absolute-upper-bound formulation given in (10), we can square both sides of the inequality constraint to receive an equivalent formulation, In other words, a field z is feasible (i.e., there exists a design θ ∈ Θ such that z and θ satisfy the diagonal physics equation ( 5)) if, and only if the following quadratic inequalities are all satisfied: (a We can, of course, multiply these inequalities by a nonnegative value λ i ≥ 0 for i = 1, . . ., n and add any number of them together to get another valid, quadratic inequality, where a T i is the ith row of A 0 .This can be more compactly expressed as where D = diag(λ) and has nonnegative entries along the diagonal.The family of inequalities (17), parametrized by the nonnegative diagonal matrices D is the same as the family given in [MCR20b,KM20], except in the case where Θ specifies a box constraint, instead of a Boolean one.Additionally, because z satisfies (16) if, and only if, it satisfies (17) for all diagonal matrices D with nonnegative diagonals, then it follows that the family of quadratic inequalities in [MCR20b, KM20] is tight, in the following sense: any field z which satisfies these power conservation laws for all nonnegative diagonal matrices D must also have a corresponding design θ ∈ Θ such that z and θ simultaneously satisfy the diagonal physics equation ( 5).(The inequalities in (17) are of a slightly different form than those presented in [MCR20b,KM20].We show their equivalence in appendix A.) Relaxed formulation.Because any feasible field z satisfies (17) for any nonnegative diagonal matrix D, we can relax the family of inequalities (17) from all nonnegative diagonal matrices D to a finite number of them, which we will write as D j for j = 1, . . ., N .The following problem can then be seen as a relaxation of (15): with the field z ∈ R n as the only variable and problem data A 0 ∈ R n×n , the excitation b 0 ∈ R n , and the diagonal matrices D j ∈ R n×n with nonnegative entries along the diagonal, for j = 1, . . ., N .Note that, if D j = E jj for j = 1, . . ., n, where E jj is the matrix with a single nonzero in the j, j entry (which is one) and is zero elsewhere, we recover the original problem (15).Additionally, while this problem is a relaxation of the original, it is still likely to be computationally hard to solve.
Quadratic objective.In general, finding lower bounds to the relaxed formulation (18) also need not be easy, except in some important special cases.For example, when the objective function f is a quadratic, for some symmetric matrix P ∈ S n , vector q ∈ R n , and r ∈ R, then problem (18) is a quadratically constrained quadratic program (QCQP, see [PB17]), and the dual function g corresponding to this problem has a closed form solution.In fact, the dual problem of ( 18) is, in general, a convex semidefinite program (SDP) which can be solved in practice for moderate values of n and N [BV04, §1], allowing us to find a lower bound to (18) and therefore to (6) efficiently.
Dual problem.To find the dual problem, we first formulate the Lagrangian of (18): where λ ∈ R N + and we have defined (This follows from expanding the expression and collecting the quadratic, linear, and constant terms in the variable z.)It is not hard to show that the dual function is, for λ ≥ 0, (See, e.g., [PB17, §3.2].)Here T (λ) + is the Moore-Penrose pseudoinverse [BV18, §11.5] of T (λ), while T (λ) ≥ 0 means that T (λ) is positive semidefinite, and R(T (λ)) is the range of T (λ).The corresponding problem of maximizing g over λ ∈ R N + can be written as a standard form SDP, which we will call the power dual bound : maximize with variables t ∈ R and λ ∈ R N .This problem can be easily specified using domain-specific languages such as CVXPY or JuMP.jl and solved using convex solvers that support SDPs, such as Mosek [ApS19], SCS [OCPB16], or COSMO.jl[GCG19].Because of the sparsity of A 0 , it is often the case that such problems have chordal structure [VA15], which can be exploited to efficiently solve problem (19) by some solvers such as COSMO.jl.
Optimal choice of D j matrices.Given the dual problem (19) of problem (18), the question remains of how to best choose the diagonal matrices D j for j = 1, . . ., N such that the optimal value of ( 19) is maximized.Historically, the initial bounds in [MPHR + 16, MVJR20, GSJC20, KZM20, MCJR20] assumed a number of fixed diagonal matrices D j .
Later, [MCR20b] and [KM20] generalized the approach to include any diagonal matrix D and [KM20] proposed an iterative algorithm which, starting with some diagonal matrix D 1 (such as D 1 = I) would solve a problem similar to (19) at iteration k and propose a new diagonal matrix D k+1 that would be appended to the constraints of (18).This method is conceptually similar to the cutting-plane method described in [PB17, §3.5] applied to the dual problem (19) with D j = E jj for j = 1, . . ., n. (More accurately, the procedure proposed in [KM20] solves an SDP relaxation of (18) instead of the dual problem (19), but it can be shown that both problems have the same optimal value by strong duality [PB17, §3.3].)It is also reasonable to ask: what are the best possible choices of D j such that the optimal value of (19) is maximized?It is not hard to show that a single (correctly chosen) diagonal matrix, D = diag(λ), suffices and that this matrix can be efficiently found.To see this, note that we can choose D j = e j e T j for j = 1, . . ., n such that, This would let us write We then note that picking λ optimal for (19), when it exists, gives a diagonal matrix D = diag(λ ).Additionally, solving (18) with this choice of matrix D is a special case where the number of quadratic constraints, N , equals 1, which can be efficiently solved and has the same optimal value as ( 19).(See, e.g., appendix B of [BV04].)

Diagonal physics dual
Another approach to computing lower bounds for the diagonal physical design problem ( 6) is by a direct application of Lagrange duality, originally given in [AVB19] and later extended in [TAS + 20] and [ZZM20].This approach gives a lower bound when the objective function is separable: and the constraints are of the following form: where S k ⊆ {1, . . ., d} for k = 1, . . ., K are disjoint sets, specifying which indices are constrained to be equal.We additionally note that a lower bound for this constraint set also yields a lower bound for the Boolean case, since this constraint set contains the Boolean one.
Problem Lagrangian.The basic trick here is to rewrite problem (6) as the following (equivalent) problem: minimize where we have pulled out the constraint θ ∈ Θ into the indicator function I : R n → R∪{+∞} for the set Θ, We can then easily formulate the Lagrangian of this problem: which we note is separable in terms of z: Dual function.As before, the dual function is defined as Computing the inner infimum is relatively simple, which gives: Here, f * i : R → R ∪ {+∞} is the convex conjugate of f i (sometimes called the Fenchel conjugate function or simply the conjugate function) defined as and is well known for a number of functions [BV04, §3.3].Additionally, we will make use of the fact that f * i is always a convex function, even when f i is not convex.To compute the outer infimum, we note that we can write where we have used the fact that − inf x h(x) = sup x −h(x) for any function h, and the fact that a scalar convex function achieves its maximum over an interval at the boundary of that interval.
Dual problem.Given the dual function g, we can find a lower bound to the original problem by evaluating g for any ν ∈ R n .We can then ask what's the best possible dual bound, which gives the following dual problem, which we will call the diagonal dual bound : This problem is a convex optimization problem with variable ν ∈ R n , whose optimal value, d , can almost always be efficiently found whenever the function f * can be efficiently evaluated.
Field bounds.In some cases, the bounds given by ( 23) can sometimes be very weak; e.g., if f i = 0 for all indices i except one.One way of improving the lower bounds is to add a redundant constraint to (22); i.e., a constraint such that, if z satisfies (A 0 + diag(θ))z = b 0 , then it also satisfies h(z) ≤ 0 for some function h : R n → R, such that the dual function of the resulting problem is still simple to evaluate.One useful example, originally presented in [TAS + 20], is to note that, if z satisfies where we have taken the norm of both sides of the expression, and the inequality follows from the triangle inequality.Using the fact that where σ 1 (A 0 ) is the smallest singular value of A 0 , we find or, after some rearrangement: Evaluating the dual for this new problem is a simple extension of the procedure given above as in [AVB19, §6].In fact, the same procedure can be extended to any norm • by replacing σ 1 (A 0 ) with 1/ A −1 0 , where A −1 0 x , but the resulting dual function need not be easy to evaluate.
Mode volume.Another extension, presented originally in [ZZM20] is for an objective function f of the form whenever z i = 0 and is +∞ otherwise.Here, i is a fixed index and we will assume there is no excitation; i.e., b 0 = 0. Note that this objective function is similar to, and can be easily extended to include, the cavity mode volume.Because b 0 = 0, then the physics equation and objective are 0-homogeneous in z, so any feasible point z with z i = 0 can be scaled by a nonzero value η ∈ R, such that ηz is also feasible with the same objective value, f (ηz) = f (z).We can then fix a normalization by setting z i = 1 to get an equivalent problem: minimize with variables z and θ and some fixed index i.The dual for problem (24) can be computed using the same method presented in this section.

Numerical results: device designs and bounds
In this section, we show a basic comparison between some of the heuristics presented in §1.3 and the performance bounds presented in §2 on a small problem of designing a Helmholtz resonator in one and two dimensions.In these cases, we can certify that the heuristics find a device whose performance is at most 2% above the global optimum, even though the original problem is likely hard.We then show examples in the literature of much larger devices for which no lower bound has been computed (and, indeed, cannot be computed with the current methods in a reasonable amount of time), but has resulted in practically useful devices whose performance is much better than that of traditional designs.

Performance bounds of small designs
In this section, we will compare the performance of three heuristics, sign-flip descent (SFD), IPOPT, and genetic algorithms (GA), and the lower bounds presented in §2.We will first compare their respective performance on a small, one-dimensional design, where both the power lower bounds and the dual lower bounds can be computed in a reasonable amount of time, and then compare their performance on a larger two-dimensional design.The times reported in this section are from a dual-core 2015 MacBook Pro laptop running at 2.9GHz with 8GB of RAM.At the time of writing, this machine is 5 years old and is roughly representative of the computational power available in more recent standard and lower-end consumer laptops.
Problem formulation.In both scenarios, our objective function f will be to best match a desired field ẑ, f (z) = z − ẑ 2 2 , where the desired field is a cosine wave with a Gaussian envelope on the left half of the domain, and is equal to zero on the right half.The excitation is a single delta function in the center of the domain.In this problem, the designer is then allowed to choose the speed of the wave at each point in space (via the design parameters θ), such that the objective function is minimized.(For more details and code, see the appendix B.) 1D problem.We compare the objective performance and computational performance of the heuristics and bounds in table 1.We also plot the corresponding fields and the desired field in figure 1.While GA gives rather poor solutions (which only somewhat match the largest features of the desired field), IPOPT and SFD give approximately-optimal fields that are essentially indistinguishable from the desired one.
We find, at least in this small scenario, that IPOPT and SFD have objective values that are extremely close to that of the power bounds and diagonal dual bounds, which means that the designs found in this scenario can, for all intents and purposes, be considered globally optimal.We also note that GA, while simple to implement, does not find a good solution even with some amount of tuning, while also having the worst performance of the available

Algorithm
Objective heuristics in terms of total time taken.Additionally, while the power bound was slightly tighter than the diagonal dual lower bound, it took nearly 90 times longer to converge for this small problem.This is due to the fact that SDPs solution times scale approximately cubically with the problem dimension (i.e., are O(n 3 )), which quickly becomes an issue for larger problems.
In fact, we note that it was difficult to find a desired field ẑ where the bounds and the performance of designs found by SFD differed significantly.We encourage readers to search for some cases where this is true by trying a few different desired fields in the code available for this paper.For more details, see appendix B.1.
2D problem.In the 2D problem, we again test the performance of GA, IPOPT, and SFD, and the diagonal design lower bound, though we do not compute the power bounds.We note that naïvely attempting to compute the power bounds in (19) by framing the problem as an SDP (using an SDP solver that does not support chordal sparsity) results in a dense matrix variable of size (251 2 ) 2 /2 = 251 4 /2, which requires approximately 16GB of memory to store, nearly double the available memory (8GB), not counting the additional memory required to perform operations on this matrix.
The results of this comparison are available in table 2 and the fields of the approximately optimized designs are shown in figure 2. We terminated any algorithm whose runtime was longer than 30 minutes on the current computer.We note that, again, SFD has surprisingly good performance, and, when combined with the lower bound, yields a design that is guaranteed to be no more than (11.9/11.7 − 1) ≈ 1.7% suboptimal, relative to the globally optimal value (in fact, as shown in figure 2, the resulting field from SFD and the desired field, ẑ, are difficult to visually distinguish).Additionally, solving the diagonal dual bound was faster than getting an approximate design from any of the heuristics, though we note that this is likely not the case with much larger designs or better-optimized heuristics.For more details, see appendix B.2.

Discussion
. The examples in this section show that modern heuristics have very good practical performance when compared to the available lower bounds.In fact, we suspect that this is true more generally: modern heuristics with reasonable initializations likely return designs that are very good, if not globally optimal, even when there exist no bounds that can certify this, or when the available bounds are very weak.Additionally, while we have made some basic performance optimizations to the available code, we opted for clarity as opposed to pure performance in the implementation of the algorithms presented, so the numbers above should only be interpreted as general guidelines.The code to generate the plots and tables is available at https://github.com/cvxgrp/pd-heuristics-and-bounds.

Practical inverse design examples
In this section, we show practical examples where the devices found by heuristics have been fabricated and experimentally verified.While the current bounds cannot be used to certify that the performance of these designs is close to globally optimal in a reasonable amount of time, the resulting designs have much better performance than that of traditional designs.For a comprehensive overview of the history of inverse design and its applications in practice, including older literature, we refer the reader to [MLP + 18].
Splitters.Perhaps the most striking applications of physical design is in the design of compact splitters-devices which, given some input in a specific scenario (for example, given an input at a specific frequency) must direct as much of the input as possible into a desired output location.Different scenarios would direct the input to different locations; i.e., the input is 'split' to different outputs depending on the scenario.Some examples of such devices and their performance can be found in, e.g., [PLL + 15, SWPM15, MSJ + 16]. Figure 3 shows an example of the splitter designed and fabricated in [SPS + 18].In this figure, three different inputs are represented by the three different colors (blue, red, and green), which are initially 'mixed' in the input on the left-hand-side of the domain.The device then separates the wavelengths into each of the three output channels.
Particle accelerators.There are other types of devices which have been designed and manufactured in practice.A recent demonstration of a laser-driven particle accelerator small enough to be placed on a chip [SYV + 20] required several components to be designed from scratch, since many did not have traditionally-designed counterparts.Additionally, the components had to be efficient enough that the input power did not destroy the material out of which they were made.The necessary components included the devices which coupled the laser to the structure, along with the actual accelerator stage.The resulting designs are   unintuitive enough that it is not clear a human could manually design a device of the same (or smaller) size that had at least the same efficiency.The accelerator from the final device of [SYV + 20] is shown in figure 4.
Lenses.There have been a few other types of devices that have been designed by these methods and created in practice.Some examples include flat lenses with a large depth of field that have the same focusing efficiency as traditional lenses [BPC + 20], and flat lenses with an ultra-wide field of view [BMM + 20].There has been some additional work in creating deformable lenses whose focal length can be controlled by stretching or contracting the material, with performances exceeding those of traditional lenses [CVJ + 18].
Fabrication constraints.We note that, while there have been many numerical demonstrations of inverse designed photonic devices, the actual fabrication of such devices has been relatively difficult until recently, when methods that could include fabrication constraints [PPSV17, VSS + 19] and could accommodate a large enough number of field variables became available and were fast enough to be used in practice.Additionally, fabrication constraints can also be applied to mass-produced photonic circuits in a foundry [PSS + 19], which makes it possible to fabricate inverse-designed photonics at scale.

Future directions
While physical design and, more specifically, inverse design, has led to some dramatic improvements in the performance of photonic devices, there are still many questions left unanswered and possible future avenues for research.
Standard benchmarks.There is a need for standardized benchmarks for both heuristic algorithms and bounds.Similar to the standard performance benchmarks for machine learning (such as MNIST [LBBH98], ImageNet [DDS + 09], CIFAR-10 [Kri09], etc.) and for optimization (such as the Maros-Mészáros test set [MM99], MIPLIB 2017 [GHG + ], etc.), which are used to compare the performance of different proposed algorithms, it is feasible to have a standard library of design specifications (objective functions and constraints) which either have known global solutions or a best known solution and best known lower bound that is updated as better ones are found.This would allow for researchers to have a concrete library which can be used to compare proposed algorithms against existing ones, in terms of both objective value and computational performance.Like in machine learning or optimization, we suspect that such benchmarks would clarify the performance tradeoffs of different algorithms and bounds.These benchmarks would help identify what problems are 'hard' for current heuristics and possibly lead to better approaches.Improved heuristics.While some of the methods presented here can scale to very large designs with millions to tens of millions of field and tens of thousands of design variables, many of the current methods do not make use of the specific structure of the physical design problem and take a long time (days or weeks) even with large computers or clusters.For example, one possible avenue is the use of a primal-dual algorithm, which can lead to drastically improved performance for some problems (see, e.g., [OV14, §3.5]).Such methods could also lead to algorithms whose natural byproduct of computing an approximately optimal design is a bound, like in §2, guaranteeing that the resulting design is close to the global optimum without requiring additional computational time.There has been additional work with combining machine learning approaches to speeding up both simulations [TSL + 19] and optimization [LTKY18, JF19, SBN + 20, CG20], which trade off training time for runtime performance.
Improved initializations.All of the heuristics used in practice often require good initializations in order to have reasonable performance.In practice, initializations that are guided by intuition appear to work well for many scenarios, but sometimes fail to reach what are known to be better designs.(See, e.g., [SVS + 20, §5].)As previously discussed, the bounds presented yield good initial designs as a by-product of the optimization procedure; but this need not be true in general, and we suspect that there might exist better methods which, while expensive to run exactly, might yield good initial designs.For example, while sign-flip descent might be potentially very expensive to run for large designs (as it requires solving a general constrained convex problem at each iteration) one could perform a small number of iterations to get a reasonably good initial field and then feed this resulting field into an optimization algorithm with faster convergence.(This latter procedure is sometimes called 'solution polishing' in combinatorial optimization.)Another avenue of similar work is the idea of 'objective-first' optimization [LV12] where the physics equation is relaxed to a penalty with a small weight and the resulting nonconvex problem is solved to get an infeasible design with good objective performance and (hopefully) small physics residual.This initial, infeasible design is then passed to any of the heuristics above to attempt to find a feasible design and field with similar performance.
Improved bounds.The small structures for which we have calculated bounds are simple cases where the bounds and the heuristics are very close; but this need not be the case.For example, the dual bounds given here can give weak performance limits when the objective function depends only on a few entries of the field variable z, even with the additional extensions presented.While we, in general, suspect the power bound yields tighter results that the diagonal dual lower bound in these cases, current off-the-shelf solvers are too slow to solve the resulting problem other than for a small number of field variables, n 10 3 .One future research avenue is to create faster solvers for problem (19), by exploiting the specific structure available in these problems.A second possibility is to find some connection between the current dual and power bounds; in particular, it is not clear how one is related to the other, if at all.We conjecture that the current power bounds (19) are always guaranteed to be at least as tight as the diagonal dual ones (23) whenever the function f is a separable convex quadratic, and where the constraint sets S k are singletons (equation ( 21) with S k = {k} for k = 1, . . ., n) but have been unable to prove this.It is also of practical interest to create lower bounds which give approximate scaling rules for designs, which would help a designer choose appropriate device sizes and materials for a desired objective.

Conclusion
There has been tremendous progress in the area of photonic inverse design, including foundrybased fabrication and the use of design software in a wide range of academic and industry labs.Still, there is a lot to be improved, including, but not limited to: designing standard benchmarks and comparisons for heuristics, overcoming computational bottlenecks in order to design larger devices than those which are currently possible, and the improvement of bounds (both theoretically and in terms of computational performance) which currently either do not apply, or are very difficult to compute, for many of the devices used in practice.Yet, despite the drawbacks of the current approaches, inverse design has yielded incredible practical benefits, and we expect that such methods will, in the near future, become a standard approach to creating practical, efficient devices that far exceed the performance of their traditionally-designed counterparts.

A Bound equivalence
We show that the bounds presented in [MCR20b,KM20] are equivalent to the ones given in section 2.2.
Derivation.First, we derive the bounds presented using this notation.We define the displacement field y i = γ i z i , for i = 1, . . ., n, where The diagonal physics equation can then be written in terms of the field z and the displacement field y: where γ ∈ [0, 1] n .We will define G = (2A 0 − I) −1 (this can be recognized, roughly speaking, as the Green's function) and b = 2Gb 0 (which is sometimes called the 'incident field') for notational convenience.We can then rewrite the diagonal physics equation using G and b : where g T i denotes the ith row of G. Multiplying on the left by y i , we find that y and z must satisfy y i (z i + g T i y) = y i b i , i = 1, . . ., n.Note that y 2 i = γ 2 i z 2 i ≤ γ i z 2 i = y i z i because γ 2 i ≤ γ i , since γ i ∈ [0, 1].Using this result, we find the following quadratic inequalities depending only on the displacement field y: y i (y i + g T i y) ≤ y i b i , i = 1, . . ., n.
Scaling each of the i inequalities by any nonnegative value λ i ≥ 0 and summing them together implies that y must satisfy: y T Dy + y T DGy ≤ y T Db , where D = diag(λ) is any matrix with nonnegative diagonal entries.These are the bounds presented in [MCR20b,KM20] in the case where Θ is a box rather than a Boolean constraint.
Tightness.As expected, these bounds are also tight in the sense that, if y satisfies (26) for all nonnegative diagonal matrices D, then there exists a feasible design θ and a field z such that z and θ satisfy the diagonal physical equation, and the displacement field y satisfies y i = (1 + θ i )z i /2 for −1 ≤ θ i ≤ 1, or, equivalently, that y i = γ i z i , where 0 ≤ γ i ≤ 1 for i = 1, . . ., n.
To show this, we will consider (as before) only the diagonal matrices D = e i e T i .The bound then implies that y i (y i + g T i y) ≤ y i b i , i = 1, . . ., n.We can then choose z i = b i − g T i y and γ i = y i /(b i − g T i y) b i − g T i y = 0 0 otherwise, for i = 1, . . ., n.Note that z i obviously satisfies z i + g T i y = b i and whenever b i − g T i y = 0, we have y i = γ i z i , by construction.On the other hand, if b i − g T i y = 0, then y 2 i = y i (y i + g T i y − b i ) ≤ 0, so y i = 0 = γ i z i , since γ i = 0 whenever b i − g T i y = 0, by definition.In other words, given a displacement field y satisfying (26) for all nonnegative diagonal matrices, we have constructed a field z and a design θ = 2γ − 1 such that y, z, and −1 ≤ θ ≤ 1 all satisfy the physics equation (25), or, equivalently, z and θ satisfy the diagonal physics equation.
Equivalence.The equivalence between the bounds follows immediately from the fact that (26) holds for all nonnegative diagonal matrices D if, and only if, the physics equation (25) holds for this choice of y, z, and θ, which, in turn, holds if, and only if, the original family of power bounds over z and θ holds for all diagonal matrices.

B Numerical examples
In this section, we focus on two basic numerical experiments for comparing both heuristics and lower bounds.Additionally, all Julia [BEKS17] code for the examples can be found at https://github.com/cvxgrp/pd-heuristics-and-bounds.The optimization problems used in computing the bounds were specified using the JuMP [DHL17] domain-specific language in Julia, and solved using Mosek [ApS19].

B.1 Small design
In this scenario, we are given a device whose field must satisfy the discretized Helmholtz equation in one dimension.At every point in the domain, the designer is allowed to choose the speed of the wave in the material within some specified range, such that the resulting field best matches some desired field.
Discretization.We can write a simple discretization of the above equation by subdividing the interval [−1, 1] into n equidistant points {x i } for i = 1, . . ., n.We then approximate the second derivative using finite-differences such that where h is the separation x i − x i−1 for any i = 2, . . ., n and is equal to h = 2/(n − 1), while z i = u(x i ) for i = 1, . . ., n.Additionally, we have defined z 0 = z n+1 = 0 for convenience.If we then define ω 2 /c(x i ) 2 = θ i , then we can approximate equation ( 27) as a diagonal physics equation of the form (5), where Problem data.In this problem, we will seek to minimize the squared Euclidean distance between the field z and some desired field ẑ; i.e., the objective is and we choose ẑ to be a cosine wave with a Gaussian envelope of width σ/ √ 2, whenever x i < 0 and 0 when x i ≥ 0: ẑi = cos(ωx i ) exp(−x 2 i /σ 2 ) x i < 0 0 x i ≥ 0.
Note that this function is discontinuous at 0. For the small numerical experiment, we will choose σ = 1/2, ω = 6π, and Θ = [1, 1.5] n .We will also set b (n+1)/2 = 2 and zero otherwise as our excitation, and set n, the number of gridpoints, to be n = 1001.This means that there are 1001 design parameters and field variables, which, in practice, is a very small number when compared to current applications.

B.2 Larger design
In this example, we will show an example of a larger physical design problem with n = d = 251 2 = 63001, that is similar in spirit to the smaller design, but is large enough that sparsity has to be exploited in order to have reasonable run time performance.

Figure 1 :
Figure 1: Approximate designs and desired field for 1D problem.

Figure 2 :
Figure 2: Approximate designs and desired field for 2D problem.

Figure 3 :
Figure 3: Scanning-electron microscope image of a three-way wavelength splitter with simulated fields overlaid.Figure and design from [SPS + 18].

Figure 4 :
Figure 4: Scanning-electron microscope image of an inverse-designed laser-driven particle accelerator with simulated fields overlaid.Figure and design from [SYV + 20].
Figure 4: Scanning-electron microscope image of an inverse-designed laser-driven particle accelerator with simulated fields overlaid.Figure and design from [SYV + 20].
Thus the Lagrange dual function gives us a performance bound, parametrized by λ. (Depending on the problem and choice of λ, it can give the trivial lower bound −∞.) Concavity.Since g(λ) is a performance bound for any λ ≥ 0, it is then natural to ask, what is the best possible performance bound?In other words, what is the largest possible value of g(λ) over the possible values of λ?This problem is called the dual problem and can be written as maximize g(λ) subject to λ ≥ 0. (

Table 1 :
Performance results for small design.

Table 2 :
Performance results for larger design.