New adjoint filtering scheme

[smartalecH]

The current filtering scheme often induces a spatial delay in each dimension. This implies the filter kernel is not a "zero-phase" filter. In theory, all of the example impulse responses currently in the codebase (top-hat, conic, and Gaussian) should be zero-phase, indicating our setup is wrong.

One way to get around this is to implement a "filt-filt" (essentially filter in both directions, and remove the phase delay). This is what is implemented in #1625. This, of course, is not what we want, as this changes the underlying kernel of the final filter operation, rendering all of our minimum length-scale constraint math invalid. It should be possible to properly implement a single-pass filter with the above kernels.

Since all of our kernels are separable, we can actually do a batched convolution (natively supported by autograd) so long as we construct our filter kernels properly and apply the right padding (to preserve the correct shape).

Right now we do an fft-based filter. We may not need to with a batched convolution (depending on the size of the kernel, an fft convolution may be slower).

(cc @hammy4815)

[stevengj]

This coordinate system looks wrong to me. For FFT-based filtering one has periodic boundary conditions, so ±L/2 are equivalent points and should not be included twice. (Then you subsequently zero-pad, but you do so by equal amounts on both sides, so you have the same problem.)

It's pretty easy to get the boundary conditions / origin slightly wrong for this sort of thing, especially in conjunction with padding and fftshift.

[stevengj]

One option would be to use np.fftfreq instead of linspace to get the coordinates, then zero-pad in the center and omit the fftshift. This might be less error-prone.

[stevengj]

One way to get around this is to implement a "filt-filt" (essentially filter in both directions, and remove the phase delay). This is what is implemented in #1625. This, of course, is not what we want, as this changes the underlying kernel of the final filter operation, rendering all of our minimum length-scale constraint math invalid

Wouldn't a simpler option just be to symmetrize the filter by adding it to its mirror flip (taking care that mirror flip is not equivalent to np.flip for FFT symmetry)? At most this changes the filter by one (material-grid) pixel.

[smartalecH]

For FFT-based filtering one has periodic boundary conditions, so ±L/2 are equivalent points and should not be included twice.

The kernel is symmetric though. So f(-L/2) = f(L/2).

One option would be to use np.fftfreq instead of linspace to get the coordinates, then zero-pad in the center and omit the fftshift. This might be less error-prone.

Note that we are first building the kernel in the spatial domain, and then running an FFT on that. We are only using the coordinates to build the kernel itself. Would fftfreq still help in this case?

(Then you subsequently zero-pad, but you do so by equal amounts on both sides, so you have the same problem.)

But if the padding is identical on both sides, it's still "periodic", right?

As you say, however, I'm sure there's a bug somewhere w.r.t. the boundary conditions or origin.

[stevengj]

The kernel is symmetric though. So f(-L/2) = f(L/2)

The point is that if the period were L (no zero padding), then you would only store f(L/2), not f(-L/2). More generally, the key point is to identify the correct "origin" for the coordinate system.

Natively in an DFT, a "symmetric" convolution kernel of length n means that x[i mod n] == x[(n-i) mod n], so that the "center" is x[0].

If n is even, then using ifftshift changes the origin (before the shift) to the index n/2. Note that this is not identical padding on both sides: there are n/2 elements before the origin, and n/2-1 elements after the origin.

[smartalecH]

If n is even, then using ifftshift changes the origin (before the shift) to the index n/2. Note that this is not identical padding on both sides: there are n/2 elements before the origin, and n/2-1 elements after the origin.

Ah this is the issue. Turns out ifftshift has logic to check if n is even or not, and does the proper shift (as you mention). The normal fftshift function does not have this logic. Currently, the codebase uses fftshift:

meep/python/adjoint/filters.py

Line 139 in 34b62e9

y = npa.fft.fftshift(npa.real(npa.fft.ifft2(Y)))

From their docs:

The inverse of fftshift. Although identical for even-length x, the functions differ by one sample for odd-length x.

[stevengj]

In particular, for a filter length n (even or odd), I think you want the filter to be centered at n/2 rounded down before the ifftshift. In particular:

1. Create the filter kernel array of length n centered at n/2 (rounded down).
2. Zero-pad it as needed by adding k zeros to both ends. (Now the center is at n/2 + k == (n+2k)/2 rounded down, which is correct for the new length n+2k.)
3. ifftshift the filter array.
4. Do the convolution: ifft(fft(ρ) * fft(paddedkernel))

[smartalecH]

Zero-pad it as needed by adding k zeros to both ends. (Now the center is at n/2 + k == (n+2k)/2 rounded down, which is correct for the new length n+2k.)

By zero-padding, doesn't that mean the new array is no longer symmetric (if the number of taps is even)?

For example, suppose x = [1,2,3,2]. x is indeed symmetric (its DFT is purely real):

> fft(ifftshift([1,2,3,2]))
[8.+0.j 2.+0.j 0.+0.j 2.+0.j]

However, as soon as I pad it (y=[0,0,0,0,1,2,3,2,0,0,0,0]) It's no longer symmetric (with or without the ifftshift):

> fft.fft(fft.ifftshift([0,0,0,0,1, 3, 2, 3,0,0,0,0]))
[ 9.        +0.00000000e+00j  7.69615242+8.66025404e-01j
  4.5       +8.66025404e-01j  1.        +0.00000000e+00j
 -1.5       -8.66025404e-01j -2.69615242-8.66025404e-01j
 -3.        -2.22044605e-16j -2.69615242+8.66025404e-01j
 -1.5       +8.66025404e-01j  1.        +0.00000000e+00j
  4.5       -8.66025404e-01j  7.69615242-8.66025404e-01j]

Or am I missing something?

[stevengj]

Yes, right, sorry. You should split the Nyquist (middle) element 50–50 between the positive and negative components when zero-padding an even-length array.

i.e. zero-pad with something like:

# zero pad by k (before ifftshift)
def symmetric_zeropad(x, k):
    y = np.pad(x, (k, k))
    if len(x) % 2 == 0: # even
        middle = x[0]
        y[k] = middle/2
        y[-k] = middle/2
    return y

[stevengj]

e.g.

>>> x = [1,2,3.0,2]
>>> y = symmetric_zeropad(x, 3)
array([0. , 0. , 0. , 0.5, 2. , 3. , 2. , 0.5, 0. , 0. ])
>>> np.fft.ifftshift(y)
array([3. , 2. , 0.5, 0. , 0. , 0. , 0. , 0. , 0.5, 2. ])
>>> np.fft.fft(np.fft.ifftshift(y))
array([8.00000000e+00+0.j, 6.54508497e+00+0.j, 3.42705098e+00+0.j,
       9.54915028e-01+0.j, 7.29490169e-02+0.j, 3.33066907e-16+0.j,
       7.29490169e-02+0.j, 9.54915028e-01+0.j, 3.42705098e+00+0.j,
       6.54508497e+00+0.j])

[stevengj]

Alternatively, it might be clearer to enforce that the convolution kernel is symmetric by only specifying half of it (from [0,L/2]), and then constructing the other half by mirror flipping (and padding as needed).