Derivatives at frequency point in inverse optimization

[berzoll]

Hello,

I am trying to design a component using inverse optimization, similar to the splitter described here.

I would like to calculate the derivative of e.g. the power about the frequency point I am solving, but the objects inside the "J" function are already evaluated at the specified frequency point only. I tried passing a pair of frequencies to the optimization problem and calculating the derivatives myself, but the adjoint run takes increasingly longer the closer the two frequencies are spaced together.

Is there a way to access the full frequency spectrum at the monitors with a single run (as in normal meep simulations), or what would be the best way to approach this?

Thank you

[smartalecH]

but the adjoint run takes increasingly longer the closer the two frequencies are spaced together.

Unfortunately, this is a consequence of the hybrid time-/frequency-domain adjoint formulation. It's hard to overcome the Fourier uncertainty principle. (although, in theory, you could probably come up with a better fitting algorithm that doesn't require indivdual basis function for each frequency point -- it just turns into a nonconvex problem.)

I think the "right way" to approach this problem is to derive the equivalent adjoint source for the derivative of the power at a particular frequency point (i.e. don't try to approximate it with finite differences in frequency space). I have quite a bit of material on how to do this in chapter 2 of my dissertation (shameless plug).

There could be some "gotchas" with a frequency derivative. But you might get lucky. For example, the group velocity (which is a frequency derivative) can be computed just from field patterns thanks to the Hellman-Feynman theorem. Maybe something similar can happen here? @stevengj any tricks that you're aware of?

[stevengj]

First, if you want to go the finite-difference route and want to make it more accurate, I wouldn't space the frequencies closer together — instead, use a higher-order finite-difference formula. (Wikipedia has a table up to 8th order). Higher-order methods should work well because the frequency dependence is a smooth (analytic) function.

If you want to do this at lots of ω points, I would ditch finite differences and use Chebyshev interpolation, e.g. via the numpy.chebyshev routines — choose your frequencies at Chebyshev points, construct the polynomial, and then take whatever derivatives you want at whatever points you want. This should converge exponentially fast with the number of ω points, so you shouldn't need points that are extremely closely spaced. (Hopefully someone has implemented AD on the numpy Chebyshev routines, I haven't checked.)

If you want to do this analytically, you need to think of Meep as a frequency-domain solver, solving M(ω)e=j, where M is the frequency-domain Maxwell operator (given by eq. 32 in our paper), e is the frequency domain fields, and j is the frequency domain sources. To take d/dω of the power, you need de/dω = d/dω [M(ω)⁻¹] j. As explained in our matrix-calculus course, this is de/dω = – M(ω)⁻¹ dM/dω M(ω)⁻¹ j = M(ω)⁻¹ (–dM/dω e). So, you can see that it will involve two Maxwell solves (two Meep runs): first you do the ordinary run to solve for the fields e, then you multiply by –dM/dω (which, from equation (32), essentially means multiplying E and H by i and using them as current sources), then doing another Maxwell solve. Then, of course, to get the derivatives with respect to the design you will need to backpropagate through this, which is another two solves. It's straightforward in principle — and, crucially, can all be implemented in Python, with no fiddling with C++ internals, but will be a bit of a headache to get all the pieces working correctly.

[berzoll]

Thank you both.
Also, no shame for the plug indeed! That's amazing work!

An issue I see is that I am trying to compute derivatives in the vicinity of a bandgap. I guess this means that I cannot use any central finite-difference formula or the Chebyshev polynomial (although maybe using the forward finite difference could work?).

Is there any example on how one would proceed to calculate the group-velocity with the Hellman-Feynman theorem?
Will this allow me to retrieve information of the velocity/delay in between the source and the monitors, i.e. where the design region is, or rather just at the monitor itself?

I am not sure whether the adjoint solver can be e.g. group-velocity driven, rather than power/intensity driven. I tried making it optimize for some phase relationship between different frequency components, but had no success. I found inverse-design algorithms that are meant for dispersion engineering, but they seem to use different techniques.

[stevengj]

An issue I see is that I am trying to compute derivatives in the vicinity of a bandgap. I guess this means that I cannot use any central finite-difference formula or the Chebyshev polynomial (although maybe using the forward finite difference could work?).

I guess you are referring to van Hove singularities? In a finite-size simulation this won't be an issue — singularities only arise in an infinite system.