Lagrange interpolation speed problem

[Popeyef5]

I think Lagrange interpolation using this library is way slower than using vanilla python (+ numpy). Below is the code I used to test this. Sorry if it's not too clean, I just grabbed a few minutes to repurpose some code from another project and put this together.

import time
import galois
import numpy as np

# BN128 prime
prime = 21888242871839275222246405745257275088696311157297823662689037894645226208583

# Pseudo random points used to interpolate
points = [
  (1, 3),
  (76, -1),
  (301, -1),
  (306, -3),
  (316, 2),
  (346, 1),
  (77, -1),
  (302, -1),
  (307, -3),
  (317, 2),
  (347, 1),
  (78, -1),
  (303, -1),
  (308, -3),
  (318, 2),
  (348, 1),
  (79, -1),
  (304, -1),
  (309, -3),
  (319, 2),
  (349, 1),
  (80, -1),
  (305, -1),
  (310, -3),
  (320, 2),
  (350, 1),
  (91, -1),
  (92, -1),
  (93, -1),
  (94, -1),
  (95, -1),
  (96, -1),
  (97, -1),
  (98, -1),
  (99, -1),
  (100, -1),
  (101, -1),
  (102, -1),
  (103, -1),
  (104, -1),
  (105, -1),
  (106, -1),
  (107, -1),
  (108, -1),
  (109, -1),
  (110, -1),
  (111, -1),
  (112, -1),
  (113, -1),
  (114, -1),
  (115, -1),
  (116, -1),
  (117, -1),
  (118, -1),
  (119, -1),
  (120, -1),
  (121, -1),
  (122, -1),
  (123, -1),
  (124, -1),
  (125, -1),
  (126, -1),
  (127, -1),
  (128, -1),
  (129, -1),
  (130, -1),
  (131, -1),
  (132, -1),
  (133, -1),
  (134, -1),
  (135, -1),
  (136, -1),
  (137, -1),
  (138, -1),
  (139, -1),
  (140, -1),
  (216, -1),
  (217, -1),
  (218, -1),
  (219, -1),
  (220, -1),
  (221, -1),
  (222, -1),
  (223, -1),
  (224, -1),
  (225, -1),
  (226, -1),
  (227, -1),
  (228, -1),
  (229, -1),
  (230, -1),
  (231, -1),
  (232, -1),
  (233, -1),
  (234, -1),
  (235, -1),
  (608, 16931903655476233095338254632264576238984860566592271344419174450729693154521),
  (609, 17146236214239816754147166163180990870230477648858748237220746264352397259281),
  (610, 6559847902765392927165377899101119595508332203342503531171719469473018943799),
  (611, 21263175092961191970729563609312792029626806314782521449370169321138435264489),
  (612, 8234573390286426051293085163061879507942235256258552856130923126023545753989),
  (613, 140024849171158027271291432481803661182650119307416967156863543882698526323),
  (614, 7225711375668528854657280807873121500722446116918743736594268363661108843873),
  (615, 14723595833217001934180066950501220341300922529823683265822680411824719048037),
  (616, 17067098378478504268539350534165895549220481996620100462708443112054585188817),
  (617, 6766451000083880288615325051605873635067013178674576210279023698938040578728),
  (618, 1850732451355444347101073140876364882826600284544979808419113817064716321813),
  (619, 9198758778983429378277385792637172896827499204255370369014849339389295241212),
  (620, 17527239903732239506881523521452399783834404526665535092768727679589907923981),
  (621, 4516550744097282765485311396867568272926768521268796161835140697953767235794),
  (622, 18390759285162489832348762787674384213745708518363481133293452477688646027132),
  (623, 14946722811212092904408944554011516250513957316282253965864455484234320488032),
  (624, 2834100621990113649740659462733521737017109282759538927809560664122282800314),
  (625, 5730640341066681336938066801956599788878104824047087731829943707659143584451),
  (626, 19100982549115761248494229384956645249034500513053223712558645898875249015294),
  (627, 20369638168768372250489646661017913465266071414690792104334557208159576130434),
  (628, 20570003701327718970829176693005700545212496026253209411345832907513067856338),
  (629, 4590161137004042632834207918652608416144186927362851968127491124759464257185),
  (630, 2196402916844031717525796063127314220830237169370629160222591752300843961674),
  (631, 2714661209920971654978942933076269549752802731184800538106681507285480659389),
  (632, 15136658254864176260192369899095401215484724479224947479759449976197074538155),
  (633, 14635834768955082886162542397479526273049138321878538905391805350837402302386),
  (634, 7910086853870684275594547992537154110394401980647780348161063986329776307621),
  (635, 15760576805200608897429478857627237758543741709780112075597302066039617415507),
  (636, 10647079538229532768090705373375971300598894011394929124617870880165404110013),
  (637, 9408088493111174442341032476524966600934814475381223071193031045194555008694),
  (638, 21362817631555286659756463208394807409761942157562773333130845140058930117760),
  (639, 1104109313804163895748716611811155089644991462892845136974224263303705136424),
  (640, 4598446496303676576331478262933730990329162005038805292742250332402381128324),
  (641, 8440418641303441336322822075286272671602015037111197583710766540972603311533),
  (642, 21482921306416717094225480031586073092993163968959209350990573920442299254519),
  (643, 5870221049464364229418443331539371655251070467668733479757269857462265311717),
  (644, 15300631826321037085121844685188879990355259901516667280424629887963739006535),
  (645, 20962675858791601447296960762315408207917520295167189611912440961527981660356),
  (646, 13988000419549441272888758369470811322169194748619918486642004046013254163174),
  (647, 21517464888385365414827043294364421760267577288372678431898476002025499448215),
  (648, 13288408766182654306479214875946903779157695664956790686495628961299430647654),
  (649, 20551973426895345046727441788422746067862489692885101934942928206847259996833),
  (650, 15844592789846983516515605889332715299792035821851924404054534596553589662370),
  (651, 14356933985599897277435099838832784301895366769815120536363043541166857033696),
  (652, 19002233373703875276358772177352739414830025344107136134224495635426514115252),
  (653, 14255638036699992960194674852661068806619126771459203305923915813697224358668),
  (654, 11461177725810942383933358476266419791702942239605170881668672622537456892310),
  (655, 15776917692974511402844607320529306346751418503610308403495205339356369904511),
  (656, 1953517658244721736255782956839611321645760203293058340973950286053044238978),
  (657, 13255711256976914018935937183016103094846783968230662139411391760316093897042),
  (658, 20153490026823100497557517417463746391764585693453798230213741259837998220915),
  (659, 11110382764805604956428288339123450246646760315385035089552463754745526795823),
  (660, 5391756991387121063197672718954829885623112343713376364996501413933051789248),
  (661, 5819371711881274789746218329396245077878312265785529744530415543729805413362),
  (662, 2557530292066403870924558973000736629675413345790756114076289155545154781963),
  (663, 20635720600940586329947927933760961842593182299926198490243114807053050739133),
  (664, 11884345252574085004330726311169768739050589462261388467532851348861421510306),
  (665, 7737417395618651376860372770952510967233179290460018011299083543680917019447),
  (666, 20227589403005159750902726044134789771669735411658551788873877762913133912490),
  (667, 2355819957807718517086546512155519380874059619930715352340520770638560549757),
  (668, 10958386503472692906597320509006450119401889122894084383650892096107402031696),
  (669, 1206780393593607081950759425975827543865021499456811172871607162179880991540),
  (670, 730753856774702382632908487600033870391982585749293833625010459048093010285),
  (671, 14293604770338907046413411468626662208600242092671070123891535270955916765716),
  (669, 0),
  (672, 16133709183648322543378792128670975459028583618420154089279244521066485287875),
  (675, 4045606801120378565512433914548501700467375859878185669837874391723571030524),
  (678, 186786316964044854555936436040180060239542074314500953092064910258057529519),
  (681, 20538094816961779044884033286127102647716226139795445777577570481064543046604),
  (684, 14162962649353112893335591059731619607351340362491753470082472140947663468272),
  (687, 3057383713887837322108513281944228331119632360528430802685646965761012221068),
  (690, 18497273444356915005100710005829911929324673997925491575326210181378368950950),
  (693, 15159645247407927847318309507596387891478538397635349081470686095546674754705),
  (696, 19702767697178831579241120482849049896840903719963267640084666869097128561184),
  (699, 21877131240975466850089711900934091622209528027073261143442762323318475837119),
  (702, 16583195271933774326399785501207785596854313053843958191198593454465154941702),
  (705, 11071289357782909654652063024437769403281078796984653231284421696216741121055),
  (708, 16773447414236905978872448950160598758263506629965242833346132267014073429373),
  (711, 429156852780434262676968308986247491491646336161290007561714672072738784549),
  (714, 8918230268796444938584693958974550234146316020491855744283571923762800970029),
  (717, 2280855086573145108086819699182727626459299096481588785779942335289319677880),
  (720, 5324022079938615267519189125159971558124012624249939627524059075831014716095),
  (723, 11107791706652849606644206765799059487378070815203907732999957194883295373292),
  (726, 7798478561246498181896849626672511877858000248104970660084608967050847540174),
  (729, 12586086184133855546331518589812856687238599280703476690207720257440746302009),
  (732, 6375774933183616626549029331170309220425052095863917524786643663928295360936),
  (735, 7379391481459794714955741061369209144174113628619985724966001815900630231223),
  (738, 6200702115948224899682204872150287519271614978242424529935529781909163502314),
  (741, 8903429865966409908293341623415702434925811340236935457546462399950938053266),
  (744, 4685698575505997742946913376385460960317548251357177171956421445861563509918),
  (747, 1979621720413076416433847421036427454624968427661579056095281495088556719323),
  (750, 11112700063426566158507161987370954969108116489341890968239384090199196191342),
  (753, 10789286966954434008460508373088674550040981599072082468580906554422365085052),
  (756, 3548681439038749671599099154637075839443699195000989731363961813776783781746),
  (759, 19253556022548416678852784426808572323040193943371388592757276106095149472178),
  (762, 19571106126238596181137148773100703191839149093503500918436851867201293058309),
  (765, 18835920608752168806035641158626828076655707998771337841363451442794525483350),
  (768, 9200152588719482716121490337539688761158649617574191055322830958432630885154),
  (771, 17087240056941382771605692866106227418836151397771918904026405586857337546313),
  (774, 3090296829639517570789795183431541852877314568957917599455353906267639470718),
  (777, 8528083195401188147378045415827956364973701062573164844655094846325537474213),
  (780, 18823677258989884515063323695768957784611474101443603979000250490629666492765),
  (783, 11377213024335459128280537060037505502340562003593812792751219248422173803119),
  (786, 12778109452474915351097274350343428244779509237792164105194512232201668088026),
  (789, 18929364173536716181271242285049437131422289816769861048803769330366249934411),
  (792, 804705889386908484460233906370859755984380946980046565984325515767014086269),
  (795, 8693229094551332817235681954132021607640482115692747651272333767243690961473),
  (798, 16839175748539200744011407356854202983660131032343293126269944827364629473521),
  (801, 9927482243707690647110940499487406378272220126157897638422020373322339006230),
  (804, 8754098645202775752964432409528228840473563451903029566731224047829829328758),
  (807, 18677783724380014553670092837321046842873677549744895276582759686699132672720),
  (810, 1022729485973043432774432640009671304156105096100141819151419572960838724478),
  (813, 861020481525123734469894521512209589083840080913334291588388296601490574981),
  (816, 15403338323865179458853606256746007053680183002963608730619470760025763996161),
  (819, 9577711637182991494004369593317413795851017562011530125394342660930391781057),
  (822, 12423459649190387394211902610953039699778838239378975303472618172368442796177),
  (825, 12019459759174777161483654794850158767769540552974523730296464950860224350771),
  (828, 6891579403412065807700730070102515838052600664961585588907627212937237925113),
  (831, 17151652880752623048694637785367185176797124158492255144889735212107169370765),
  (834, 1767217940876143551933815695100867176763679419480788964761836556629928211667),
  (837, 4346113949279595622029914426135773883087645469015648113390978406692481192966),
  (840, 14857887677542248043416592787290475095732537922292075709873035019518468699315),
  (843, 6257862018539659392966471042776268825669355975888962694625759909109564670726),
  (846, 21129041182910710787148377343801009976757369485137830147274868848291718412175),
  (849, 18355461139232681416100453168981970215886245258068959310155593850688015252943),
  (852, 1910651628904534112343965659424742861243187389154571574355347589400952741333),
  (855, 3070082047638309603519269790808658519320716998052311241645637545833751299695),
  (858, 14414009007687342025526645003307639786191886886413750648631138442071909631647)
]


def mod_inv(x, p):
  x = x % p
  t = 0
  nt = 1
  r = p
  nr = x
  while nr != 0:
    q = r // nr
    t, nt = nt, t - q * nt
    r, nr = nr, r - q * nr
  if r != 1:
    raise ValeError("number has no inverse")
  return t % p


def polydivmod(u, v, p):
  m = len(u)
  n = len(v)
  s = mod_inv(v[0], p)
  q = np.array([0]*max(m-n+1, 1), dtype=object)
  r = np.copy(u)

  for i in range(m - n +1):
    d = s * r[i] % p
    q[i] = d
    r[i:i+n] = np.mod(r[i:i+n] - v * d, p)

  first = 0
  for i in range(len(r)-1):
    if r[i] != 0:
      break
    first += 1

  r = r[first:]

  return q, r


class fpoly1d:

  def __init__(self, c,  prime=prime):
    self.prime = prime

    if isinstance(c, fpoly1d):
      self.coeffs = c.coeffs
      return

    c = np.atleast_1d(c).astype(object)
    c = np.trim_zeros(c, 'f')
    self.coeffs = np.mod(c, self.prime)

  def __mul__(self, other):
    if isinstance(other, fpoly1d):
      return fpoly1d(np.convolve(self.coeffs, other.coeffs))
    else:
      return fpoly1d(self.coeffs * other)

  def __rmul__(self, other):
    if np.isscalar(other):
      return fpoly1d(other * self.coeffs)
    else:
      other = fpoly1d(other)
      return fpoly1d(np.convolve(self.coeffs, other.coeffs))

  def __add__(self, other):
    other = fpoly1d(other)
    return fpoly1d(np.polyadd(self.coeffs, other.coeffs))

  def __radd__(self, other):
    other = fpoly1d(other)
    return fpoly1d(np.polyadd(self.coeffs, other.coeffs))

  def __div__(self, other):
    if np.isscalar(other):
      other_inv = mod_inv(other, self.prime)
      return fpoly1d(self.coeffs * other_inv)
    q, r = polydivmod(self.coeffs, other.coeffs, self.prime)
    return fpoly1d(q), fpoly1d(r)

  __truediv__ = __div__

  def __rdiv__(self, other):
    if np.isscalar(other):
      other_inv = mod_inv(other. self.prime)
      return fpoly1d(self.coeffs * other_inv)
    q, r = polydivmod(self.coeffs, other.coeffs, self.prime)
    return fpoly1d(q), fpoly1d(r)


def vanilla_interpolate(x, y):
  # Emulating SciPy's Lagrange implementation
  poly = fpoly1d(0)
  for i in np.flatnonzero(y):
    pt = fpoly1d(y[i])
    for j in range(len(x)):
      if i == j:
        continue
      fac = x[i] - x[j]
      pt *= fpoly1d([1, -x[j]])/fac
    poly += pt
  return poly


def main():
  # Settings
  q = 10
  r = 360

  # Galois
  t0 = time.monotonic()
  
  GF = galois.GF(prime)
  
  t1 = time.monotonic()
  print(f"Finished creating GF(p). Took {t1-t0}s.")

  x = np.arange(1, r+1, dtype=object).view(GF)
  y = np.zeros(r, dtype=object).view(GF)
  for px, py in points[:q]:
    y[px-1] = py % prime


  p_g = galois.lagrange_poly(x, y)
  
  t2 = time.monotonic()
  print(f"Finished interpolating with Galois. Took {t2-t1}s.")

  # Vanilla python
  x = np.arange(1, r+1, dtype=object)
  y = np.zeros(r, dtype=object)
  for px, py in points[:q]:
    y[px-1] = py % prime
  
  p_v = vanilla_interpolate(x, y)

  t3 = time.monotonic()
  print(f"Finished interpolating with vanilla python. Took {t3-t2}s. Factor: {(t2-t1)/(t3-t2)}")

  # Consistency
  print(f"Consistency check (should be 0): {min(abs(p_g.coefficients()-p_v.coeffs.view(GF)))}")


if __name__ == "__main__":
  main()

This worked with multiple combinatios of (q, r).
(1, 3) had ~20x difference.
(10, 360) 150x difference.
(224, 1110) took Galois almost two hours.

Let me know if you need anything else.

[mhostetter]

I will start looking into this.

One tip you may not be aware of. The primary time hog in generating the large GF(p) field is factoring p - 1 (which is needed for finding the primitive root). This is a current limitation of the factorization algorithms, see #187 and #103. The current workaround is to cache the primitive element after a single construction and use verify=False.

Before

In [1]: import galois

In [2]: %time GF = galois.GF(21888242871839275222246405745257275088696311157297823662689037894645226208583)
CPU times: user 1min 6s, sys: 0 ns, total: 1min 6s
Wall time: 1min 6s

In [3]: print(GF.properties)
Galois Field:
  name: GF(21888242871839275222246405745257275088696311157297823662689037894645226208583)
  characteristic: 21888242871839275222246405745257275088696311157297823662689037894645226208583
  degree: 1
  order: 21888242871839275222246405745257275088696311157297823662689037894645226208583
  irreducible_poly: x + 21888242871839275222246405745257275088696311157297823662689037894645226208580
  is_primitive_poly: True
  primitive_element: 3

After

In [1]: import galois

In [2]: %time GF = galois.GF(21888242871839275222246405745257275088696311157297823662689037894645226208583, primitive_element=3, verify=False)
CPU times: user 578 µs, sys: 872 µs, total: 1.45 ms
Wall time: 1.45 ms

In [3]: print(GF.properties)
Galois Field:
  name: GF(21888242871839275222246405745257275088696311157297823662689037894645226208583)
  characteristic: 21888242871839275222246405745257275088696311157297823662689037894645226208583
  degree: 1
  order: 21888242871839275222246405745257275088696311157297823662689037894645226208583
  irreducible_poly: x + 21888242871839275222246405745257275088696311157297823662689037894645226208580
  is_primitive_poly: True
  primitive_element: 3

[Popeyef5]

That's definitely helpful, thanks!

Though, just to clarify, the time spent generating GF wasn't taken into the performance comparisons. It was mostly lagrange_poly vs vanilla_interpolate