matrixfunctions / GraphMatFun.jl

Computation graphs for matrix functions
MIT License
13 stars 0 forks source link

Sastre code gen #30

Closed jarlebring closed 3 years ago

jarlebring commented 3 years ago

Although generators/sastre.jl does contain some code corresponding to the paper https://doi.org/10.1016/j.laa.2017.11.010 it is not yet complete.

I don't think the MEP-code should go into the package / paper - too complicated to make it complete. If possible we can provide the coefficients that were not reported in the paper by using the MEP-approach. This is currently done for the exponential with 6 multiplications.

Details to be specified.

eringh commented 3 years ago

I started on this by introducing tests, and I cannot get gen_sastre_basic_exp(6) to work. I have looked, and the coefficients match those in Table 7, but the output is not a Taylor approximation of the exponential.

(graph,cref) = gen_sastre_basic_exp(6);
using Polynomials
x=Polynomial("x");
coeffs(eval_graph(Compgraph(Any,graph),x))
21-element Vector{Float64}:
 2.755731922398589e-7
 2.505210838544172e-8
 2.08767569878681e-9
 1.6059043836821616e-10
 1.1470745597729725e-11
 7.647163731819817e-13
 4.7794773323873865e-14
 2.8114572543455214e-15
 1.561920696858623e-16
 8.220635246624331e-18
 4.1103176233121653e-19
 1.9572941063391266e-20
 8.896791392450576e-22
 3.868170170630684e-23
 1.6117375710961186e-24
 6.446950284384475e-26
 2.4795962632247983e-27
 9.183689863795547e-29
 3.2798892370698385e-30
 1.1309962886447718e-31
 3.769987628815906e-33
eringh commented 3 years ago

However, it seems to be the evaluation of equation (52) in the paper that is missing:

1 ./ factorial.(Int128.(10:30))
21-element Vector{Float64}:
 2.755731922398589e-7
 2.505210838544172e-8
 2.08767569878681e-9
 1.6059043836821613e-10
 1.1470745597729725e-11
 7.647163731819816e-13
 4.779477332387385e-14
 2.8114572543455206e-15
 1.5619206968586225e-16
 8.22063524662433e-18
 4.110317623312165e-19
 1.9572941063391263e-20
 8.896791392450574e-22
 3.8681701706306835e-23
 1.6117375710961184e-24
 6.446950284384474e-26
 2.4795962632247972e-27
 9.183689863795546e-29
 3.279889237069838e-30
 1.1309962886447718e-31
 3.7699876288159054e-33
jarlebring commented 3 years ago

It might be that that the data in table 7 was not correctly used. It should be combined with (52) which it is not.

I thought there were also other procedures in that paper with numbers which were not tabulated.

eringh commented 3 years ago

I thought there were also other procedures in that paper with numbers which were not tabulated.

There seems to be a lot of versions where the coefficients are not tabulated. And a lot of versions that are based on slightly different evaluation schemes. Currently the function gen_sastre_basic_exp(k) is differentiating methods depending on the number of multiplications. But that is not uniquely determining which of the methods described in the paper that is used (although, at least for now, it is uniquely determining method combined with coefficients described in the paper).

How much is it wort to have all aspects of the paper implemented? It can get quite time consuming.

jarlebring commented 3 years ago

Currently the function gen_sastre_basic_exp(k) is differentiating methods depending on the number of multiplications

Then add a kwarg?

How much is it wort to have all aspects of the paper implemented? It can get quite time consuming.

It's better if you decide this yourself since you have the best overview right now. Start somewhere and see how far you want to go.

jarlebring commented 3 years ago

Just to be clear here: I think it's okay if the set of graphs is only a subset of everything in that paper considering that many values are not tabulated. We already have some graphs which are not tabulated (yes, I added the coefficients from an MEP-computation).

eringh commented 3 years ago

Sounds good. I was actually thinking of calling it done with possibly one exception: Do we want a conversion of the general framework (62)-(65) into Degopt? It might look good in the remark about related notation if we have a automatic conversion from y_ks to degopt. bild

jarlebring commented 3 years ago

Indeed, if we implement (62)-(65) it would be good to have it on degopt form.

eringh commented 3 years ago

Now the SID-paper can also be evaluated with the function gen_sastre_yks_degopt. Test for higher values of k gives further confidence in the implementation.