Move graph generators

[jarlebring]

When we transfer graph generators to github, I suggest we use the opportunity to improve the code.

• All generators should return (graph,cref)
• We can this function name convention:
  • Function gen_X(): generates a graph corresponding to technique X, or
  • Function gen_X_F(): generates a graph corresponding to technique X for matrix function F
• The type of the coeffs can be specified by kwarg T=Float64. (see example in Denman-Beavers c433773 ) unless given implicitly by other parameters.
• All generator functions should have julia style comment, including a reference. See example in denman beavers c433773.
• Every generator should have at least one unit test.

To-do list for generators:

• [x] gen_denman_beavers (from dropbox denman_beavers.jl. remove cref_mode-kwarg)
• [x] monomial.jl: gen_monomial from monomial_graph_gen(). Replace explicit summing with call to add_sum! (Maybe there is a better name than monomial, e.g. monomial_poly)
• [x] pade.jl: gen_pade. Replace explicit summing to a call to one of the polynomial functions (user can specify gen_monomial, gen_ps, etc)
• [x] ps_poly.jl: gen_ps
• [x] newton_schulz_graph_gen.jl: gen_newton_schulz. Should return (graph,cref).
• [x] horner_graph_gen.jl -> gen_horner. Make kwargs consistent with gen_monomial. Return two (graph,cref)
• [x] expm: julia native. call add_sum!.
• expm2009: call add_sum!.
• sastre exponential: The exponential of sastres paper, with hard-coded coeffs. Call gen_general_poly_recursion.
• sastre general. The general function for sufficiently small p. Call gen_general_poly_recursion.
• [x] bbc general. Clean up function names and arguments. Now: gen_general_poly_recursion.
• [x] bbc exponential. Hard-coded coefficients as in the paper. Call gen_general_poly_recursion.

[jarlebring]

Concerning the expm versions: The comparison in the BBC-paper: https://doi.org/10.3390/math7121174 contains these expm-versions:

The notation expm2005, expm2009, expm2016, seems natural.

Question (maybe for @mfasi): Are there better/other versions now? Which version is in Julia?

As far as I understand, tweaks like in expm2016 are for problems which are not the main the main target class for us. I think we can probably do these tweaks but nothing is gained from the graph viewpoint.

[mfasi]

Are there better/other versions now?

There is a multiprecision version published in 2019, but even though we run some experiments in double precision there I doubt that algorithm would be of interest in this context. A more relevant recent implementation for double precision is in this paper, which is already in the list above.

Which version is in Julia?

Julia implements the version in the book Functions of Matrices: Theory and Computation, which is essentially expm2005 in this paper. For reference, the code is here.

[eringh]

I implemented a function for generating rational functions, rather than calling it Padé, but ticked to pade.jl-box anyway.

[eringh]

How do we want scaling to work for the polynomial generators? Do we want the behavior to align across different methods?

Currently the scaling works a bit different among them. For ps_poly and monomial the scaling is purely applied to the coefficients (in cref). Hence, if new coefficients are set (using the returned cref), then the scaling disappear. However, in horner the scaling is mostly (but not completely!) applied to the other coefficients (not part of cref). Hence, if new coefficients are set (using cref), then it is difficult to overlook what polynomial is evaluated.

[jarlebring]

How do we want scaling to work for the polynomial generators? Do we want the behavior to align across different methods?

I would say no. Different methods have different ways of naturally doing a scaling. Not all methods need to have scaling. If scaling can be achieved by preprocessing, it seems too basic to have as a feature.

[eringh]

What do we want in the cref for Newton–Schultz? I currently left it empty.

[jarlebring]

What do we want in the cref for Newton–Schultz? I currently left it empty.

Let's decide when we reconnect with simulations based on that function.

[eringh]

I am currently looking into the expm2009 from this paper . It seems to need some 1-norm estimation. The reference code from BBC uses the matlab-function normest1. Does anyone know if there is such functionality in Julia? Or should we implement something like this? I can only find old references to norm-estimation in the Julia standard library. Perhaps @mfasi knows, you seemed to have been involved in some of the discussions.

[mfasi]

I can confirm that the functionality is not in Julia. Several people tried to get that algorithm into Base, but we all failed.

The most recent implementation, I think, is the one you can find in this PR, specifically here. This function estimates the 1-norm of a matrix given explicitly, whereas for expm we need to estimate the 1-norm of a matrix given only a function that computes the action of that matrix (and its conjugate transpose) on a vector, but it's close enough.

If you're OK with that, I can implement that function, as I've already worked with that algorithm and have an idea of how it should work. If you want to use the functionality in the mean time, you can probably just define a one-liner that computes the 1-norm of the pth power of A instead of estimating it:

# Estimate 1-norm of A^p.
opnormAp(A::Matrix,p::Integer)=opnorm(A^p,1)

[eringh]

If you're OK with that, I can implement that function, as I've already worked with that algorithm and have an idea of how it should work. If you want to use the functionality in the mean time, you can probably just define a one-liner that computes the 1-norm of the pth power of A instead of estimating it:

That would be great! I did a line-by-line conversion of the expm2009 code and simply added a placeholder here.

[eringh]

I had a talk with @jarlebring and we thought that all the glue-code for expm2009 is possibly not needed in the generators. Later we might try to add it to generated code, but for now we halt this. @mfasi no need to implement an 1-norm estimator.

[jarlebring]

Sorry for not following this thread carefully. I didn't realize the 1-norm estimator work was initiated. As I mentioned to Emil, I think the glue code for graph selection needs separate attention (e.g. we think through how it best fits with the generated code), and is maybe less central for this package which is focused on graphs rather than selecting the correct graph.

[mfasi]

I implemented the code for 1-norm estimation, but won't commit it to the repository as it is not clear whether we're gonna need it or not. For now, it's in this gist.

[jarlebring]

Thanks! It will probably be useful for glue-code later

[jarlebring]

Moved the sastre code gen to #30. Closes this issue.