P-S rule as a general poly recursion call

[jarlebring]

If we only count computation cost in terms of :mul-operations, the general polynomial recursion implemented in gen_general_poly_recursion is a framework including all polynomial based algorithms. Therefore, the P-S rule can also be reformulated in this form. It would be useful to have a variant of the P-S rule graph generator, where we generate the graph by a call to gen_general_poly_recursion. It can be useful, since then we can use the P-S rule as a starting value for optimization of the graph.

A proof-of-concept is in bbc_exp.jl: For 2 matrix-matrix multiplies it carries out the P-S rule with k=4, s=2, \nu=2.

The current P-S rule implementation can remain where it is, since it is more efficient in terms of memory and if we also count :lincomb's.

[eringh]

I did an initial implementation of this. Now it generates the correct polynomial, but the number of multiplications seems to be increased. (Perhaps it is multiplications with I, when we do linear combintations.) Can anyone take a look and see if it seems reasonable?

[jarlebring]

Yes. Your first elements in x are okay (you are constructing the monomials), after that you have a couple of ([XXXXX],[1,0,....]) which are multiplications with identity and can be avoided. The conversion of PS to this form needs to be documented / explained, also in the manuscript so we anyway need to work out an example by hand. This is quite tricky to get right.

More precisely these lines generate x[s],... x[s+v-1] which are formed from identity multiplications https://github.com/jarlebring/GraphMatFun.jl/blob/91082164fad7c348af8c3a85802a661ba1a7571f/src/generators/ps_poly.jl#L170-L177

These can be merged with the later x-values. Roughly speaking, wherever one of the x[s],... x[s+v-1] is used, you should instead add the coefficients that are currently in x[s][1],... x[s+v-1][1].

I would first store the coefficients currently in x[s][1],... x[s+v-1][1] as a temporary variable xt:

     xt=Vector{Vector{T}}(undef,v);
     for i = 0:(v-1)
        lower_idx = i*s+1
        upper_idx = lower_idx+s-1
        idx = lower_idx:upper_idx
        c = view(a, idx)
        xt[i+1] = vcat(c,zeros(T,i+1))
    end

In this way, for instance this line https://github.com/jarlebring/GraphMatFun.jl/blob/91082164fad7c348af8c3a85802a661ba1a7571f/src/generators/ps_poly.jl#L185 can be converted to something like

  xt_1=vcat(c,zeros(diff),zeros(T,i-s+1));
  xt_2=vcat(one(T),zeros(T,i))
  x[i]=(xt_1[1:s],xt_2[1:s]);
  for j=s+1:size(xt_1,1)
      x[j][1] += xt_1[j]*xt[j-s];
      x[j][2] += xt_2[j]*xt[j-s];
  end

Hopefully the general idea can be understood from the above. I can describe more later.

[eringh]

Thanks! Improved in db6527b0ca41965855aa0ea3ac27a1a651c6d7cd, but still some more to do.

[eringh]

To easier see what is happening, one can paste the code

for i = 1:length(x)
    println( x[i][1], "     ", x[i][2])
end
println(z)

just before the return-statement in gen_ps_poly_recursive, and then run

import Base.zero
zero(::Type{Any}) = 0
import Base.one
one(::Type{Any}) = 1
coeff = Vector{Any}([:C0, :C1, :C2, :C3, :C4, :C5, :C6, :C7, :C8, :C9, :C10, :C11, :C12])
(graph,cref) = gen_ps_recursive(coeff);

The output is

Any[0, 1]     Any[0, 1]
Any[0, 0, 1]     Any[0, 1, 0]
Any[0, 0, 0, 1]     Any[0, 1, 0, 0]
Any[0, 0, 0, 0, 1]     Any[:C8, :C9, :C10, :C11, :C12]
Any[0, 0, 0, 0, 1, 0]     Any[:C4, :C5, :C6, :C7, 0, 1]
Any[:C0, :C1, :C2, :C3, 0, 0, 1] 

This is the case of s dividing k, i.e., the highest degree achievable for a fixed number of multiplications. For lower degrees, the vector Any[:C8, :C9, :C10, :C11, :C12] will start to fill out with zeros.

[jarlebring]

Cool! The pattern doesn't seem so complicated. Maybe the generation code can now be completely rewritten in a simpler way? A clear formulation of PS in terms of the recursive form is needed for the manuscript anyway.

[eringh]

I cannot find a simpler formulation. Even the original PS-implementation has some special handling for when the degree of the polynomial, k, is divisible by s=ceil(sqrt(k)).