Opened 21 months ago
Closed 14 months ago
#21131 closed enhancement (fixed)
Interpolation, Minimal Vanishing Polynomial and Multi Point Evaluation of Skew Polynomials
Reported by:  arpitdm  Owned by:  

Priority:  major  Milestone:  sage7.6 
Component:  coding theory  Keywords:  sd75 
Cc:  tscrim, caruso, jsrn, dlucas  Merged in:  
Authors:  Arpit Merchant  Reviewers:  Johan Rosenkilde 
Report Upstream:  N/A  Work issues:  
Branch:  4a9e4b0 (Commits)  Commit:  4a9e4b097681ffa0534918c01b4d31b41b38cfbf 
Dependencies:  #13215  Stopgaps: 
Description
Support for computing:
 Interpolation Polynomial
 Minimal Vanishing Polynomial
 Multi Point Evaluation
in Skew Polynomial Rings.
Change History (45)
comment:1 Changed 21 months ago by
 Branch set to u/arpitdm/mvp_mpe
comment:2 Changed 21 months ago by
 Commit set to c568c420198a1013d05eddd134f2fabb3bd8cdac
comment:3 Changed 21 months ago by
 Commit changed from c568c420198a1013d05eddd134f2fabb3bd8cdac to 282cfb1bf5a92744918b56d0190f6789295ca532
Branch pushed to git repo; I updated commit sha1. New commits:
282cfb1  added check for linear independence of evalpts in MVP

comment:4 followup: ↓ 6 Changed 21 months ago by
I've added methods for minimal_vanishing_polynomial
and multi_point_evaluation
. These work on skew polynomials over finite fields, both with sigma being the Frobenius or a power of Frobenius. MVP evaluates to 0 on all its eval_pts
and for MPE, evaluations match the individual calls. However, there is no improvement in running time of MPE as compared to the individual calls.
I've added interpolation_polynomial
. But its not correct yet. For instance, for any finite field, if I give 3 or more eval_pts and values, it fails.
Also, I am not sure that these same algorithms are valid for skew polynomials over general rings. I ran some examples and found that the MVP does not evaluate to zero, MPE does not match individual calls.
comment:5 followup: ↓ 16 Changed 21 months ago by
You could be more helpful by explaining *which* examples you ran, and tried to make them small and easy to debug. As it stands, I'm confused about what examples you did run, because MVP just threw an exception when I tried the "usual" nonfinite field example ZZ[t][x; t > t+1]
(see below).
minimal_vanishing_polynomial
generally outputs something in the skew polynomial over the fraction field of the base ring. I.e. ifS = R[x; sigma]
for some ringR
, then the minimal vanishing polynomial is defined overQ[x; sigma]
whereQ
is the fraction field ofR
. This is clear from the linex  (sigma(eval_pts[0]) / eval_pts[0])
, where the constant coefficient is clearly in the fraction field of the base ring.
I hadn't thought of that before. Add a check at the beginning of the call. If
R
is not a field, then construct the skew polynomial ring over the fraction field and return the result of the call from there. Do the same with the two other methods.
 In
multi_point_evaluation
you were doing something really strange in the case of one element. Just do the one evaluation (I already fixed it).
 In
interpolation_polynomial
for the 1point case you were doing quorem where you should just do the actual fraction. If the division doesn't go well, there is no way just taking the floored quotient will work! The result simply lives in the skew polynomial ring over the fraction field.
 Add in the description of the three methods that the evaluation points must be linearly independent over the fixed field of the twist map. Your current description of
eval_pts
is incorrect (the evaluation points are all in the base ring of self, and are therefore (almost) *never* linearly independent over that same ring).
 Add a doctest for a nonfinite field skew polynomials. Because of the above restriction, use e.g. a skew polynomial over the fraction field of ZZ[t] instead of ZZ[t].
 Add a doctest for
minimal_vanishing_polynomial
that throws the linearindependence error. Do the same formulti_point_evaluation
andinterpolation
.
 I removed a check from
interpolation_polynomial
. That check did not make any sense in most cases (the characteristic ofZZ[t]
is 0 for instance).
 In
interpolation_polynomial
you've now added a check to see if the evaluations match. Two comments: 1) why are you now using multipoint evaluation for this? 2) the check is currently being done in every recursive call ofinterpolation_polynomial
which is completely wasteful. Make the check only at the top level by creating a new inner function which is the actual recursive function. This top level can also be in charge of all the checks on the input lengths etc.
That's it for now, I think.
Best, Johan
comment:6 in reply to: ↑ 4 Changed 21 months ago by
Replying to arpitdm: However, there is no improvement in running time of MPE as compared to the individual calls.
What examples did you try? Which sizes of input? Which number of evaluation points? Which base fields? Which skew rings? There's many parameters here, and the speed profiles might be completely different over different choices. Please post your test code so your tests can be replicated.
If we can't find a single set of parameters for which the new multipoint evaluation is fastest, then we should just replace the method with the oneliner return [ p(e) for e in eval_pts ]
. We could still keep the code for posterity, since e.g. after the Karatsuba implementation, a strategy such as this *should* give an improvement for large enough skew polynomials (over finite fields).
Best, Johan
comment:7 Changed 21 months ago by
 Branch changed from u/arpitdm/mvp_mpe to u/jsrn/mvp_mpe
comment:8 Changed 21 months ago by
 Commit changed from 282cfb1bf5a92744918b56d0190f6789295ca532 to 6618340ece629db87b92d77b0c7ae0e482480791
forgot to push my changes.
New commits:
3166c6d  Merge branch 'u/arpitdm/mvp_mpe' of git://trac.sagemath.org/sage into 21131_interpolation_skew_poly

7806d86  Fixed bug in interpolation

e1cce49  Improved doc test for multipoint evaluation

b683f55  dividebytwo potential error. Improved an error message

012c9d3  rm nongeneral sanity check

9064e0e  Fix division bug in interpolation

6618340  Fix dividebytwo potential problem in interpolation

comment:9 Changed 21 months ago by
Here are some mathematical comments:
 I would say (but I'm not sure) that the benefit of using a divideandconquer method for multipoint evaluation only appears if you are using fast multication algorithms
 the problem of "minimum subspace polynomial" reduces to the computation of some lcm: the minimum subspace polynomial of a family (a_1, ..., a_n) (with a_i nonzero for all i) is the left lcm of the degree 1 skew polynomials (X  sigma(a_i)/a_i). Again I think that computing it using a divideandconquer method is only interesting when fast multiplication and halfgcd are available (but these should be hidden at some point in WachterZeh algorithm)
 the interpolation problem is just a noncommutative CRT: finding P such that P(a_i)=b_i is equivalent to finding a skew polynomial P which is congruent to b_i modulo (X  sigma(a_i)/a_i). I think that it would be interesting in any case to implement a method CRT (and then possibly let interpolation call it, but of course it is not mandatory).
comment:10 followup: ↓ 11 Changed 21 months ago by
 Commit changed from 6618340ece629db87b92d77b0c7ae0e482480791 to 0635ec993a6434362e13c6d2b617670f11c6a6c8
Branch pushed to git repo; I updated commit sha1. New commits:
362fff7  Merge branch 'u/arpitdm/skew_polynomials' of git://trac.sagemath.org/sage into 13215_skew_polynomials

4a33c5f  mv skew_polynomial_element in module_list

a921040  Changed unfortunate parameter name

f43767a  Some minor doc and error message modifications

7c5c511  Modelled SkewPolynomialBaseringInjection more closely on PolynomialBaseringInjection

d6fd904  Fix Cython warning

7134b7e  Removed _call_with_args

291c28c  Fixed construct0 bug

ca15975  Fix doctest

0635ec9  Merge branch '13215_skew_polynomials' into 21131_interpolation_skew_poly

comment:11 in reply to: ↑ 10 Changed 21 months ago by
Merged the newest 13215 into this branch.
comment:12 Changed 21 months ago by
Hi Xavier,
Great to see that you're listening in :)
Replying to caruso:
Here are some mathematical comments:
 I would say (but I'm not sure) that the benefit of using a divideandconquer method for multipoint evaluation only appears if you are using fast multication algorithms
I believe you're right, but I only thought about it after Arpit had already implemented it. So now I think we might as well look at its running time rather than guessing.
 the problem of "minimum subspace polynomial" reduces to the computation of some lcm: the minimum subspace polynomial of a family (a_1, ..., a_n) (with a_i nonzero for all i) is the left lcm of the degree 1 skew polynomials (X  sigma(a_i)/a_i). Again I think that computing it using a divideandconquer method is only interesting when fast multiplication and halfgcd are available (but these should be hidden at some point in WachterZeh algorithm)
Hmm, yes, I hadn't thought about that. However, I don't see why doing successive left_lcm
wouldn't be cubic complexity? It is easy to do a divide & conquerversion of the successive lcm (see below) which should be quadratic however:
def mvp_lcm(S, pts): x = S.gen() p = S.one() for a in pts: q = x  sigma(a)/a p = p.left_lcm(q) return p def mvp_lcm_dc(S, pts): x = S.gen() def rec(pts): if len(pts) > 1: t = len(pts)//2 left = mvp_lcm_dc(S, pts[:t]) right = mvp_lcm_dc(S, pts[t:]) return left.left_lcm(right) else: return x  sigma(pts[0])/pts[0] return rec(pts)
I did some experiments. Over finite fields, the D&C lcm wins:
m = 400 k.<a> = GF(2^m,'a') sigma = k.frobenius_endomorphism() S.<x> = k['x', sigma] pts = [ k.random_element() for i in range(m//2) ] assert not(0 in set(pts)) , "zero not allowed in the set (try again)" assert len(set(pts)) == len(pts) , "the set contains duplicates (try again)" %timeit S.minimal_vanishing_polynomial(pts) 1 loop, best of 3: 756 ms per loop %timeit mvp_lcm(S, pts) 1 loop, best of 3: 773 ms per loop %timeit mvp_lcm_dc(S, pts) 1 loop, best of 3: 369 ms per loop
However, over a field with coefficientgrowth, the implemented version wins:
m = 10 R.<t> = QQ[] QR = R.fraction_field() t = QR(t) sigma = QR.hom([t+1]) S.<x> = QR['x', sigma] pts = [ QR.random_element(degree=3) for i in range(m//2) ] assert not(0 in set(pts)) , "zero not allowed in the set (try again)" assert len(set(pts)) == len(pts) , "the set contains duplicates (try again)" %timeit S.minimal_vanishing_polynomial(pts) 10 loops, best of 3: 108 ms per loop %timeit mvp_lcm(S, pts) 1 loop, best of 3: 386 ms per loop %timeit mvp_lcm_dc(S, pts) 1 loop, best of 3: 835 ms per loop
The above timings vary a lot with the size of the fractions that were randomly drawn, but the relative ranking between the methods seems to stay put.
 the interpolation problem is just a noncommutative CRT: finding P such that P(a_i)=b_i is equivalent to finding a skew polynomial P which is congruent to b_i modulo (X  sigma(a_i)/a_i). I think that it would be interesting in any case to implement a method CRT (and then possibly let interpolation call it, but of course it is not mandatory).
That's a good idea. Feel free to do that. We are on a pretty tight schedule with Arpit's Google Summer of Code, so that won't be a priority for us now.
Best, Johan
comment:13 followup: ↓ 14 Changed 21 months ago by
I just made some more experiments: changing multi_point_evaluation
into the trivial return [ p(e) for e in eval_pts ]
then minimal_vanishing_polynomial
is much faster. It beats the other two alternatives in both domains.
I think for now, multi_point_evaluation
should be changed to just that simple function. As an orthogonal concern, I think that method should be a method on skew polynomials, not on the ring.
comment:14 in reply to: ↑ 13 ; followup: ↓ 15 Changed 21 months ago by
Replying to jsrn:
I just made some more experiments: changing
multi_point_evaluation
into the trivialreturn [ p(e) for e in eval_pts ]
thenminimal_vanishing_polynomial
is much faster. It beats the other two alternatives in both domains.
But maybe the fast multiplication and related methods once available, will improve the speed. I am not sure it makes sense to have a method to simply call the evaluate method n
times. Users can do that themselves. I think we should keep the D&C implementation and use the trivial calling wherever we are using it right now, with a note that says call MPE once that is fast.
I think for now,
multi_point_evaluation
should be changed to just that simple function. As an orthogonal concern, I think that method should be a method on skew polynomials, not on the ring.
Agreed. It is a method every skew polynomial should have if it is fast.
comment:15 in reply to: ↑ 14 Changed 21 months ago by
Replying to arpitdm:
But maybe the fast multiplication and related methods once available, will improve the speed. I am not sure it makes sense to have a method to simply call the evaluate method
n
times. Users can do that themselves. I think we should keep the D&C implementation and use the trivial calling wherever we are using it right now, with a note that says call MPE once that is fast.
I disagree: if the method is there, it should in all cases be at least as fast
as what a user could simply do himself. It would be nice if we could find a skew polynomial ring for which the current implementation is fastest, but it probably doesn't exist. Therefore, we have, IMHO, two options: 1) not having the method at all; or 2) having the method with the trivial implementation. I vote for the latter because of three things: A) a fast implementation will come at one point (hopefully; B) minimal_vanishing_polynomial
wouldn't have to change once a fast implementation of multipoint evaluation is there; and C) it is good to advertise multipoint evaluation to users, so that their code uses that, and will therefore automatically be fast with fast implementations.
I think for now,
multi_point_evaluation
should be changed to just that simple function. As an orthogonal concern, I think that method should be a method on skew polynomials, not on the ring.Agreed. It is a method every skew polynomial should have if it is fast.
OK
Best, Johan
comment:16 in reply to: ↑ 5 ; followup: ↓ 17 Changed 21 months ago by
Replying to jsrn:
I hadn't thought of that before. Add a check at the beginning of the call. If
R
is not a field, then construct the skew polynomial ring over the fraction field and return the result of the call from there. Do the same with the two other methods.
if not isinstance(R, Field): Q = R.fraction_field() t = Q(R.gen()) S = Q[self.variable_name(), sigma]
For example if we have some twist map of the original skew polynomial ring based on R sigma
, what's the proper way of converting to the equivalent twist map over the fraction field of R?
comment:17 in reply to: ↑ 16 Changed 21 months ago by
Replying to arpitdm:
Replying to jsrn:
if not isinstance(R, Field): Q = R.fraction_field() t = Q(R.gen()) S = Q[self.variable_name(), sigma]For example if we have some twist map of the original skew polynomial ring based on R
sigma
, what's the proper way of converting to the equivalent twist map over the fraction field of R?
That's a good question. Off the top of my hat, I don't immediately see how to do that in a bulletproof fashion. However, the following should work for all finitely generated rings: construct a sigma_frac
from sigma
as the homomorphism whose image on the the generators is as sigma
. Something like (untested):
Q = R.fraction_field() gens = R.gens() sigma_frac = Q.hom([ Q(sigma(g)) for g in gens ])
You'd better wrap that in a tryexcept and throw a ValueError: Unable to lift the twist map to a twist map over <insert Q here>
or something.
Best, Johan
comment:18 Changed 21 months ago by
 Branch changed from u/jsrn/mvp_mpe to u/arpitdm/mvp_mpe
comment:19 Changed 21 months ago by
 Commit changed from 0635ec993a6434362e13c6d2b617670f11c6a6c8 to 8f42a338565e35ce2210216c9ec9e7a240f97fe9
I've made the following changes based on discussions above:
 created twist map over the fraction field. fixed methods to accommodate fraction fields.
 refactored mvp and interpolation to have inner functions that are called recursively so that the checks on the arguments are made only at the top level.
 changed method
MPE
to the trivial repeated calling of the evaluate method and moved this toclass Skew_Polynomial
.
I'm opening the ticket for review now.
Best, Arpit.
Last 10 new commits:
c034df8  removed unused private methods and imports

ac109c8  converted to standard copyright template

5137ef9  removed deprecated method getslice

b1e42c3  Merge commit 'aca3398a81b3688e1f2e69c5910b8214d13be925' of git://trac.sagemath.org/sage into skew_polynomials

f563aba  mod and floordiv modified to use coercion model. replaced type with isinstance.

26e4689  modified copyright banner

6c86537  merged updated code from 13215. fixed interpolation to accomodate fraction fields. updated documentation.

50f9bdd  refactored mvp and interpolation to have inner functions so that checks on the arguments are made only at the top level.

dd667ed  changed MPE to the trivial repeated calling of the evaluation function. added a todo block to the documentation.

8f42a33  made MPE a method on skew polynomials.

comment:20 Changed 21 months ago by
 Status changed from new to needs_review
comment:21 Changed 21 months ago by
 Status changed from needs_review to needs_work
The three methods are mutually recursive, so you're still for instance recreating the same polynomial ring over and over again in the recursive calls. You should make three toplevel private functions, without the checks, and all the mutually recursive calls will call those toplevel private functions, taking the correct target ring as an input.
The creation of the fractionfield skewpolyring should be a helper function since it's duplicated code. But make it a private function.
Best, Johan
comment:22 Changed 21 months ago by
 Commit changed from 8f42a338565e35ce2210216c9ec9e7a240f97fe9 to 34bebbc6eb951d7e459ee2ee9988461d89af27d3
Branch pushed to git repo; I updated commit sha1. New commits:
f25f886  created top level private function for MVP

88196bf  some more fixes to the documentation.

b9f8239  merging latest changes

34bebbc  created top level functions for interpolation and fraction field skew polynomial ring. added documentation and indirect doctests.

comment:23 followup: ↓ 24 Changed 21 months ago by
I created three new private functions namely _gen_one_sigma
, _create_mvp
and _interpolate
. Are these names alright?
Best, Arpit.
comment:24 in reply to: ↑ 23 Changed 21 months ago by
Replying to arpitdm:
I created three new private functions namely
_gen_one_sigma
,_create_mvp
and_interpolate
. Are these names alright?
_gen_one_sigma
is a pretty terrible name. I previously suggested _base_ring_to_fraction_field
, but perhaps that mail was accidentally only "Reply" to David instead of "Reply all".
_create_mvp
is suggest _minimal_vanishing_polynomial
, i.e. exactly the same name as the main function, but with an underscore.
Best, Johan
comment:25 Changed 21 months ago by
 Commit changed from 34bebbc6eb951d7e459ee2ee9988461d89af27d3 to b523806a461e61036de748c7d7598044acc85892
Branch pushed to git repo; I updated commit sha1. New commits:
b523806  changed names of private methods

comment:26 Changed 21 months ago by
 Status changed from needs_work to needs_review
comment:27 Changed 21 months ago by
Hello,
I just have a few comments:
interpolation_polynomial
:
 It does not respect the docstring format: there should be a line break after "Return the interpolation polynomial"
 It's not much, but I think it would be a bit better if the description of
eval_pts
was the same in the documentation of_interpolation
` andinterpolation_polynomial
minimal_vanishing_polynomial
:
 Same remark as above regarding the docstring format...
 And regarding "consistency" between description of input arguments.
Best,
David
comment:28 Changed 21 months ago by
 Commit changed from b523806a461e61036de748c7d7598044acc85892 to d403fae361316c8da23e388381e69148863ae48b
Branch pushed to git repo; I updated commit sha1. New commits:
eb78e69  small fixes to docstrings

2d67e0e  merging updates

a2c4f06  removed unused imports, signal statements. small fixes to documentation.

e435f56  Merge branch 'u/arpitdm/mvp_mpe' of git://trac.sagemath.org/sage into mvp_mpe

d403fae  fixes to documentation

comment:29 Changed 21 months ago by
Hi David,
I've made the changes you proposed. Opening it for review again.
Best, Arpit.
comment:30 Changed 20 months ago by
 Branch changed from u/arpitdm/mvp_mpe to u/jsrn/mvp_mpe
comment:31 followup: ↓ 33 Changed 20 months ago by
 Commit changed from d403fae361316c8da23e388381e69148863ae48b to 4bb82b8c5ffff8e442ef5aad6a224b590c3f9331
 Status changed from needs_review to needs_work
I made a number of modifications:
 Renamed
interpolation
tolagrange_polynomial
to matchPolynomialRing
. Also changed the calling convention for this.
 Simplified
_base_ring_to_fraction_field
: it is much simpler and nicer to just return the constructed ring.
 Moved the helper functions outside the class.
 Allow MVP to get input which is linearly dependent over the fixed field of the twist map. In this case, the degree will simply be lower. Removed the
check
argument.
 Removed the
check
argument fromlagrange_polynomial
and simply detect the linear dependence (an evaluation point will become 0 during a recursive call).
 General improvements to docstrings and tests
I tried for a while to make lagrange_polynomial
accept linearly dependent evaluation points, and then only fail if the values were not similarly dependent (recall that if e.g. the last evaluation point is linearly dependent on the rest, then the value that *any* skew polynomial takes on that last evaluation point is the corresponding linear combination of the value it takes on the previous evaluation points). But the algorithm doesn't seem to be able to detect *what* the linear dependence is in a convenient way, so I gave up and decided to simply throw an exception.
Please check the compiled doc. I have no more time today.
Best, Johan
New commits:
5041176  Merge branch 'u/arpitdm/mvp_mpe' of git://trac.sagemath.org/sage into 21131_interpolation_skew

ce8b6b4  mpe: modified TODO and improved doc test

256052a  minimum > minimal

2707768  clarify a doctest

3a1b992  Moved _functions outside class. Simplified call conventions

0ed0179  Allow minimal_vanishing_polynomial on input which is linearly dependent

db72757  interpolation > lagrange_polynomial and make call convention as PolynomialRing

115cf04  Improved doc and example of _base_ring_to_fraction_field

ff371ab  lagrange: improve doc. Rm check param and instead discover linear dependence

4bb82b8  Improve a bit the code and doc

comment:32 Changed 20 months ago by
 Keywords sd75 added
comment:33 in reply to: ↑ 31 Changed 19 months ago by
Based on discussions, we are putting the SkewPolynomial_finite_field class on hold. The multi_point_evaluation
method is in the skew_polynomial_element.pyx
file and the interpolation
and minimum_vanishing_polynomial
methods are in the skew_polynomial_ring.py
file. This ticket does not depend on the finite fields ticket and the last discussion (during SD75) was that the history would have to be rewritten. I just want to confirm that that's still the direction to take and then I will push the changes.
comment:34 Changed 19 months ago by
Agreed, let's go ahead with it!
comment:35 Changed 19 months ago by
 Branch changed from u/jsrn/mvp_mpe to u/arpitdm/mvp_mpe
comment:36 Changed 19 months ago by
 Commit changed from 4bb82b8c5ffff8e442ef5aad6a224b590c3f9331 to 0572e852c1be546810e3f53ef08b569cb592415a
comment:37 Changed 19 months ago by
 Dependencies changed from #21088, #13215 to #13215
comment:38 Changed 19 months ago by
 Status changed from needs_work to needs_review
comment:39 Changed 18 months ago by
Coming back to it now, I feel sort of silly finalising the review myself, since the latest version was written by me. I feel that the last step is that someone should review my modifications.
I've tried to understand what you did in commit 0572e85, but it's not clear to me. All the finite field stuff is still showing up in the log, and have just been removed again  did you just remove them in the merge commit 0572e85?
Best, Johan
comment:40 Changed 18 months ago by
 Branch changed from u/arpitdm/mvp_mpe to u/jsrn/mvp_mpe
comment:41 Changed 18 months ago by
 Commit changed from 0572e852c1be546810e3f53ef08b569cb592415a to 4a9e4b097681ffa0534918c01b4d31b41b38cfbf
I fixed some small doc markup issues. Otherwise the code looks good to me.
New commits:
4a9e4b0  Doc markup fixes

comment:42 Changed 15 months ago by
 Status changed from needs_review to positive_review
comment:43 Changed 15 months ago by
 Milestone changed from sage7.3 to sage7.6
 Reviewers set to Johan Rosenkilde
comment:44 Changed 15 months ago by
Thanks
comment:45 Changed 14 months ago by
 Branch changed from u/jsrn/mvp_mpe to 4a9e4b097681ffa0534918c01b4d31b41b38cfbf
 Resolution set to fixed
 Status changed from positive_review to closed
Branch pushed to git repo; I updated commit sha1. New commits:
added interpolation polynomial method