Interpolation, Minimal Vanishing Polynomial and Multi Point Evaluation of Skew Polynomials

[arpitdm]

Support for computing:

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

in Skew Polynomial Rings.

[git]

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

​c568c42

added interpolation polynomial method

[git]

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

​282cfb1

added check for linear independence of evalpts in MVP

[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.

[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

[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

[jsrn]

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

[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).

[git]

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

[jsrn]

Merged the newest 13215 into this branch.

[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

[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.

[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.

[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

[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?

[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

[arpitdm]

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:

​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.

[jsrn]

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

[git]

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.

[arpitdm]

I created three new private functions namely _gen_one_sigma, _create_mvp and _interpolate. Are these names alright?

Best, Arpit.

[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

[git]

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

​b523806

changed names of private methods

[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

[git]

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

[arpitdm]

Hi David,

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

Best, Arpit.

[jsrn]

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:

​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

[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.

[jsrn]

Agreed, let's go ahead with it!

[arpitdm]

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:

​0572e85

Removed dependency on the SkewPolynomial_finite_field class.

[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

[jsrn]

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

---

New commits:

​4a9e4b0

Doc markup fixes

[jsrn]

Thanks