RIF/CIF: search and subintervals by a bisection-algorithm

[dkrenn]

This ticket implements a bisection algorithm used with interval field elements.

CC: @cheuberg

Component: numerical

Author: Daniel Krenn

Branch/Commit: u/dkrenn/rif/bisect @ 4f4b2c9

Issue created by migration from https://trac.sagemath.org/ticket/19033

[dkrenn]

Branch: u/dkrenn/rif/bisect

[dkrenn]

Commit: 788fb33

[dkrenn]

New commits:

457fe87	write bisect function
61b296b	write RIF/CIF.find_subintervals
788fb33	small improvement in docstring

[videlec]

comment:3

Hello Daniel,

This is cool!

I guess that the line

verbose('iteration %s with results in %s of %s cells' %
   (iteration, len(result), len(open)), level=2)

is quite bad for performances. The string formatting is done whatever the verbose function is doing. Did you do some profiling?

Just for curiosity, do you know how good is fast_callable with interval fields? It looks pretty bad on the following example

sage: f(x) = log(exp(x*sin(x)) + exp(x*cos(x)))
sage: F = fast_callable(f, domain=RIF)
sage: r = RIF(0.1,0.2)
sage: %timeit F(r)
10000 loops, best of 3: 54.6 µs per loop
sage: %timeit ((r*r.sin()).exp() + (r*r.cos()).exp()).log()
10000 loops, best of 3: 32.8 µs per loop

compared to reals

sage: F = fast_callable(f, domain=RR)
sage: r = 0.1
sage: %timeit F(r)
100000 loops, best of 3: 12.4 µs per loop
sage: %timeit ((r*r.sin()).exp() + (r*r.cos()).exp()).log()
10000 loops, best of 3: 16.5 µs per loop

Vincent

[dkrenn]

comment:4

Hi Vincent,

Replying to @videlec:

I guess that the line

verbose('iteration %s with results in %s of %s cells' %
   (iteration, len(result), len(open)), level=2)

is quite bad for performances. The string formatting is done whatever the verbose function is doing. > Did you do some profiling?

No, I didn't do any with this line included (I did on my old code, which used if and print für debugging).

Just for curiosity, do you know how good is fast_callable with interval fields?

Don't know. I used it with some larger symbolic expressions, and what I remember, I had a significant speedup there.

It looks pretty bad on the following example

sage: f(x) = log(exp(x*sin(x)) + exp(x*cos(x)))
sage: F = fast_callable(f, domain=RIF)
sage: r = RIF(0.1,0.2)
sage: %timeit F(r)
10000 loops, best of 3: 54.6 µs per loop
sage: %timeit ((r*r.sin()).exp() + (r*r.cos()).exp()).log()
10000 loops, best of 3: 32.8 µs per loop

compared to reals

sage: F = fast_callable(f, domain=RR)
sage: r = 0.1
sage: %timeit F(r)
100000 loops, best of 3: 12.4 µs per loop
sage: %timeit ((r*r.sin()).exp() + (r*r.cos()).exp()).log()
10000 loops, best of 3: 16.5 µs per loop

Not so good...

Any suggestions on the bisect-code using fast_callable?

Daniel

[videlec]

comment:5

Replying to @dkrenn:

Hi Vincent,

Replying to @videlec:

I guess that the line

verbose('iteration %s with results in %s of %s cells' %
   (iteration, len(result), len(open)), level=2)

is quite bad for performances. The string formatting is done whatever the verbose function is doing. > Did you do some profiling?

No, I didn't do any with this line included (I did on my old code, which used if and print für debugging).

If you really want these verbose you should use lazy formatting. I do not remember what is the state of the art, but you can have a look at sage.misc.lazy_string and sage.misc.lazy_format. But you are really creating a string within the critical loop. I do find it not very reasonable (the other is fine though).

Just for curiosity, do you know how good is fast_callable with interval fields?

Don't know. I used it with some larger symbolic expressions, and what I remember, I had a significant speedup there.

I am pretty sure that fast_callable is better than symbolic and you seem to confirm that. However, you can see in my small snippet above that fast_callable looks worse than using Python. I should try on a larger expression but then you have to translate into an object call oriented expression and I am lazy...

Any suggestions on the bisect-code using fast_callable?

Check what fast_callable is doing with domain=RIF and domain=CIF (quickly looking at ext/fast_callable.pyx I would say nothing). It can be done independently of this ticket.

Other comments:

• What is your variable result for?

• did you try to more typing. In other words open = [] -> cdef list open = [] and iteration = 0 -> cdef size_t iteration = 0, etc. Cython does good job with Python list when it knows that these are lists.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

a97ddce

reintroduced flag bisect_on_success

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 788fb33 to a97ddce

[dkrenn]

comment:7

Replying to @videlec:

• What is your variable result for?

It collects the already finished cells. I've rewritten the code and introduced a new option (this was already there in a pre-version of this code).

[fchapoton]

comment:8

Bad syntax for INPUT::, it should be INPUT:

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from a97ddce to ac3428c

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

dfb229a	doc-build: fix INPUT-block
de2633d	fix broken links in existing code
bd9cfa1	fix broken links in bisect
ac3428c	increase performance

[dkrenn]

comment:10

Replying to @videlec:

If you really want these verbose you should use lazy formatting. I do not remember what is the state of the art, but you can have a look at sage.misc.lazy_string and sage.misc.lazy_format. But you are really creating a string within the critical loop. I do find it not very reasonable (the other is fine though).

LazyFormat did not bring something on the performance. Now checking verbosity level directly; this is faster now.

[dkrenn]

comment:11

Replying to @videlec:

• did you try to more typing. In other words open = [] -> cdef list open = [] and iteration = 0 -> cdef size_t iteration = 0, etc. Cython does good job with Python list when it knows that these are lists.

Done. Is faster now.

[dkrenn]

comment:12

Replying to @fchapoton:

Bad syntax for INPUT::, it should be INPUT:

Corrected. Also fixed some broken links in the files (in existing code).

Again, needs_review :)

[jdemeyer]

comment:13

Unless I'm missing something, the case f is None is completely undocumented and untested.

[jdemeyer]

comment:14

Also, I guess the test is meant to satisfy the assumption that if an interval A passes, any interval containing A should also pass. This should be documented.

[jdemeyer]

comment:15

Another detail: you can replace

from sage.rings.real_mpfi import bisect
return bisect(f, self, test, **kwds)

by

return real_mpfi.bisect(f, self, test, **kwds)

[jdemeyer]

comment:16

Replying to @videlec:

• did you try to more typing. In other words open = [] -> cdef list open = [] and iteration = 0 -> cdef size_t iteration = 0, etc. Cython does good job with Python list when it knows that these are lists.

I agree, there should be a lot more typing. All variables which are boolean should be typed bint and all integers a suitable integer type (Py_ssize_t for lengths).

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

565c482	Merge tag '6.9' into t/19033/rif/bisect
dfbbccb	Trac #19033, comment 15: simplify import
1381c29	Trac #19033, comment 16: types specified
6ccfcf0	Trac #19033, comment 13: remove case f=None
dd7cb13	Trac #19033. comment 14: add note on monotonicity of parameter test

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from ac3428c to dd7cb13

[dkrenn]

comment:18

Replying to @jdemeyer:

Unless I'm missing something, the case f is None is completely undocumented and untested.

Is now removed (since not needed).

[dkrenn]

comment:19

Replying to @jdemeyer:

Also, I guess the test is meant to satisfy the assumption that if an interval A passes, any interval containing A should also pass. This should be documented.

I've added a note block explaining this.

[dkrenn]

comment:20

Replying to @jdemeyer:

Another detail: you can replace

from sage.rings.real_mpfi import bisect
return bisect(f, self, test, **kwds)

by

return real_mpfi.bisect(f, self, test, **kwds)

Done.

[dkrenn]

comment:21

Replying to @jdemeyer:

Replying to @videlec:

• did you try to more typing. In other words open = [] -> cdef list open = [] and iteration = 0 -> cdef size_t iteration = 0, etc. Cython does good job with Python list when it knows that these are lists.

I agree, there should be a lot more typing. All variables which are boolean should be typed bint and all integers a suitable integer type (Py_ssize_t for lengths).

I've added cdef bint success; this saves about 5 percent in time in my tested examples. I've tried to add Py_ssize_t and other bint...it does not seem that this brought something. Maybe I am using it wrong...

[dkrenn]

comment:22

I've also merged in 6.9 (I am sorry for this, since it was not needed, but I didn't have the previous version available and did not want to want until it was compiled).

Set it back to needs_review.

[fchapoton]

comment:23

one failing doctest

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from dd7cb13 to 99c3e9d

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

99c3e9d

Revert "Trac #19033, comment 15: simplify import" since imports changed on 6.10.beta0

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 99c3e9d to faacb5f

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

faacb5f

Merge tag '6.10.beta0' into t/19033/rif/bisect

[dkrenn]

comment:26

Replying to @fchapoton:

one failing doctest

Was not there on 6.9.

Merged 6.10.beta0 and reverted the changes from comment 15 of this ticket (change of imports). Back to needs_review.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

0dfca39

Merge tag '7.1.beta3' into t/19033/rif/bisect

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from faacb5f to 0dfca39

[dkrenn]

comment:28

Merged 7.1.beta3.

[fchapoton]

comment:29

does not apply

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

d378b1f

Merge tag '7.1.beta4' into t/19033/rif/bisect

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 0dfca39 to d378b1f

[dkrenn]

comment:31

Replying to @fchapoton:

does not apply

Merged in 7.1.beta4.

[fchapoton]

comment:32

does not apply

[videlec]

comment:33

Note that arb already offers a bissection/newton scheme for root finding of analytic functions, see arb_calc.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

4f4b2c9

Merge tag '8.6' into u/dkrenn/rif/bisect

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from d378b1f to 4f4b2c9

[dkrenn]

comment:35

Replying to @videlec:

Note that arb already offers a bissection/newton scheme for root finding of analytic functions, see arb_calc.

Good to know, thanks.

[dkrenn]

comment:36

Merged in 8.6, so back to needs review again.

When I see this correctly, all comments given on this ticket so far have been incorporated or answered, but noone gave a positive review three years ago.....I guess, however, we are close to a positive review after all. Can someone have a look again, please.

[videlec]

comment:38

[comment:15] is not addressed. real_mpfi is already globally cimported in complex_interval.pyx.

[videlec]

comment:39

I don't think that "bisection on success"/"bisection on failure" is enough. In many instances you have better than a boolean test. Let me consider the situation where you want to isolate all roots of the function f (e.g. a square-free polynomial). Interval arithmetic allows you to:

• prove that a cell contains a single zero (if f(left) and f(right) have different signs and f'(cell) does not vanish)
• prove that a cell does not contain a zero (if f(cell) does not contain zero) In this situation, you want to perform bisection if the two above certificate failed. It is a True/False/Unknown situation.

I think that the bisection should be general enough to handle this search because it does provide a proof of what you are looking for. The return value would be a pairs of lists: the one you are interested in and the pieces for which bisection failed to determinate their status.

Also, isolating roots in an example such as

sage: def contains_zero(fct, cell):
....:     return fct(cell).contains_zero()
sage: result = bisect(lambda x: x^3-4*x+2, RIF(-10, 10), contains_zero)

would be faster since you would not refine the already valid cells.

[videlec]

comment:40

(for all this, you should have a look at arb_calc_isolate_roots that was already mentioned)

[videlec]

comment:41

The following is a big waste.

if join_neighboring_cells:
     verbose('joining neighboring cells...', level=1)
     cells = result
     result = []

     while len(cells) > 0:
         joined = cells.pop(0)
         k = 0
         while k < len(cells):
             try:
                 joined.intersection(cells[k])
             except ValueError:
                 k += 1
             else:
                 joined = joined.union(cells.pop(k))
         result.append(joined)

The intervals should be stored in order (a binary tree seems the most appropriate).

[dkrenn]

comment:43

Replying to @videlec:

[comment:15] is not addressed. real_mpfi is already globally cimported in complex_interval.pyx.

Doesn't work:

        return real_mpfi.bisect(f, self, test, **kwds)
                       ^
------------------------------------------------------------

sage/rings/complex_interval.pyx:2201:24: cimported module has no attribute 'bisect'