Opened 12 months ago

Closed 5 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: 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) Commit: 4a9e4b097681ffa0534918c01b4d31b41b38cfbf
Dependencies: #13215 Stopgaps:

Description

Support for computing:

  1. Interpolation Polynomial
  2. Minimal Vanishing Polynomial
  3. Multi Point Evaluation

in Skew Polynomial Rings.

Change History (45)

comment:1 Changed 12 months ago by arpitdm

  • Branch set to u/arpitdm/mvp_mpe

comment:2 Changed 12 months ago by git

  • Commit set to c568c420198a1013d05eddd134f2fabb3bd8cdac

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

c568c42added interpolation polynomial method

comment:3 Changed 12 months ago by git

  • Commit changed from c568c420198a1013d05eddd134f2fabb3bd8cdac to 282cfb1bf5a92744918b56d0190f6789295ca532

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

282cfb1added check for linear independence of evalpts in MVP

comment:4 follow-up: Changed 12 months ago by arpitdm

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: Changed 12 months ago by jsrn

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. if S = R[x; sigma] for some ring R, then the minimal vanishing polynomial is defined over Q[x; sigma] where Q is the fraction field of R. This is clear from the line x - (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 for multi_point_evaluation and interpolation.
  • I removed a check from interpolation_polynomial. That check did not make any sense in most cases (the characteristic of ZZ[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 of interpolation_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 12 months ago by jsrn

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 12 months ago by jsrn

  • Branch changed from u/arpitdm/mvp_mpe to u/jsrn/mvp_mpe

comment:8 Changed 12 months ago by jsrn

  • Commit changed from 282cfb1bf5a92744918b56d0190f6789295ca532 to 6618340ece629db87b92d77b0c7ae0e482480791

forgot to push my changes.


New commits:

3166c6dMerge branch 'u/arpitdm/mvp_mpe' of git://trac.sagemath.org/sage into 21131_interpolation_skew_poly
7806d86Fixed bug in interpolation
e1cce49Improved doc test for multi-point evaluation
b683f55divide-by-two potential error. Improved an error message
012c9d3rm non-general sanity check
9064e0eFix division bug in interpolation
6618340Fix divide-by-two potential problem in interpolation

comment:9 Changed 12 months ago by 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
  • 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: Changed 12 months ago by git

  • Commit changed from 6618340ece629db87b92d77b0c7ae0e482480791 to 0635ec993a6434362e13c6d2b617670f11c6a6c8

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

362fff7Merge branch 'u/arpitdm/skew_polynomials' of git://trac.sagemath.org/sage into 13215_skew_polynomials
4a33c5fmv skew_polynomial_element in module_list
a921040Changed unfortunate parameter name
f43767aSome minor doc and error message modifications
7c5c511Modelled SkewPolynomialBaseringInjection more closely on PolynomialBaseringInjection
d6fd904Fix Cython warning
7134b7eRemoved _call_with_args
291c28cFixed construct-0 bug
ca15975Fix doctest
0635ec9Merge branch '13215_skew_polynomials' into 21131_interpolation_skew_poly

comment:11 in reply to: ↑ 10 Changed 12 months ago by jsrn

Merged the newest 13215 into this branch.

comment:12 Changed 12 months ago by jsrn

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: Changed 12 months ago by jsrn

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: Changed 12 months ago by arpitdm

Replying to jsrn:

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.

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 12 months ago by jsrn

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: Changed 12 months ago by arpitdm

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 12 months ago by jsrn

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 12 months ago by arpitdm

  • Branch changed from u/jsrn/mvp_mpe to u/arpitdm/mvp_mpe

comment:19 Changed 12 months ago by arpitdm

  • Commit changed from 0635ec993a6434362e13c6d2b617670f11c6a6c8 to 8f42a338565e35ce2210216c9ec9e7a240f97fe9

I've made the following changes based on discussions above:

  1. created twist map over the fraction field. fixed methods to accommodate fraction fields.
  2. 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.
  3. changed method MPE to the trivial repeated calling of the evaluate method and moved this to class Skew_Polynomial.

I'm opening the ticket for review now.

Best, Arpit.


Last 10 new commits:

c034df8removed unused private methods and imports
ac109c8converted to standard copyright template
5137ef9removed deprecated method getslice
b1e42c3Merge commit 'aca3398a81b3688e1f2e69c5910b8214d13be925' of git://trac.sagemath.org/sage into skew_polynomials
f563abamod and floordiv modified to use coercion model. replaced type with isinstance.
26e4689modified copyright banner
6c86537merged updated code from 13215. fixed interpolation to accomodate fraction fields. updated documentation.
50f9bddrefactored mvp and interpolation to have inner functions so that checks on the arguments are made only at the top level.
dd667edchanged MPE to the trivial repeated calling of the evaluation function. added a todo block to the documentation.
8f42a33made MPE a method on skew polynomials.

comment:20 Changed 12 months ago by arpitdm

  • Status changed from new to needs_review

comment:21 Changed 12 months ago by jsrn

  • 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 12 months ago by git

  • Commit changed from 8f42a338565e35ce2210216c9ec9e7a240f97fe9 to 34bebbc6eb951d7e459ee2ee9988461d89af27d3

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

f25f886created top level private function for MVP
88196bfsome more fixes to the documentation.
b9f8239merging latest changes
34bebbccreated top level functions for interpolation and fraction field skew polynomial ring. added documentation and indirect doctests.

comment:23 follow-up: Changed 12 months ago by arpitdm

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 12 months ago by jsrn

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 12 months ago by git

  • Commit changed from 34bebbc6eb951d7e459ee2ee9988461d89af27d3 to b523806a461e61036de748c7d7598044acc85892

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

b523806changed names of private methods

comment:26 Changed 12 months ago by arpitdm

  • Status changed from needs_work to needs_review

comment:27 Changed 12 months ago by dlucas

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` and interpolation_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 12 months ago by git

  • Commit changed from b523806a461e61036de748c7d7598044acc85892 to d403fae361316c8da23e388381e69148863ae48b

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

eb78e69small fixes to docstrings
2d67e0emerging updates
a2c4f06removed unused imports, signal statements. small fixes to documentation.
e435f56Merge branch 'u/arpitdm/mvp_mpe' of git://trac.sagemath.org/sage into mvp_mpe
d403faefixes to documentation

comment:29 Changed 12 months ago by arpitdm

Hi David,

I've made the changes you proposed. Opening it for review again.

Best, Arpit.

comment:30 Changed 11 months ago by jsrn

  • Branch changed from u/arpitdm/mvp_mpe to u/jsrn/mvp_mpe

comment:31 follow-up: Changed 11 months ago by jsrn

  • Commit changed from d403fae361316c8da23e388381e69148863ae48b to 4bb82b8c5ffff8e442ef5aad6a224b590c3f9331
  • Status changed from needs_review to needs_work

I made a number of modifications:

  • Renamed interpolation to lagrange_polynomial to match PolynomialRing. 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 from lagrange_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:

5041176Merge branch 'u/arpitdm/mvp_mpe' of git://trac.sagemath.org/sage into 21131_interpolation_skew
ce8b6b4mpe: modified TODO and improved doc test
256052aminimum -> minimal
2707768clarify a doctest
3a1b992Moved _-functions outside class. Simplified call conventions
0ed0179Allow minimal_vanishing_polynomial on input which is linearly dependent
db72757interpolation -> lagrange_polynomial and make call convention as PolynomialRing
115cf04Improved doc and example of _base_ring_to_fraction_field
ff371ablagrange: improve doc. Rm check param and instead discover linear dependence
4bb82b8Improve a bit the code and doc

comment:32 Changed 11 months ago by jsrn

  • Keywords sd75 added

comment:33 in reply to: ↑ 31 Changed 10 months ago by arpitdm

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 10 months ago by jsrn

Agreed, let's go ahead with it!

comment:35 Changed 10 months ago by arpitdm

  • Branch changed from u/jsrn/mvp_mpe to u/arpitdm/mvp_mpe

comment:36 Changed 10 months ago by arpitdm

  • Commit changed from 4bb82b8c5ffff8e442ef5aad6a224b590c3f9331 to 0572e852c1be546810e3f53ef08b569cb592415a

I've removed the dependency on the finite field class. Multi-point evaluation, minimum vanishing polynomial and Lagrange interpolation are now available in the original files from #13215.


New commits:

0572e85Removed dependency on the SkewPolynomial_finite_field class.

comment:37 Changed 10 months ago by arpitdm

  • Dependencies changed from #21088, #13215 to #13215

comment:38 Changed 10 months ago by arpitdm

  • Status changed from needs_work to needs_review

comment:39 Changed 9 months ago by jsrn

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 9 months ago by jsrn

  • Branch changed from u/arpitdm/mvp_mpe to u/jsrn/mvp_mpe

comment:41 Changed 9 months ago by jsrn

  • Commit changed from 0572e852c1be546810e3f53ef08b569cb592415a to 4a9e4b097681ffa0534918c01b4d31b41b38cfbf

I fixed some small doc markup issues. Otherwise the code looks good to me.


New commits:

4a9e4b0Doc markup fixes

comment:42 Changed 6 months ago by arpitdm

  • Status changed from needs_review to positive_review

comment:43 Changed 6 months ago by tscrim

  • Milestone changed from sage-7.3 to sage-7.6
  • Reviewers set to Johan Rosenkilde

comment:44 Changed 6 months ago by jsrn

Thanks

comment:45 Changed 5 months ago by vbraun

  • Branch changed from u/jsrn/mvp_mpe to 4a9e4b097681ffa0534918c01b4d31b41b38cfbf
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.