odlgroup / odl

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

Integration operator #440

Open ozanoktem opened 8 years ago

ozanoktem commented 8 years ago

We need to implement an integration functional which accepts a real value discretized function and evaluates, using some quadrature rule, its integral. It should account for boundary behavior and support integration along axes in multidimensional functions. Such functionals are needed when defining many regularization functionals and currently we use the norm of the underlying space, which can lead to working results. One idea would be to use the scipy.integrate package. One may also check out the Python wrapper for the Cubature algorithm for multidimensional integration (cubature) of vector-valued functions over hypercubes.

adler-j commented 8 years ago

One obvious comment is that the integral needs to match the inner product and the likes, that is (on unweighted spaces):

x.norm()**2 == (a ** 2).integrate()
x.inner(y) == (x * y).integrate()

Another is that we will likely need the gradient of the integral (that is, the weights).

kohr-h commented 8 years ago

One obvious comment is that the integral needs to match the inner product and the likes

That actually depends. Perhaps you want to have the vanilla integral even if your reconstruction space is weighted? Or, in other words, sometimes you want to manually set the weighting in the integral functional.

adler-j commented 8 years ago

That actually depends. Perhaps you want to have the vanilla integral even if your reconstruction space is weighted? Or, in other words, sometimes you want to manually set the weighting in the integral functional.

Well, if the integral is a property of the vector, then it should match, if it is a freestanding functional, then we can ofc give options. Anyway, my comment was specifically on unweighted spaces. For weighted spaces I expect the result to be:

x.norm()**2 == (w * a ** 2).integrate()
x.inner(y) == (w * x * y).integrate()
kohr-h commented 8 years ago

One idea would be to use the scipy.integrate package.

That one is intended for integration of a continuously given function where you can freely select the quadrature points. Possible in principle, but arbitrarily slow.

One may also check out the Python wrapper for the Cubature algorithm for multidimensional integration (cubature) of vector-valued functions over hypercubes.

Same here, although it may be faster by a more sophisticated handling of the integrand function. But why would we want to depend on such a package?

kohr-h commented 8 years ago

We need to implement an integration functional which accepts a real value discretized function and evaluates, using some quadrature rule, its integral.

If it's already discretized (in the way we understand now), there is no choice of quadrature rule.

ozanoktem commented 8 years ago

I would not claim that integral is a property of the vector space because it is not. In fact, the notion of integral is not necessarily related to norm. Furthermore, I would not go as far as to claim that the quadrature rule is uniquely given by the discretization. At least not when the discretization lacks an interpolation operator.

adler-j commented 8 years ago

I would not claim that integral is a property of the vector space because it is not. In fact, the notion of integral is not necessarily related to norm.

While this is true in the most general setting, it depends on how you define your vector space. For Lp type spaces most authors define the norm by

(\int |x|^p d \mu)^(1/p)

Here the measure \mu plays an integral role in the definition and this role would probably be expected to follow through to other integrals (at least as the default value). This should (imo) be a guideline in how we chose defaults etc. This should only be complicated for points on the edge.

Furthermore, I would not go as far as to claim that the quadrature rule is uniquely given by the discretization. At least not when the discretization lacks an interpolation operator.

If no interpolation is given the choice is clearly open. However, we dont currently have any discretization where this is the case.

aringh commented 8 years ago

TODO: when this is implemented, the evaluation of the L1-norm in default_functionals should be changed to use this.