odlgroup / odl

Operator Discretization Library https://odlgroup.github.io/odl/
Mozilla Public License 2.0
375 stars 105 forks source link

WaveletTransform adjoint is wrong #733

Open adler-j opened 8 years ago

adler-j commented 8 years ago

By adding this test:


def test_wavelet_adjoint(wave_impl, shape_setup, floating_dtype):
    """Verify that the adjoint of the wavelet operator is correct.

    This is done by checking the definition of the adjoint::

        <Ax, y> = <x, Aty>
    """
    # Verify that the operator works as expected
    wavelet, pad_mode, nlevels, shape, _ = shape_setup
    ndim = len(shape)

    space = odl.uniform_discr([-1] * ndim, [1] * ndim, shape,
                              dtype=floating_dtype)

    wave_trafo = odl.trafos.WaveletTransform(
        space, wavelet, nlevels, pad_mode, impl=wave_impl)

    Axy = wave_trafo(x).inner(y)
    xAty = x.inner(wave_trafo.adjoint(y))

    assert Axy == pytest.approx(xAty)

we find that the wavelet transform has a wildly wrong adjoint, off by a scaling factor.

By adding a scaling factor to the last lines, i.e

    scale = 1/wave_trafo.domain.partition.cell_volume
    xAty = scale * x.inner(wave_trafo.adjoint(y))

most tests start working, but certainly not all.

kohr-h commented 8 years ago

It's a known issue that the wavelet adjoint is wrong, at least for non-orthogonal wavelets. See #400, #412 and PR #421 for more info. That PR is incomplete, though, and needs some more work. It has overlap with boundary issues for real-space convolutions. I have the expressions written down but implementation is not done. Surprising though that this volume scaling isn't taken care of.

adler-j commented 8 years ago

Most surprising is that we didn't have tests for this, frankly. Also added tests for fourier transform now and they also fail. All of them.