Opened 6 years ago
Closed 5 years ago
#21131 closed enhancement (fixed)
Interpolation, Minimal Vanishing Polynomial and Multi Point Evaluation of Skew Polynomials
Reported by: | arpitdm | Owned by: | |
---|---|---|---|
Priority: | major | Milestone: | sage-7.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, GitHub, GitLab) | 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 6 years ago by
- Branch set to u/arpitdm/mvp_mpe
comment:2 Changed 6 years ago by
- Commit set to c568c420198a1013d05eddd134f2fabb3bd8cdac
comment:3 Changed 6 years 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 follow-up: ↓ 6 Changed 6 years 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 follow-up: ↓ 16 Changed 6 years 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" non-finite 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 1-point case you were doing quo-rem 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 non-finite 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 linear-independence 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 multi-point 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 6 years 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 multi-point evaluation is fastest, then we should just replace the method with the one-liner 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 6 years ago by
- Branch changed from u/arpitdm/mvp_mpe to u/jsrn/mvp_mpe
comment:8 Changed 6 years 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 multi-point evaluation
|
b683f55 | divide-by-two potential error. Improved an error message
|
012c9d3 | rm non-general sanity check
|
9064e0e | Fix division bug in interpolation
|
6618340 | Fix divide-by-two potential problem in interpolation
|
comment:9 Changed 6 years ago by
Here are some mathematical comments:
- I would say (but I'm not sure) that the benefit of using a divide-and-conquer method for multi-point 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 divide-and-conquer method is only interesting when fast multiplication and half-gcd are available (but these should be hidden at some point in Wachter-Zeh 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 follow-up: ↓ 11 Changed 6 years 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 construct-0 bug
|
ca15975 | Fix doctest
|
0635ec9 | Merge branch '13215_skew_polynomials' into 21131_interpolation_skew_poly
|
comment:11 in reply to: ↑ 10 Changed 6 years ago by
Merged the newest 13215 into this branch.
comment:12 Changed 6 years 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 divide-and-conquer method for multi-point 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 divide-and-conquer method is only interesting when fast multiplication and half-gcd are available (but these should be hidden at some point in Wachter-Zeh 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 & conquer-version 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 coefficient-growth, 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 follow-up: ↓ 14 Changed 6 years 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 ; follow-up: ↓ 15 Changed 6 years 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 6 years 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 multi-point evaluation is there; and C) it is good to advertise multi-point 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 ; follow-up: ↓ 17 Changed 6 years 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 6 years 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 bullet-proof 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 try-except 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 6 years ago by
- Branch changed from u/jsrn/mvp_mpe to u/arpitdm/mvp_mpe
comment:19 Changed 6 years 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 6 years ago by
- Status changed from new to needs_review
comment:21 Changed 6 years ago by
- Status changed from needs_review to needs_work
The three methods are mutually recursive, so you're still for instance re-creating the same polynomial ring over and over again in the recursive calls. You should make three top-level private functions, without the checks, and all the mutually recursive calls will call those top-level private functions, taking the correct target ring as an input.
The creation of the fraction-field skew-poly-ring should be a helper function since it's duplicated code. But make it a private function.
Best, Johan
comment:22 Changed 6 years 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 follow-up: ↓ 24 Changed 6 years 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 6 years 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 6 years 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 6 years ago by
- Status changed from needs_work to needs_review
comment:27 Changed 6 years 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 6 years 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 6 years ago by
Hi David,
I've made the changes you proposed. Opening it for review again.
Best, Arpit.
comment:30 Changed 6 years ago by
- Branch changed from u/arpitdm/mvp_mpe to u/jsrn/mvp_mpe
comment:31 follow-up: ↓ 33 Changed 6 years 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 6 years ago by
- Keywords sd75 added
comment:33 in reply to: ↑ 31 Changed 6 years 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 6 years ago by
Agreed, let's go ahead with it!
comment:35 Changed 6 years ago by
- Branch changed from u/jsrn/mvp_mpe to u/arpitdm/mvp_mpe
comment:36 Changed 6 years ago by
- Commit changed from 4bb82b8c5ffff8e442ef5aad6a224b590c3f9331 to 0572e852c1be546810e3f53ef08b569cb592415a
comment:37 Changed 6 years ago by
- Dependencies changed from #21088, #13215 to #13215
comment:38 Changed 6 years ago by
- Status changed from needs_work to needs_review
comment:39 Changed 6 years 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 6 years ago by
- Branch changed from u/arpitdm/mvp_mpe to u/jsrn/mvp_mpe
comment:41 Changed 6 years 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 5 years ago by
- Status changed from needs_review to positive_review
comment:43 Changed 5 years ago by
- Milestone changed from sage-7.3 to sage-7.6
- Reviewers set to Johan Rosenkilde
comment:44 Changed 5 years ago by
Thanks
comment:45 Changed 5 years 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