bugGNU Octave - Bugs: bug #62468, quadgk.m should be able to do...


bug #62468: quadgk.m should be able to do array-valued integration

Submitter:  Michael Leitner <mleitner>
Submitted:  Sun 15 May 2022 05:03:41 PM UTC
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  None
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment (Rich Markup)

Comment Type & Canned Response:


( Jump to the original submission )

Thu 02 Jun 2022 11:59:26 PM UTC, comment #28: 

i pushed http://hg.savannah.gnu.org/hgweb/octave/rev/451fb63a10a0

updates integral to call quadgk instead of quadv. simplified some of the selection logic, updated docstring, BISTs and NEWS, and a typo in quadgk. passes all tests.

Nicholas Jankowski <nrjank>
Project Member
Thu 02 Jun 2022 04:32:33 PM UTC, comment #27: 

I can take a look at that if you’re not on it already

Nicholas Jankowski <nrjank>
Project Member
Thu 02 Jun 2022 04:30:33 PM UTC, comment #26: 

No further comments so I checked in the new quadgk here: http://hg.savannah.gnu.org/hgweb/octave/rev/12f8fb75fc30.

Marking as fixed and closing report.

Oh, I guess integral.m should be updated to use the new quadgk.m rather than quadv.m for "ArrayValued" integrations.

Rik <rik5>
Project Administrator
Thu 26 May 2022 03:02:47 PM UTC, comment #25: 

I'm uploading what I believe is a now complete quadgk.m file.  I implemented the "ArrayValued" input option and updated the documentation for both "ArrayValued" and for "Waypoints" which I thought could be clearer.  Lots of small fixups in the code such as using rows() and columns() instead of size (x,1) or size (x,2) for clarity.  There are a few BIST tests for the vectorized code path.

This is ready for final review and then it can be checked in.

(file #53257)

Rik <rik5>
Project Administrator
Wed 25 May 2022 04:02:19 PM UTC, comment #24: 

I like Nicholas's suggestion that "ArrayValued" be an option to pass to quadgk to explicitly inform the function to expect a vectorized integrand.

Rik <rik5>
Project Administrator
Fri 20 May 2022 06:11:51 PM UTC, comment #23: 

since this quadgk function is an extension beyond what matlab's quadgk does (they seem to have extended quadgk to use in integral, but not actually modified/extended their quadgk), I think it might be best to make use of the same "'arrayvalued',true" method of manually selecting the alternate path. This will prevent an incorrect detection either causing unnecessary performance degradation one way or function errors the other way. that same straightforward selection would also apply when we later modify integral to call it.

Nicholas Jankowski <nrjank>
Project Member
Fri 20 May 2022 08:21:57 AM UTC, comment #22: 

I completely agree with your comments. In the non-arrayvalued code paths, both octave and matlab require the integrands to be vectorized, meaning you can pass in arrays x of any size, they return a with size(a)==size(x), and for all i in [1..numel(x)] it is fulfilled f(x(i))==f(x)(i). This is to allow vectorized evaluations, which in these interpreted language are much more efficient than element-wise.

And yes, in my octave-implementation of the arrayvalued codepath in quadgk, as well as in octave's existing quadv, the price you pay for being able to integrate array-valued functions is that you can only evaluate them point-wise. And as it was to be expected, your test how matlab behaves implies that it does it also like that. But matlab has JIT, meaning that it probably doesn't hurt so much on their side.

The only difference is that my implementation of quadgk tries to decide on its own whether the integrand is array-valued or not. I think there should be no problems with that, but it is to be decided whether we want this behaviour, or whether we do it like matlab and require the user to explicitly specify array-valued integration if it should be done.

Michael Leitner <mleitner>
Thu 19 May 2022 08:45:34 PM UTC, comment #21: 

Okay, I looking back I think my size comment doesn't apply to 'arrayvalued'. My understanding of Matlab's integrator requiremnet for other functions is that for any f such that a = f(x), no matter what goes on inside f, size(a)=size(x).  that way it can iterate on f with vectorized expansion of x of any size, and not have 'a' grow or shrink in size. 

I guess this requirement is disabled for 'arrayvalued', where they probably chose a less efficient codepath to allow arbitrary resizing.  e.g.,

integral(@(x) x*x, 0, 1,'arrayvalued',1)

ans =


>> integral(@(x) x*x, 0, 1)
Error using  *
Incorrect dimensions for matrix multiplication. Check that the
number of columns in the first matrix matches the number of rows
in the second matrix. To operate on each element of the matrix
individually, use TIMES (.*) for elementwise multiplication.

Error in @(x)x*x

Nicholas Jankowski <nrjank>
Project Member
Thu 19 May 2022 07:44:25 PM UTC, comment #20: 

I don't think that the double loop is slowing it down so much. At least not when the dimension of the array returned by a single call to f is massive, as in your examples -- 150 calls that operate on 1e5 numbers each do not really take longer than one call that operates on 1.5e7 numbers. But yes, for more realistic cases it is not efficient, that's why I said it is only the dumbest solution. Please see my thoughts on possibilities to improve in the original post: if f uses internally only element-wise operations, it would be possible, but there is no way to find that out with certainty that I see. If f uses matrix arithmetic to arrive at its array-valued output and is written expecting only scalar inputs, then there is definitely no general way to do this vectorized. And no, that it is a double loop is not a problem, with some effort one could code that as a single loop, but one wouldn't gain anything.

And I don't really understand what you mean with size(out)==size(in): in the array-valued code path, out is definitely larger (it has even additional dimensions) than in.

Michael Leitner <mleitner>
Thu 19 May 2022 07:33:04 PM UTC, comment #19: 

ok, actually looking at quadgkv_eval_array and the two level loop that is probably slowing it down the most.  is it possible to redefine f(x) = f(x)(:) before passing to the subfunction such that f(t(:,i,j)) could be a single vectorized call? or at least one loop level done in a shot? Since there's a requirement that size(out)=size(in) for arbitrary vector size of in, it seems like that should be possible? but maybe I'm just not thinking of enough things that would break it.

Nicholas Jankowski <nrjank>
Project Member
Thu 19 May 2022 07:05:59 PM UTC, comment #18: 

yeah, I don't think we need to worry too much about optimizing performance further here, since like many other things speed has always been worse on this side, and it's not like quadv was much faster.  I just thought it worth checking that it wasn't an obvious problem with an algorithm choice in these changes.

from here, it looks like the function passes all of the quadv and integral arrayvalued tests I'd expect.

Other than the docstring, can anyone think of anything else it should do before rolling it back over to quadgk and updating integral accordingly?

Nicholas Jankowski <nrjank>
Project Member
Thu 19 May 2022 06:27:47 PM UTC, comment #17: 

Next version: It is now a bit faster (not much, and will only be visible for really trivial integrands) due to how the array in __quadgkv_eval_array__ is filled. But yes, it takes about 20 times as long for the execution of __quadgkv_eval_array__ as it would take for purely evaluating the function. This has to be due to the transformation that is applied to it.

(file #53233)

Michael Leitner <mleitner>
Thu 19 May 2022 06:00:04 PM UTC, comment #16: 

copy/paste error in last comment. intmax line with the error was obv not a 1x500000000 (37.3GB) array. should have said
"Requested 1x2147483646 (16.0GB) array"

Nicholas Jankowski <nrjank>
Project Member
Thu 19 May 2022 05:58:36 PM UTC, comment #15: 

regarding my previous question - so intmax ~ 2.1e9 . calling something silly like:

tic;a = integral(@(x) repelem(x,5e8), 1, 1,'arrayvalued',true);toc
Elapsed time is 5.079442 seconds.
>> tic;a = integral(@(x) repelem(x,intmax), 1, 1,'arrayvalued',true);toc
Error using repelem
Requested 1x5000000000 (37.3GB) array exceeds maximum array size
preference (15.9GB). This might cause MATLAB to become

Error in @(x)repelem(x,5e9)
>> tic;a = integral(@(x) repelem(x,intmax-5e8), 1, 1,'arrayvalued',true);toc
Elapsed time is 56.228911 seconds.

while in octave,

tic;quadgkv(@(x) repelems(x,[1;5000000]), 1, 1);toc
Elapsed time is 43.6698 seconds.

(profiler shows most time spent in _quadgkv_eval_array_)
didn't take it up any higher, and my system can just barely fit an intmax sized double array anyway.  I had just been wondering if octave's linear indexing limit would be a problem anywhere, but seems performance limits will hit first.

Nicholas Jankowski <nrjank>
Project Member
Thu 19 May 2022 05:48:43 PM UTC, comment #14: 

Yes, there is a (:) hidden inside. Which trivial integrand did get really slow? If it is really trivial, then already the first stage of the refining should return, where you evaluate the function at 650 points. And 650*5e7 is only 3e10. Multiplying this number of entries by weighting factors and summing should not take so extraordinarily long.

Michael Leitner <mleitner>
Thu 19 May 2022 05:41:32 PM UTC, comment #13: 

playing with test cases including ndim arrays, noticed that neither octave nor matlab's quadv can handle something like

quadv(@(x) cat(3,x,x), 0, 1)

but matlab's integral can and this quadgkv can. i suspect because of reshaping to a single column for processinng.  that is going to limit any large data processing to numel=intmax, no?  I don't know if that's something we hould be concerned about.

(I started testing and noticed that a trivial integrand got really slow for n>50M, so I didn't push it to intmax)

Nicholas Jankowski <nrjank>
Project Member
Thu 19 May 2022 05:17:13 PM UTC, comment #12: 

Fixed another small bug if the integrand is specified like

f = (x) 1;

so that it can always be called in a vectorized way. This was only relevant in the complex integration path case, as otherwise the transformation had already vectorized it.

(file #53232)

Michael Leitner <mleitner>
Thu 19 May 2022 05:01:09 PM UTC, comment #11: 

Why the recursive call: because otherwise you have to set a flag to multiply the result at the end. But yes, this seems better to me.

For testing, I attach now a version that incorporates all the points mentioned in the meantime:
- I deleted the (as yet inactive) test that would set too_close
- I deleted the lines that set h0 (which is never used)
- For real-valued integration paths (and only those) with b<a, the two ends are exchanged, waypoints are sorted, and the result is multiplied by -1 at the end
- And the two documentation inaccuracies are fixed.

What is not done is the documentation of the additional features this bug report is about: the capability to do array-valued integration.

(file #53231)

Michael Leitner <mleitner>
Thu 19 May 2022 04:22:02 PM UTC, comment #10: 

thanks for the complex analysis reminder. it's been a decade or so since i had any reason to use it and wasn't thinking about what I was actually doing with the waypoints. :)

confirming your last waypoint examples below:

>> quadgk(@(x) 1./x, -1+i,1+i, 'waypoints', [-1-.5i,1-.5i])
ans =
  -0.0000 + 4.7124i

>> quadgk(@(x) 1./x, -1+i,1+i, 'waypoints', [1-.5i,-1-.5i])
ans =
  -0.0000 - 7.8540i

i personally can't think of any reason without look at the code why we would need a<b always. maybe it simplifies the underlying algorithm? even if so, why not just swap a and b and sort the waypoints accordingly rather than recursively call quadgk and run through all the overhead twice.

Nicholas Jankowski <nrjank>
Project Member
Thu 19 May 2022 04:12:56 PM UTC, comment #9: 

No, as far as I can tell, Shampine does not use h0. It is always just (b-a), as there the test is written in untransformed coordinates. In transformed coordinates it simply becomes 1 or 2, which is here called h. One would have to do code archaeology to see when this definition of h0 was added to octave's code, and whether it was used for something at that time.

Michael Leitner <mleitner>
Thu 19 May 2022 03:51:04 PM UTC, comment #8: 

For the purely real case, waypoints are misleadingly named, they are just used as dividing points for the first division of the total interval into subintervals. The purpose is that if the integrand has integrable singularities in the interior, it will work much better if these are at the boundaries of the subintervals, so this is an opportunity for the user to give their positions and have better efficiency. Thus, it is sensible to use them in sorted order. It also would in principle work if the algorithm used them in any random order (also if some of them are outside [a,b]), as the corresponding contributions would just cancel (in exact arithmetics). In finite-precision arithmetics, it is of course sensible to have them always sorted and between a and b.

For the complex case, it is different. Here their purpose is to specify the integration path through the complex plane. Due to the residue theorem, the specific path you go is not relevant, all that matters is whether the poles are on the one or the other side of the path. So if you want to do some integral over the complex plane, you have to set the waypoints so that the connecting line from one to the next passes on the correct side (the one you want) of each pole. Specifically, if you do a closed contour integral around a pole, if you flip the sense of rotation, you get a minus in the result. So no, matlab's results are perfectly understandable, you would have seen a difference in your non-closed contour if you had had a choice of waypoints where you went under the pole at zero, but you didn't. Therefore you can contract your integration path to the direction connection between a and b in all cases by the residue theorem, and you get the same result.

For instance these two should give different results -- once you go effectively around the pole anticlockwise but miss the way from -1+i to 1+i, and once you effectively go once around the pole clockwise plus from -1+i to 1+i

quadgk(@(x) 1./x, -1+i,1+i, 'waypoints', [-1-.5i,1-.5i])
quadgk(@(x) 1./x, -1+i,1+i, 'waypoints', [1-.5i,-1-.5i])

Rik: waypoints are completely orthogonal to arrayvalued or not. We should add the sorting of waypoints in the purely real case, this should make the integration much more efficient (if they are not already sorted). And with respect to h0: I am still reading.

Michael Leitner <mleitner>
Thu 19 May 2022 03:25:08 PM UTC, comment #7: 

re waypoints:

according to the matlab 2022a help: 

When A, B,and the waypoints are all real, only the waypoints between A and B are used, and they are used in sorted order.
If A, B, or any entry of the waypoints vector is complex, the integration is performed over a sequence of straight line paths in the complex plane, from A to the first waypoint, from the first waypoint to the second, and so forth, and finally from the last waypoint to B.


>> quadgk(@(x) 1./x, 1,2)
ans =
>> quadgk(@(x) 1./x, 2,1)
ans =
>> quadgk(@(x) 1./x, 1,2,'waypoints',[1.25,1.5,1.75])
ans =
>> quadgk(@(x) 1./x, 2,1,'waypoints',[1.25,1.5,1.75])
ans =
>> quadgk(@(x) 1./x, 2,1,'waypoints',[1.75,1.5,1.25])
ans =


>> quadgk(@(x) 1./x, 1+i,2+i)
ans =
   0.4581 - 0.3218i
>> quadgk(@(x) 1./x,2+i,1+i)
ans =
  -0.4581 + 0.3218i
>> quadgk(@(x) 1./x, 1+i,2+i, 'waypoints', [1.25+1,1.5+i, 1.75+i])
ans =
   0.4581 - 0.3218i
>> quadgk(@(x) 1./x, 2+i,1+i, 'waypoints', [1.25+1,1.5+i, 1.75+i])
ans =
  -0.4581 + 0.3218i
>> quadgk(@(x) 1./x, 1+i,2+i, 'waypoints', [1.75+1,1.5+i, 1.25+i])
ans =
   0.4581 - 0.3218i
>> quadgk(@(x) 1./x, 2+i,1+i, 'waypoints', [1.75+1,1.5+i, 1.25+i])
ans =
  -0.4581 + 0.3218i

closed contour, using matlab help example:

>> f = @(z) 1./(2.*z-1);
contour_segments = [1+1i 0+1i 0-1i 1-1i];
q = quadgk(f,1,1,'Waypoints',contour_segments)
q =
  -0.0000 + 3.1416i
>> q = quadgk(f,1,1,'Waypoints',fliplr(contour_segments))
q =
  -0.0000 - 3.1416i

it seems waypoints ordering only matters for closed intervals. I would have thought the complex piecewise integration would have mattered, but it must still sort first, which seems odd, as it should cause a counting of doubly traversed paths. (but maybe those go negative and cancel? )

Nicholas Jankowski <nrjank>
Project Member
Thu 19 May 2022 03:09:50 PM UTC, comment #6: 

I tested the second version of quadgkv.m and it successfully handles the case I found in comment #1.

Is the "h0" variable used in the original Shampine paper?  I agree that it is unused in the code now, but maybe there was some reason for it that was not correctly translated from the academic paper when the original quadgk.m was written.

Rik <rik5>
Project Administrator
Thu 19 May 2022 03:05:37 PM UTC, comment #5: 

Does Matlab support reversed integration with Waypoints in integral()?  And does it also support Waypoints for ArrayValued functions?

It seems to me that maybe some of these combinations are esoteric and could simply be detected in the input validation and an error message delivered to the user to re-frame their problem.

Rik <rik5>
Project Administrator
Thu 19 May 2022 07:13:41 AM UTC, comment #4: 

Nicholas, that seems correct, however this made me see some other problems: first, if waypoints are given, they would have to be reversed also. This would probably make quite ugly code at this point, but I ask myself: is it necessary at all to have always a<b?  Apart from line 193, it seems to me that all the rest can equally well handle reversed subintervals (perhaps the isinf(a) (analogously isinf(b)) transformation implicitly assumes that if isinf(a) it is minus infinity, this would have to be looked into). A further problem: line 192 should test also whether a and b are real, not only waypoints.

And h0 is always only set and never used. This can definitely go.

And two documentation bugs: only when no waypoints and the limits are real are given, the interval is divided into 10 subintervals, and if the limits or any waypoint is complex, it is a piecewise linear (or straight) path, not only continuous.

Michael Leitner <mleitner>
Thu 19 May 2022 04:02:55 AM UTC, comment #3: 

I think since we made some warnings dependent on nargout, the b<a call to quadgkv around line 150 should probably become:

  if (b < a)
    ## Reverse integration
    if (nargout < 2)
      q = quadgkv (f, b, a, varargin{:});
      [q, err] = quadgkv (f, b, a, varargin{:});
    q = -q;

otherwise the error tol warning doesn't trigger:

>> q= quadgkv(@(t) -1 ./ t.^1.1,1,9999999999999)
warning: quadgkv: maximum interval count (650) exceeded
warning: called from
    quadgkv at line 407 column 9

warning: quadgkv: Error tolerance not met.  Estimated error 7.23189e-05
warning: called from
    quadgkv at line 418 column 7

q = -9.5319
>> q= quadgkv(@(t) -1 ./ t.^1.1,9999999999999,1)
warning: quadgkv: maximum interval count (650) exceeded
warning: called from
    quadgkv at line 407 column 9
    quadgkv at line 153 column 16

q = 9.5319

Nicholas Jankowski <nrjank>
Project Member
Tue 17 May 2022 09:27:54 AM UTC, comment #2: 

The problem should now be fixed. While doing that, I also simplified the body of the function by converting the output of the integrand first to a one-dimensional array, always work with this one-dimensional representation, and only at the end reshaping back. Also, I fixed the error message if the tolerance is not met.

(file #53218)

Michael Leitner <mleitner>
Tue 17 May 2022 03:16:00 AM UTC, comment #1: 

I think a discussion about how to handle the "too_close" condition belongs either on Discourse or in a new bug report.

I made a few changes for Octave coding conventions such as wrapping conditional in parentheses.  I also added the two vector tests from quadv.m.

For debug purposes, I have temporarily renamed the function quadgkv.m.  There seems to be a problem for some input functions.  I tried

octave:1> f = @(x) sin ((1:100) * x)
f =

@(x) sin ((1:100) * x)

octave:2> quadgkv (f, 0, 1)
error: product: nonconformant arguments (op1 is 2x100, op2 is 28x1)
error: called from
    quadgkv at line 446 column 11

Both quadv.m in Octave and integral in Matlab can integrate this function.

(file #53216)

Rik <rik5>
Project Administrator
Sun 15 May 2022 05:03:41 PM UTC, original submission:  

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.

Michael Leitner <mleitner>


(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

Attached Files
file #53257:  quadgk.m added by rik5 (25KiB - text/x-matlab)
file #53233:  quadgkv.m added by mleitner (23KiB - text/x-matlab)
file #53232:  quadgkv.m added by mleitner (23KiB - text/x-matlab)
file #53231:  quadgkv.m added by mleitner (23KiB - text/x-matlab)
file #53218:  quadgkv.m added by mleitner (24KiB - text/x-matlab)
file #53216:  quadgkv.m added by rik5 (24KiB - text/x-matlab)
file #53205:  quadgk.m added by mleitner (24KiB - text/x-matlab)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by rik5 (Updated the item)
  • -email is unavailable- added by mleitner (Submitted the item)

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only project members can vote.





    Follow 10 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-06-02 rik5 StatusPatch Reviewed Fixed
        Open/ClosedOpen Closed
    2022-05-26 rik5 Attached File- Added quadgk.m, #53257
    2022-05-19 mleitner Attached File- Added quadgkv.m, #53233
    2022-05-19 mleitner Attached File- Added quadgkv.m, #53232
    2022-05-19 mleitner Attached File- Added quadgkv.m, #53231
    2022-05-17 mleitner Attached File- Added quadgkv.m, #53218
    2022-05-17 rik5 StatusNone Patch Reviewed
    2022-05-17 rik5 Attached File- Added quadgkv.m, #53216
    2022-05-15 mleitner Attached File- Added quadgk.m, #53205

    Back to the top

    Powered by Savane 3.9