improve-stream-division

[Ralf Hemmecke]

I have some big computation that fails after about 5 hours claiming that
there is not enough stack space. That FriCAS was compiled with
"--dynamic-space-size 23906".

I tried to figure out where all that space was used and hit recip in
sttaylor.spad. In particular, the function iDiv looked suspicious to me.

iDiv : (ST A, ST A, A) -> ST A
iDiv(x, rst_y, ry0) == delay
empty? x => empty()
c0 := frst x * ry0
concat(c0, iDiv(rst x - c0*rst_y, rst_y, ry0))

In my understanding that means that for every new computed coefficient
of the result of recip(x), there are two new streams created, namely
"c0*rst_y" and "rst x - c0*rst_y".

I might be wrong, but that looked to me like wasting memory, because the
initial segments of these streams are not needed anymore, once the
respective coefficient is computed, but I suspect(ed) that the memory is
not freed.

Anyway, I suggest another implementation, that uses the inversion
formula directly, and only keeps a list of the computed coefficients in
reverse order for further computation.

With that implementation, my lengthy computation finished successfully.

The implementation can be found in my branch
"wip/improve-stream-division".

https://github.com/hemmecke/fricas/tree/wip/improve-stream-division
(patch attached)

I took care that the timings are basically the same.
See attached testfile.

Is there a chance that this could make it into FriCAS?

Suggestions for further improvements are welcome.

Ralf

[Waldek Hebisch]

On Sat, Jul 27, 2024 at 12:38:30AM +0200, Ralf Hemmecke wrote:
> I have some big computation that fails after about 5 hours claiming that
> there is not enough stack space. That FriCAS was compiled with
> "--dynamic-space-size 23906".
>
> I tried to figure out where all that space was used and hit recip in
> sttaylor.spad. In particular, the function iDiv looked suspicious to me.
>
> iDiv : (ST A, ST A, A) -> ST A
> iDiv(x, rst_y, ry0) == delay
> empty? x => empty()
> c0 := frst x * ry0
> concat(c0, iDiv(rst x - c0*rst_y, rst_y, ry0))
>
> In my understanding that means that for every new computed coefficient of
> the result of recip(x), there are two new streams created, namely
> "c0*rst_y" and "rst x - c0*rst_y".
>
> I might be wrong, but that looked to me like wasting memory, because the
> initial segments of these streams are not needed anymore, once the
> respective coefficient is computed, but I suspect(ed) that the memory is not
> freed.

AFAICS this is essentially the same problem that we had with
multiplication of power series (see doc/sttaylor.txt). Actual
reasons, look more subtle. Unreferenced things are supposed to be
garbage collected. Once 'iDiv' is done computing element number n
in the answer it should not look at earlier elements. Basically
after that there should be something like 2*n streams, but they
should take O(n) space because nobody is supposed to look at earlier
elements and streams should just reduce to functions which deliver
rest of the stream. Computing element n + 1 should adance those
2*n streams by one step, effectively discarding intermediate results
and only delivering the resulting coefficient. And of course,
2 more streams are created. At least when using sbcl empirical
observation indicates that memory use increases quadratically.

> Anyway, I suggest another implementation, that uses the inversion formula
> directly, and only keeps a list of the computed coefficients in reverse
> order for further computation.
>
> With that implementation, my lengthy computation finished successfully.
>
> The implementation can be found in my branch
> "wip/improve-stream-division".
>
> https://github.com/hemmecke/fricas/tree/wip/improve-stream-division
> (patch attached)
>
> I took care that the timings are basically the same.
> See attached testfile.
>
> Is there a chance that this could make it into FriCAS?

Yes, this looks reasonable. ATM main trouble I see are longish lines
and differen code style. For example instread of

restDiv(a: ST A, u: A, n: I, k: I, ra: ST A, b: ST A, revc: List A): ST A ==
delay
...

I would write something like

restDiv(a : ST A, u : A, n : I, k : I, ra : ST A, b : ST A,
revc : List A) : ST A == delay
...

It makes sense to run a bit more tests, I will do this and come back to
this in a day ar two.

--
Waldek Hebisch

[Ralf Hemmecke]

> AFAICS this is essentially the same problem that we had with
> multiplication of power series (see doc/sttaylor.txt).

Oh, cool. That is a good piece of explanation. Unfortunately, I haven't
seen it before.

And yes, the already used initial part of these streams might get
garbage collected, but I guess, it is better, to simply not create them
in the first place. Then also the garbage collector has nothing to do.

>> Is there a chance that this could make it into FriCAS?
>
> Yes, this looks reasonable.

Oh, that makes me happy.

> ATM main trouble I see are longish lines and differen code style.

Well, I usually go < 80, but if you prefer <72 that is also fine for me.
And yes, I like "a: T" much more than "a : T", but if you insist on the
latter, then it's also fine for me.

> It makes sense to run a bit more tests, I will do this and come back to
> this in a day ar two.

OK, I wait with the style modifications until then.

Ralf

[Martin R]

I'd be very interested in the actual computation.  I suspect that the sage implementation is much slower, but I'd be curious to know!

Best wishes,

Martin

[Ralf Hemmecke]

On 7/27/24 19:06, 'Martin R' via FriCAS - computer algebra system wrote:
> I'd be very interested in the actual computation. I suspect that the sage
> implementation is much slower, but I'd be curious to know!

Hi Martin,

I would bet that Fredrik Johansson has something faster, but that would
be for truncated power series.

Why do you claim that the sage implementation is slower? Did you mean
the pure-python implementation of infinite power series in python that
you (or your student) did?

Ralf

[Martin R]

It is only a suspicion, because FriCAS is generally very fast because it essentially compiles code, and python is generally rather slow.  However, I (and Travis Scrimshaw, who wrote that code with me), took some care to make it reasonable.

In any case, I'd like to play with it.

Best wishes,

Martin

[Waldek Hebisch]

On Sat, Jul 27, 2024 at 12:38:30AM +0200, Ralf Hemmecke wrote:
>

> Anyway, I suggest another implementation, that uses the inversion formula
> directly, and only keeps a list of the computed coefficients in reverse
> order for further computation.
>
> With that implementation, my lengthy computation finished successfully.
>
> The implementation can be found in my branch
> "wip/improve-stream-division".
>
> https://github.com/hemmecke/fricas/tree/wip/improve-stream-division
> (patch attached)
>
> I took care that the timings are basically the same.
> See attached testfile.
>
> Is there a chance that this could make it into FriCAS?
>
> Suggestions for further improvements are welcome.

With the patch applied the following gives wrong result (the
last thing should be stream of zero matrices, I get 5 zero
terms and then nonzero terms):

rI := SquareMatrix(2, FRAC(INT))
sT := Stream(rI)
sTO := StreamTaylorSeriesOperations(rI)
m1 := matrix([[0, 1],[0, 0]])$rI
m2 := matrix([[0, 0],[1, 0]])$rI
s1 := construct([3::rI, 0, m1, m2, 1])$sT
rs1 := recip(s1)$sTO;

s2 := construct([1::rI, 1::rI + m1, m2, m2])$sT
rs2 := recip(s2)$sTO;rI := SquareMatrix(2, FRAC(INT))

s3 := (*$sTO)(rs2, s1)
((/$sTO)(s3, s1) - rs2)$sTO


--
Waldek Hebisch

[Ralf Hemmecke]

Before I am going to fix this, I guess I need a discussion of what
x/y or 1/x is actually supposed to means in a non-commutative ring.

You probably realized that

((*$sTO)(rs1, s3) - rs2)$sTO

gives the same "wrong" result.

I have nothing against introducing right-/left- division, but / should
be reserved for the case when x * y^(-1) = y^(-1) * x otherwise I would
leave it undefined.

In other words we should have

if A has CommutativeRing then
"exquo" : (ST A,ST A) -> Union(ST A,"failed")
"/" : (ST A,ST A) -> ST A
recip : ST A -> UN

Maybe we can a bit weaker for recip, but I would still require that
x * recip(x) = recip(x) * x. Otherwise recip should fail.

Ralf

[Ralf Hemmecke]

On 7/31/24 16:05, Ralf Hemmecke wrote:
> You probably realized that
>
> ((*$sTO)(rs1, s3) - rs2)$sTO
>
> gives the same "wrong" result.

I meant "without applying the patch.

Ralf

[Waldek Hebisch]

On Wed, Jul 31, 2024 at 04:05:07PM +0200, Ralf Hemmecke wrote:
> Before I am going to fix this, I guess I need a discussion of what
> x/y or 1/x is actually supposed to means in a non-commutative ring.
>
> You probably realized that
>
> ((*$sTO)(rs1, s3) - rs2)$sTO
>
> gives the same "wrong" result.

It should be

((*$sTO)(s3, rs1) - rs2)$sTO

and this works without the patch.

> I have nothing against introducing right-/left- division, but / should be
> reserved for the case when x * y^(-1) = y^(-1) * x otherwise I would leave
> it undefined.

Well, division solves equation x = d*y, so we have d = x * y^(-1).
This definition is necessary to have sane interaction with multiplication.
To do it on the other side is a different operation. If you look
at definitions in the algebra, FreeDivisionAlgebra, Group and
XPolynomialRing are quite explicit and define x/y as x * y^(-1).
In most other cases order does not matter as corresponding
multiplication is commutative. Possibly the only unclear definition
is the one in StreamTaylorSeriesOperations.

> In other words we should have
>
> if A has CommutativeRing then
> "exquo" : (ST A,ST A) -> Union(ST A,"failed")
> "/" : (ST A,ST A) -> ST A
> recip : ST A -> UN
>
> Maybe we can a bit weaker for recip, but I would still require that
> x * recip(x) = recip(x) * x. Otherwise recip should fail.

We have:

recip : % -> Union(%,"failed")
++ recip(a) returns an element, which is both a left and a right
++ inverse of \spad{a},
++ or \spad{"failed"} if such an element doesn't exist or cannot
++ be determined (see unitsKnown).

So yes, x * recip(x) = recip(x) * x = e where e is the unit element
for multiplication.

And once you have recip "/" follows logically.

And we also have 'leftRecip' and 'rightRecip' for cases when inverse
works only on one side.

--
Waldek Hebisch

[Ralf Hemmecke]

On 7/31/24 16:49, Waldek Hebisch wrote:>> ((*$sTO)(rs1, s3) - rs2)$sTO

>>
>> gives the same "wrong" result.
>
> It should be
>
> ((*$sTO)(s3, rs1) - rs2)$sTO
>
> and this works without the patch.

I know, but why do you claim that it must be the way you say in a
non-commutative ring? I do not think that it is only my opinion that
this is ambiguous.

>> I have nothing against introducing right-/left- division, but /
should be
>> reserved for the case when x * y^(-1) = y^(-1) * x otherwise I
would leave
>> it undefined.

> Well, division solves equation x = d*y, so we have d = x * y^(-1).

How can you say so in a non-commutative ring? I guess some people claim
that it solves x=y*d.

> This definition is necessary to have sane interaction with
multiplication.
> To do it on the other side is a different operation. If you look
> at definitions in the algebra, FreeDivisionAlgebra, Group and
> XPolynomialRing are quite explicit and define x/y as x * y^(-1).
> In most other cases order does not matter as corresponding
> multiplication is commutative. Possibly the only unclear definition
> is the one in StreamTaylorSeriesOperations.

OK, you convinced me with our definitions in FriCAS. I still see danger
that our formatters transform x/y to \frac{x}{y} also for a
non-commutative domain (I haven't actually checked) and that would be
certainly wrong.

And yes, the docstring of / in StreamTaylorSeriesOperations must
definitely reflect compatibility with FreeDivisionAlgebra and all other
(non-commutative) places.

>> In other words we should have
>>
>> if A has CommutativeRing then
>> "exquo" : (ST A,ST A) -> Union(ST A,"failed")
>> "/" : (ST A,ST A) -> ST A
>> recip : ST A -> UN
>>
>> Maybe we can a bit weaker for recip, but I would still require that
>> x * recip(x) = recip(x) * x. Otherwise recip should fail.
>
> We have:
>
> recip : % -> Union(%,"failed")
> ++ recip(a) returns an element, which is both a left and a right
> ++ inverse of \spad{a},
> ++ or \spad{"failed"} if such an element doesn't exist or cannot
> ++ be determined (see unitsKnown).
>
> So yes, x * recip(x) = recip(x) * x = e where e is the unit element
> for multiplication.

Yes, but then we should also update the docstring in
StreamTaylorSeriesOperations

recip : ST A -> UN

++ recip(a) returns the power series reciprocal of \spad{a}, or
++ "failed" if not possible.

and make this explicit. I do not like that users (even though pretty
intuitive) have to guess that the word "reciprocal" actually refers to
the docstring of recip: % -> Union(%,"failed") (somewhere else) as you
cited above, i.e. that is a left+right inverse.

Should we introduce also a function "\" into StreamTaylorSeriesOperations?

Then I will come up with then necessary modifications in a new patch.

Ralf

[Waldek Hebisch]

On Wed, Jul 31, 2024 at 05:55:30PM +0200, Ralf Hemmecke wrote:
> On 7/31/24 16:49, Waldek Hebisch wrote:>> ((*$sTO)(rs1, s3) - rs2)$sTO
> >>
> >> gives the same "wrong" result.
> >
> > It should be
> >
> > ((*$sTO)(s3, rs1) - rs2)$sTO
> >
> > and this works without the patch.
> I know, but why do you claim that it must be the way you say in a
> non-commutative ring? I do not think that it is only my opinion that this is
> ambiguous.

Well, people use various weird conventions. But mainstream practice
to is chose conventions in a way that make computations simpler.
Normal rule is s*(x/y) = (s*x)/y. Another rule is (x/y)*y = x.
Both rules would be broken if you define division differently.
Or to put is differenly, if you want simple rules do not move things
to the other side on the whim.

> >> I have nothing against introducing right-/left- division, but / should be
> >> reserved for the case when x * y^(-1) = y^(-1) * x otherwise I would
> leave
> >> it undefined.
>
> > Well, division solves equation x = d*y, so we have d = x * y^(-1).
> How can you say so in a non-commutative ring? I guess some people claim that
> it solves x=y*d.
>
> > This definition is necessary to have sane interaction with multiplication.
> > To do it on the other side is a different operation. If you look
> > at definitions in the algebra, FreeDivisionAlgebra, Group and
> > XPolynomialRing are quite explicit and define x/y as x * y^(-1).
> > In most other cases order does not matter as corresponding
> > multiplication is commutative. Possibly the only unclear definition
> > is the one in StreamTaylorSeriesOperations.
> OK, you convinced me with our definitions in FriCAS. I still see danger that
> our formatters transform x/y to \frac{x}{y} also for a non-commutative
> domain (I haven't actually checked) and that would be certainly wrong.

If they do, we would need to fix this. But it is likely that formatters
do not produce 'x/y'.

> And yes, the docstring of / in StreamTaylorSeriesOperations must definitely
> reflect compatibility with FreeDivisionAlgebra and all other
> (non-commutative) places.

Yes.

OK. However, operations normally should be defined in categories
and users should look for definitions in categories. In some cases
like StreamTaylorSeriesOperations we want to keep functionality
as a package. Technically we can not reuse definition from categories
in such case, but users should understand that the operation are
the same or at least as close as possible do definitions in
categories. The definition above is from MagmaWithUnit, and is
inherited by Ring. StreamTaylorSeriesOperations in general
description says that it 'implements Taylor series arithmetic'.
I am not sure how explict we should be in the description, but
Taylor series _domains_ are declared to be a Ring, and "arithmetic"
operations defined in StreamTaylorSeriesOperations should satisfy
Ring properties. So, any deviations from definition in Ring either
are bugs or rewording wiht the same meaning or maybe a specialization
(ring of power series have some special properties compared to
deneral rings). From my point of view clearest docsting would say
"see corresponding operation in Ring".

To put this differently, at programming language level programmers
can implement quite a different property for operations with the same
name. But in many cases including StreamTaylorSeriesOperations
users should regard operations as "the same" as operations defined
in relevant category, that is Ring in this case.

> Should we introduce also a function "\" into StreamTaylorSeriesOperations?

If you wish. For me it is important to avoid breaking things
that used to work, so '/' is important. To say the truth, 'exquo'
is more important as is has uses in other places and is not the
the same as multiplication by reciprocal (which in case of 'exquo'
may be non-existent). Also, if you want tha add new operation
than something like 'leftExquo' probably would fit.

--
Waldek Hebisch

[Ralf Hemmecke]

Sorry that it took so long.

I hope that I did the changes to your liking.

The last 3 commits of

https://github.com/fricas/fricas/compare/master...hemmecke:fricas:wip/improve-stream-division

are the changes that I added on top of what I had in July.

All 4 commits in one patch is attached.

Ralf

[Waldek Hebisch]

On Sun, Oct 20, 2024 at 11:58:18PM +0200, 'Ralf Hemmecke' via FriCAS - computer algebra system wrote:
> Sorry that it took so long.
>
> I hope that I did the changes to your liking.
>
> The last 3 commits of
>
> https://github.com/fricas/fricas/compare/master...hemmecke:fricas:wip/improve-stream-division
>
> are the changes that I added on top of what I had in July.
>
> All 4 commits in one patch is attached.

Thanks, please commit.

--
Waldek Hebisch