[Octave-bug-tracker] [bug #62468] quadgk.m should be able to do array-valued integration

[Michael Leitner]

URL:
  <https://savannah.gnu.org/bugs/?62468>

                 Summary: quadgk.m should be able to do array-valued
integration
                 Project: GNU Octave
            Submitted by: mleitner
            Submitted on: Sun 15 May 2022 05:03:41 PM UTC
                Category: Octave Function
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: None
                  Status: None
             Assigned to: None
         Originator Name: 
        Originator Email: 
             Open/Closed: Open
                 Release: dev
         Discussion Lock: Any
        Operating System: Any

    _______________________________________________________

Details:

As discussed over at bug #62412, it would be nice if quadgk could compute
vector-valued integrals. Note that in the context of this discussion it is
always understood that the quadrature algorithm evaluates the different
entries at exactly the same points, specifically it adaptively refines the
points for all entries in exactly the same way.

First why it is worthwhile to be able to directly compute vector-valued
integrals: Of course, if you have a function like
++
f=@(x) cos(x.*(1:1000))
--
then you could just as well compute the 1000 entries separately directly.
However, if f should be a black box, or if it is something like
++
f=@(x) cos(very_complicated_function(x).*(1:1000))
--
then it could be a factor 1000 more efficient to do the vector-valued
integral. Naturally, this is only the best case, as the different entries can
have numerically difficult behaviour in different parts of the integration
range, so that entry-wise integration would adaptively refine differently (and
thus more efficiently).

At present, quadv is the only routine that can do vector-valued integrals.
However, it does so by the Simpson rule, which is of quite low order and thus
for sufficiently smooth functions is expected to be less efficient than
algorithms of higher order, and further it is undesirable to have
unnecessarily different calling conventions (and returned parameters) for the
different routines, as well as that one is forced to use different quadrature
rules for scalar and array-valued integrations.

I have now given quadgk the ability to compute vector-valued integrals. For
now, it is only the dumbest solution, namely the function is evaluated always
only point-wise, which means that you have as many function calls (of the
interpreter) as (algorithmic) function evaluations, while in the scalar-valued
code path you will have a number of function calls that goes probably with the
logarithm of the requested precision. For an interpreted language like octave,
the former behaviour can cost significantly in efficiency, but this is also
the present state with quadv.

It would be conceivable to get better behaviour for cases like 
++
f=@(x) [cos(x); sin(x)]
--
where f([1 2]) gives a 2x2-result, where the second dimension (or in general,
the last dimension) varies of the evaluation points. However, this does not
work if the function uses anything beyond elementwise operations, such as
matrix operations -- one could try to call it first with a vector-valued
argument, and if that fails, go back to the present code path, but this does
not seem robust -- it is conceivable to have functions that work for
vector-valued inputs, but where it is not necessarily true that 
++
cat(n+1,f(1),f(2))==f(shiftdim([1; 2],-n))
--
where n is the number of dimensions of the output of f. 

It is not yet finished -- the documentation has to be updated, BISTs for the
new functionality are missing (they can be taken from quadv, but also there
there are only two), and I have not done much testing, specifically with
respect to whether it honors the requested precision. 

Finally: already in the previous version there was a code block in the
subfunction __quadgk_eval__ that sets too_close to true if the subinterval has
become very small, sets the contribution of the interval to zero, and does not
even emit a warning. The motivation is seemingly that such small intervals can
only happen due to numerical noise, and that this has to be suppressed. In my
view, a function in a numerical framework is just as definitive as a function
in abstract mathematics -- on a given hardware, it gives a given value, and if
it shows noise, then this noise belongs to the function, and should contribute
to the integral. Perhaps a warning should be emitted, with the specific point
in the integration range where this happens.




    _______________________________________________________
File Attachments:


-------------------------------------------------------
Date: Sun 15 May 2022 05:03:41 PM UTC  Name: quadgk.m  Size: 24KiB   By:
mleitner

<http://savannah.gnu.org/bugs/download.php?file_id=53205>

    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?62468>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/