Skip to main content

symbolic - Improving the performance of a package for working with rational functions



As Mathematica gets slow for large symbolic calculations, the cost of putting terms over a common denominator (Together), in particular, gets too high. It occurred to me that, if one has a small number of variables, it should be much faster if we represent our expressions by arrays of their coefficients.


I have implemented a small package that does this with an improvement of a factor of 10 in both speed and memory (in putting Together a sum of rational terms, more or less independent of the size of the sum), but it seems to me that I should be able to do better. I post my package here with the hope that someone might point out glaring performance bottlenecks.


First, the public initializations, with their usability comments.


AA::usage = "DATA STRUCTURE AA[num,den] that stores polynomials as their tensor coefficients of a basis, multiplication is the tensor product"
PolytoAA::usage = "PolytoAA[poly_,vars_List] returns an AA[num_SparseArray,1]"
AAtoPoly::usage = "AAtoPoly[aa_AA,vars_] gives back a rational function of the polynomials num/den"

The actual code begins here, PolProd takes care of multiplying polynomials


SetAttributes[PolProd, Orderless]
PolProd[i__Integer] := Times @ i;

PolProd[sa__SparseArray] := Outer[Times, sa];
  (*multiplying 2 polynomials is the tensor product in their vector space*)
PolProd[i__Integer, sa__SparseArray] := Times[i] Outer[Times, sa];

(*It would be cleaner and more object oriented to call two instances of polprod 
  on i and sa but it would force more pattern matching so it would be less efficient*)

This code performs the arithmetic operations on the AA[num,den] data structure


AA /: i_Integer*AA[num_, den_] := AA[i*num, den]
 (*absorb all integer into AA, this means that the sum rule has to do less 

   pattern matching and is more efficient, overall more efficient?*)
AA /: AA[a_, a_] := 1
AA /: AA[0, _] := 0
AA /: AA[num1_, den1_]*AA[num2_, den2_] := 
  AA @@ (PolProd[#1, #2]& @@@ {{num1, num2}, {den1, den2}})
AA /: AA[num1_, den_] + AA[num2_, den_] := 
  AA[num1 + num2, den]
  (*can't check for a SparseArray cos it might be an Integer, 
    by dimensional analysis we should get two arrays here though*)
AA /: AA[num1_, den1_] + AA[num2_, den2_] := 

  AA[PolProd[num1, den2] + PolProd[num2, den1], PolProd[den1, den2]]
  (*together rule, this a recursive implementation that doesn't treat the 
    sum of many AA all at once, this is intentional cause it minimises the 
    number of array outer products that are calculated and I believe it's 
    the more efficient but maybe not*)
AA /: aa_AA^n_Integer := 
  If[n < 0, Reverse @ #, #]& @ (PolProd[Sequence @@ ConstantArray[#, Abs @ n]]& /@ aa)
  (*this takes care of both efficiency of recursively pattern matching the 
    product rule on large n and dividing by AA*)
(*NOTICE THE WEIRD SYNTAX/BUG If[n<0, Reverse @ #, #]& @ THERE ARE TWO INTANCES OF @*)


I also post the the functions that go between AA and the polynomials for the sake of clarity and completeness


AAtoPoly[aa_AA,vars_] := #1/#2& @@ (aa /. l_SparseArray :> 
  Dot[l, Sequence @@ ConstantArray[vars, Depth[l] - 1 (*gives rank of the tensor*)]])
PolytoAA[expr_^n_, vars_List] := PolytoAA[expr, vars]^n
PolytoAA[expr : (_Times | _Plus), vars_List] := PolytoAA[#, vars]& /@ expr
PolytoAA[expr_,vars_List] /; PolynomialQ[expr, vars] := 
  AA[Last[CoefficientArrays[expr, vars]](*num*), 1 (*den*)]

BENCHMARKS



I define some random, small polynomials


vars = a /@ Range[10]
pols = Table[
   a[RandomInteger[{1, 10}]] + a[RandomInteger[{1, 10}]], {20}] // 
  DeleteDuplicates
% // Length


{a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8], a[9], a[10]}


{a[1] + a[7], a[4] + a[6], a[2] + a[4], a[5] + a[9], 2 a[8], a[5] + a[6], a[2] + a[9], a[4] + a[7], a[4] + a[9], a[7] + a[9], a[3] + a[7], a[2] + a[6], a[1] + a[9], a[2] + a[3], a[7] + a[10], 2 a[2], a[3] + a[9]}



17



First I put them Together using my code


aa = PolytoAA[#, vars] & /@ pols;
togaa = Sum[aa[[i]]/aa[[i + 1]], {i, Length@pols - 1}] // 
  AbsoluteTiming
% // ByteCount


{0.8960513,AA[SparseArray[<274432>,{10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10}],SparseArray[<16384>,{10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10}]]}



18614128



and now without my implementation


togpols = (Sum[pols[[i]]/pols[[i + 1]], {i, Length@pols - 1}] // 
    Together) // AbsoluteTiming//
% // ByteCount


{8.2054693,...}


34567480




I check that both give the same result


togpols2 = AAtoPoly[togaa // Last, vars](*//Simplify*);
Last[togpols] - togpols2 // Together


0





Comments

Post a Comment

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...
Post a Comment
Read more

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...
Post a Comment
Read more