evaluating the order of accuracy of the adjoint gradient

[oskooi]

While we have already demonstrated second-order accuracy for subpixel smoothing of the MaterialGrid in a "forward" simulation (see results in #1503), following #1780 which enabled support for backpropagating the gradients from the Yee grid to the MaterialGrid, we have yet to validate the second-order accuracy of the gradients computed using an adjoint simulation.

ToDos

• using subpixel smoothing, verify that the gradient via the directional derivative converges with second-order accuracy with the Meep grid resolution
• using no smoothing, verify that the directional directive converges with first-order accuracy

The following 2d test computes the directional derivative using the adjoint gradient over a design region consisting of a ring structure which connects an input and output port/waveguide. (This is a similar setup to the existing unit tests for the adjoint solver in which the design region contains random pixel values and the directional derivative is computed in a random direction.) The objective function is a single DFT field point in the output waveguide at a single frequency (the red dot in the schematic below). (We would also eventually like to adapt this test to the other two objective functions: mode coefficients and near-to-far fields.) The ring structure is a MaterialGrid at a fixed resolution of 640 pixels/μm while the Meep grid resolution is repeatedly doubled: 20, 40, 80, and 160 pixels/μm. The directional derivative should be converging to a fixed value with increasing resolution but this does not appear to be the case. The output from running the script below shows the resolution in the second column and directional derivative in the third column.

grad:, 20.0, -8.3007345310497
grad:, 40.0, -9.140322608657096
grad:, 80.0, -8.32203124442412
grad:, 160.0, -7.0622182341423985

Before we can validate the accuracy of the gradients, we need to ensure that the results are converging to a fixed value.

import meep as mp
import meep.adjoint as mpa
import numpy as np
from autograd import numpy as npa

silicon = mp.Medium(epsilon=12)

sxy = 5.0
cell_size = mp.Vector3(sxy,sxy,0)

dpml = 1.0
pml_xy = [mp.PML(thickness=dpml)]

eig_parity = mp.EVEN_Y + mp.ODD_Z

design_region_size = mp.Vector3(1.5,1.5)
design_region_resolution = 640
Nx = int(design_region_resolution*design_region_size.x) + 1
Ny = int(design_region_resolution*design_region_size.y) + 1

w = 1.0
waveguide_geometry = [mp.Block(material=silicon,
                               center=mp.Vector3(),
                               size=mp.Vector3(mp.inf,w,mp.inf))]

fcen = 1/1.55
df = 0.23*fcen
wvg_source = [mp.EigenModeSource(src=mp.GaussianSource(fcen,fwidth=df),
                                 center=mp.Vector3(-0.5*sxy+dpml,0),
                                 size=mp.Vector3(0,sxy-2*dpml),
                                 eig_parity=eig_parity)]

x = np.linspace(-0.5*design_region_size.x,0.5*design_region_size.x,Nx)
y = np.linspace(-0.5*design_region_size.y,0.5*design_region_size.y,Ny)
xv, yv = np.meshgrid(x,y)

rad = 0.538948295
wdt = 0.194838432
ring_weights = np.logical_and(np.sqrt(np.square(xv) + np.square(yv)) > rad,
                              np.sqrt(np.square(xv) + np.square(yv)) < rad+wdt,
                              dtype=np.double)
ring_weights = ring_weights.flatten()

def adjoint_solver(resolution):
    matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
                              mp.air,
                              silicon,
                              weights=np.ones((Nx,Ny)))

    matgrid_region = mpa.DesignRegion(matgrid,
                                      volume=mp.Volume(center=mp.Vector3(),
                                                       size=mp.Vector3(design_region_size.x,
                                                                       design_region_size.y,
                                                                       0)))

    matgrid_geometry = [mp.Block(center=matgrid_region.center,
                                 size=matgrid_region.size,
                                 material=matgrid)]

    geometry = waveguide_geometry + matgrid_geometry

    sim = mp.Simulation(resolution=resolution,
                        cell_size=cell_size,
                        boundary_layers=pml_xy,
                        sources=wvg_source,
                        geometry=geometry)

    obj_list = [mpa.FourierFields(sim,
                                  mp.Volume(center=mp.Vector3(1.25),
                                            size=mp.Vector3()),
                                  mp.Ez)]

    def J(mode_mon):
        return npa.power(npa.abs(mode_mon),2)

    opt = mpa.OptimizationProblem(
        simulation=sim,
        objective_functions=J,
        objective_arguments=obj_list,
        design_regions=[matgrid_region],
        frequencies=[fcen])

    f, dJ_du = opt([ring_weights])

    sim.reset_meep()

    return f, dJ_du



for r in [20., 40., 80., 160.]:
    adjsol_obj, adjsol_grad = adjoint_solver(r)
    print("grad:, {}, {}".format(r,np.sum(adjsol_grad)))

[smartalecH]

Even though smoothing is "on", it's not going to work unless the input is smoothed and beta=np.inf.

[ianwilliamson]

It looks like your design resolution (set to 640?) is much higher than even your highest simulation resolution (160). Is this intended? I thought the correct approach was to use a comparable or lower design resolution to ensure a slowly varying latent design over each Yee grid cell?

One issue with this approach is that the bilinear interpolation from a high resolution design array to a relatively lower resolution Yee grid will lead to a very sparse gradient, where may of the gradient values are zero (by definition of the bilinear interpolation). The locations of the pixels which are non-zero will probably shift significantly with changes to the Yee grid resolution, which could partially explain the behavior you're seeing. It would be useful to plot the gradient field to check this.

[oskooi]

It looks like your design resolution (set to 640?) is much higher than even your highest simulation resolution (160). Is this intended?

This is intentional because we want to ensure that, as part of the convergence study, the MaterialGrid is always being down-sampled to Meep's Yee grid. Note that in this example, the inner ring radius is ~0.5 μm and the ring width is ~0.2 μm while the largest Meep grid voxel dimension is 0.05 μm (at the smallest resolution of 20 pixels/μm). This means the MaterialGrid geometry should be resolvable on even the coarsest Yee grid used in the convergence study.

The locations of the pixels which are non-zero will probably shift significantly with changes to the Yee grid resolution, which could partially explain the behavior you're seeing.

The gradient for this geometry is indeed sparse (only the values at the dielectric interfaces are non-zero) which could therefore cause large changes in the directional derivative with "small" changes in the Meep grid resolution. Will post results on this topic shortly.

@smartalecH I did try smoothing the ring geometry of the MaterialGrid and using the built-in projection feature with beta=np.inf. However, the directional derivative still does not appear to be converging to a fixed value with increasing Meep grid resolution.

[smartalecH]

I thought the correct approach was to use a comparable or lower design resolution to ensure a slowly varying latent design over each Yee grid cell?

What really matters is that the first-order Taylor approximation used to define the interfaces of the new subpixel smoothing routine is satisfied. The intended way to do this is using a smoothing filter (in fact, @stevengj derived this approach with the assumption the input variables would be smoothed, as the original adjoint solver was built on a density-based framework).

One issue with this approach is that the bilinear interpolation from a high-resolution design array to a relatively lower resolution Yee grid will lead to a very sparse gradient

All the more reason to first smooth using a filter. Backpropagating through this filter is just a transposed convolution, which will "smear" the gradient energy to all surrounding pixels so that the resulting gradient is no longer sparse. The filter will also "regularize" the problem (although that doesn't matter here).

The gradient for this geometry is indeed sparse (only the values at the dielectric interfaces are non-zero) which could therefore cause large changes in the directional derivative with "small" changes in the Meep grid resolution. Will post results on this topic shortly.

This is actually a different kind of sparsity than what @ianwilliamson is referring to. This kind of sparsity is a result of levelset "shape optimization" methods (i.e. when beta=np.inf), where only the boundaries contain gradient information.

I did try smoothing the ring geometry of the MaterialGrid and using the built-in projection feature with beta=np.inf. However, the directional derivative still does not appear to be converging to a fixed value with increasing Meep grid resolution.

@oskooi can you post the code and the results? I don't see any filtering or projection in the above script.

How are you computing the directional derivative? The above code is just computing the sum of the gradient. You're not scaling with resolution (i.e. doing a proper integral). In that case, we shouldn't expect the sum to converge to a finite value with increasing resolution (although the integral should converge...I think...). Edit: although if the design region resolution is consistent, the sum should also work.

But the directional derivative should converge, assuming you're acting in the same direction with increasing resolution (which means your design grid resolution must remain the same throughout all runs).

[ianwilliamson]

This is intentional because we want to ensure that, as part of the convergence study, the MaterialGrid is always being down-sampled to Meep's Yee grid.

This is what I was trying to point out. By definition, using bilinear interpolation to downsample from high resolution to low resolution will always create "sparsity" in the gradient. This is the correct gradient for the operation being performed, but I'm not sure if it's what you want for the test you're creating.

What really matters is that the first-order Taylor approximation used to define the interfaces of the new subpixel smoothing routine is satisfied. The intended way to do this is using a smoothing filter (in fact, @stevengj derived this approach with the assumption the input variables would be smoothed, as the original adjoint solver was built on a density-based framework).

I think even with a filter applied, a sufficiently high resolution design array will still have sparsity in the gradient. It should depend on the ratio of the resolutions and the size of the filter kernel. If your filter doesn't cover enough pixels in the high resolution design array, then some design pixels will just never influence the simulation and, as a result, will have zero gradient. Just want to point out that having a filter doesn't guarantee the intended behavior.

[oskooi]

Here are some results for the adjoint gradient map obtained using subpixel smoothing (script below) for a fixed MaterialGrid resolution (160 pixels/μm) and different combinations of (1) Meep resolution (40 and 80 pixels/μm), (2) the tanh projection bias β (1.0 and 10.0), and (3) the radius of the conic filter (16 and 64 pixels).

Some comments:

• increasing the filter radius does seem to mitigate the sparsity even when β is large (compare each inset in the two figures below)
• the larger the down sampling (i.e., the ratio of the MaterialGrid resolution to the Meep grid resolution), the greater the sparsity
• increasing β increases the overall sparsity of the gradient map by forcing the gradient to be non-zero mainly at or near the dielectric interfaces

The next thing to check based on the first comment is to redo the convergence study with a larger filter radius.

1. filter radius: 16 pixels

2. filter radius: 64 pixels

The script used for subpixel smoothing is as follows:

weights = np.logical_and(np.sqrt(np.square(xv) + np.square(yv)) > rad,
                         np.sqrt(np.square(xv) + np.square(yv)) < rad+wdt,
                         dtype=np.double)

filter_radius = 16/design_region_resolution
filtered_weights = mpa.conic_filter(weights,
                                    filter_radius,
                                    design_region_size.x,
                                    design_region_size.y,
                                    design_region_resolution)

ring_weights = filtered_weights.flatten()

def adjoint_solver(resolution):
    matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
                              mp.air,
                              silicon,
                              weights=np.ones((Nx,Ny)),
                              do_averaging=True,
                              beta=np.inf)

    f, dJ_du = opt([ring_weights])

[oskooi]

Quick update: increasing the filter radius from 16 pixels to 32 and then 64 produces results for the directional derivative which appear to be converging to a fixed value. However, as shown below, the convergence rate appears to be linear (or slightly greater) but not quadratic. It's likely that the convergence rate is indeed quadratic and we just need to extend the data out to larger resolutions to see the trend more clearly...

[smartalecH]

Are you backpropagating the gradient through the filter? How are you computing the directional derivative?

[oskooi]

Based on recent discussions, I modified the test above to remove the input and output waveguides (since the terminations contain sharp corners and thus introduce field discontinuities which can affect the results) so there is now just the ring structure. I extended the Meep grid resolution range to 20, 40, 80, 160, and 320 using a MaterialGrid with a fixed resolution of 1280. The "exact" result is computed at a Meep grid resolution of 640. The gradient computed on the MaterialGrid is backpropagated through the conic filter using the vJp. The directional derivative is the sum of all the values of the gradient map (equivalent to a direction vector of [1, 1, 1, ... ]). The full script is here.

The figure below shows the results for three different filter radii. The convergence rate is (nearly) second order. Note that increasing the filter radius does not seem to significantly change the results. For completeness, we can try to continue extending the resolution range even further but these results seem sufficient.

Next steps are (1) to demonstrate that reducing the filter radius to zero will eventually lead to a first-order convergence rate and (2) show results for no smoothing.

[oskooi]

As expected, gradually decreasing the filter radius to zero seems to be causing the convergence rate to switch to linear from quadratic for large filter radii. This is shown in the figure below for the results of the smallest filter radius of 2 pixels.

[stevengj]

Maybe a clearer way to present this would be to plot the |gradient - FD gradient| vs resolution, showing that gradient always matches the forward solver.

Ideally, this error should be small regardless of the resolution. That way, if the forward solver is converging quadratically, the gradient necessarily is as well.