Built-in graphs

 (graph,crefs)=graph_monomial(a; input=:A, polyname=:P)

Generates the graph for the polynomial using the monomial basis coefficients. More precisely, it corresponds to the evaluation of the polynomial

p(s) = a[1] + a[2]*s + ... + a[n]*s^(n-1),

where s^k is naively evaluated as s^k=s*s^(k-1) for k=2,3,...,n-1. The kwarg polyname specifies the name of intermediate variables.

 (graph,crefs)=graph_monomial_degopt(a; input=:A)

Generates the same polynomial as graph_monomial, in the monomial basis. However, it does so by wrapping a call to graph_degopt, resulting in more degrees of freedom in crefs.

 (graph,crefs)=graph_horner(a; input=:A, B=:B, C=:C, scaling=1.0)

Generates the graph for the polynomial using Horner's method. More precisely, it corresponds to the evaluation of the polynomial

p(s) = a[1] + a[2]*(αs) + ... + a[n-1]*(αs)^(n-2) + a[n]*(αs)^(n-1),

as the recursion

Cj=s*B{j+1}
Bj=a[j]*I+α*Cj,

where α=scaling.

The kwargs B and C specifies the base-names of these intermediate variables.

 (graph,crefs)=graph_horner_degopt(a; scaling=1.0, input=:A)

Generates a polynomial using Horner's evaluation scheme. The polynomial

 p(s) = a[1] + a[2]*(αs) + ... + a[n-1]*(αs)^(n-2) + a[n]*(αs)^(n-1),

is evaluated as

 p(s) = a[1] + (αs)*(a[2] + (αs)*(... + (αs)*(a[n-1] + a[n]*(αs))...)),

where α=scaling. However, the function uses a call to graph_degopt, resulting in more degrees of freedom in crefs. See also graph_horner.

 (graph,crefs)=graph_ps(a; input=:A,
                      B_base=:B, C_base=:C, P_base=:P)

Generates the graph for the Paterson–Stockmeyer procedure with monomial basis coefficieents. More precisely, it corresponds to evaluation of the polynomial

p(s) = a[1] + a[2]*s + ... + a[n]*s^(n-1).

The code follows the description in the second of the papers referenced below.

References

1. M. Paterson and L. Stockmeyer. "On the number of nonscalar multiplications necessary to evaluate polynomials". SIAM Journal on Scientific Computing, 2(1):60-66, 1973. DOI: 10.1137/0202007

2. M. Fasi. "Optimality of the Paterson–Stockmeyer method for evaluating matrix polynomials and rational matrix functions". Linear Algebra and its Applications, 574, 2019. DOI: 10.1016/j.laa.2019.04.001

 (graph,crefs)=graph_ps_degopt(a; input=:A)

Generates the same polynomial as graph_ps, i.e., the graph for the Paterson–Stockmeyer procedure with monomial basis coefficieents. However, it does so by wrapping a call to graph_degopt, resulting in more degrees of freedom in crefs.

(graph,cref)=graph_sastre_poly(b)

Computes the degree-8 polynomial

p(z)=b[1]+z*b[2]+z^2*b[3]+...+z^8*b[9]

according to Example 3.1 in the reference.

Reference

1. J. Sastre. "Efficient evaluation of matrix polynomials". Linear Algebra and its Applications, 539:229-250, 2018. DOI: 10.1016/j.laa.2017.11.010

degopt=Degopt(x,y)
degopt=Degopt(HA,HB,y)
degopt=Degopt(graph)

struct Degopt{T}
    x::Vector{Tuple{Vector{T},Vector{T}}}
    y::Vector{T}
end

Creates an object representing a degree-optimal polynomial, i.e., the coefficients in an evaluation scheme which maximizes the degree for a fixed number of multiplications. The object represents the coefficients in

A0=I
A1=A
A2=(x *A0+x *A1)(x *A0+x *A1)
A4=(x *A0+x *A1+x*A2)(x *A0+x *A1+x *A2)
A8=(x *A0+x *A1+x*A2)(x *A0+x *A1+x *A2)
 ..

and

Out=y*A0+y*A1+y*A4+y*A8+y*B4...

The x-values are given in the argument x, which is a Vector{Tuple{Vector{T},Vector{T}}}, containing the elements of each sum. The y-vector contains the elements to form the output. In the constructor with HA, HB and y the elements of x are stored as matrices.

The coefficients in graphs generated by graph_degopt can be recovered by using Degopt(graph).

(graph,crefs)=graph_degopt(k;T=ComplexF64,input=:A)
(graph,crefs)=graph_degopt(x,z;input=:A)
(graph,crefs)=graph_degopt(d::Degopt;input=:A)

Corresponds to the (for a fixed number of multiplications) degree-optimal polynomial

B1=A
B2=(x *I+x *A)(x *I+x *A)
B3=(x *I+x *A+x*B2)(x *I+x *A+x *B2)
 ..

and

Out=z*I+z*A1+z*B2+z*B3+z*B4...

The x-values are given in the argument x, which is a Vector{Tuple{Vector,Vector}}, containing the elements of each sum. The z-vector contains the elements to form the output, and input determines the name of the matrix A above. If the parameter k is supplied instead of the coefficients, all coeffs will be set to one. The general recursion is given in (9) in the paper referenced below.

Reference

1. P. Bader, S. Blanes, and F. Casas, "Computing the matrix exponential with an optimized Taylor polynomial approximation", Mathematics, 7(12), 2019. DOI: 10.3390/math7121174

grow!(degopt::Degopt)

Increases the Degopt by one multiplication without modifying the function values.

scale!(degopt::Degopt,α)

Effectively change a Degopt such that the input is scaled by α. If p is the original function, p(α x) will be the new function.

square!(degopt::Degopt)

Effectively square a Degopt in the sense that the output is square. If p is the original function, p(x)^2 will be the new function.

normalize!(degopt::Degopt,tp=:row1) -> degopt

Normalizes the Degopt coefficients, in the way specified by tp. If the tp==:row1 the degopt will be transformed to an equivalent Degopt with first row equal to (0 1) (0 1). If tp==:col1 the first column in the Ha and Hb matrices will be transformed to zero.

get_degopt_coeffs(degopt) -> (HA,HB,y)
get_degopt_coeffs(graph) -> (HA,HB,y)

Returns the coefficients of the degree-optimal polynomial using two matrices and one vectors as described in Degopt.

(x,z)=get_degopt_crefs(k)
(x,z)=get_degopt_crefs(graph)

Returns linear combination references (crefs) related to graph_degopt. Specifically x is a Vector{Tuple{Vector{Tuple{Symbol,Int}},Vector{Tuple{Symbol,Int}}}} such that x[2][1] corresponds to the coefficients of the left hand side of the multiplication

B2=(α_2_1 *I + α_2_2 *A)(β_2_1 *I + β_2_2 *A)

i.e., the crefs corresponding to [α_2_1, α_2_2]. See graph_degopt. Hence, get_coeffs(graph,x[2][1]) returns the corresponding numerical values of the coefficients.

(graph,cref)=graph_sastre_yks_degopt(k,s,c)

Transforms the polynomial evaluation format given by equations (62)-(65) in the reference to degop-format. The graph is a representation of y_{k,s}. Input c is a grouping of the coefficients as given by the representation (62)-(65). c is a Vector{Vector{Vector}} of length k+1, representing

[
[c_i^{(0,1)}, c_i^{(0,2)}]
[c_i^{(1,1)}, c_i^{(1,2)}, c_i^{(1,3)}, c_i^{(1,4)}, c_i^{(1,5)}, c_i^{(1,6)}]
...
[c_i^{(k,1)}, c_i^{(k,2)}, c_i^{(k,3)}, c_i^{(k,4)}, c_i^{(k,5)}, c_i^{(k,6)}]
]

Hence, c[1] contains two vectors, the first of length s-1 and the second of length s. (Note: In the first vector the constant for I is set to zero) and c[2] up to c[k+1] conatins six vectors: Even-numbered vectors have length s+1 Odd-numbered vectors have length j-1, where j is the intex in c, e.g., c[2][1] has one element and c[3][5] have two. For example, equations (57)-(59) are implemented as:

c = [
[[c15, c16], [0.0, 0, 0]], # (57)
[[1.0], [0, c13, c14], [1.0], [c11, 0, c12], [c10], [0.0, 0, 0]], # (58)
[[0.0, 1], [0, c8, c9], [c7, 1], [0, c6, 0], [c4, c5], [c1, c2, c3]], # (59)
]

Reference

1. J. Sastre. "Efficient evaluation of matrix polynomials". Linear Algebra and its Applications, 539:229-250, 2018. DOI: 10.1016/j.laa.2017.11.010

 (graph, cref) = graph_rational(den_coeffs, num_coeffs, poly_gen=graph_ps)
 (graph, cref) = graph_rational(
    den_graph,
    num_graph;
    den_cref = Vector{Tuple{Symbol,Int}}(),
    num_cref = Vector{Tuple{Symbol,Int}}(),
)

Generates the graph for the rational approximation

r(A)=q(A)^{-1}p(A)

where p(A) and q(A) are polynomials defined by the coefficients den_coeffs and num_coeffs, and generated by the function poly_gen, which is called as (graph,cref)=poly_gen(coeffs), see, e.g., graph_monomial and graph_ps.

The alternative call-signature involves the graphs for p and q directly, as den_graph and num_graph. The corresponding den_cref and num_cref can also be passed to be modified accordingly, otherwise the return value cref is empty for this call.

 (graph,crefs)=graph_exp_native_jl(A; input=:A)

Creates a graph for the native scaling-and-squaring for the matrix exponential, as implemented in Julia. The matrix A is taken as input to determine the length of the Padé approximant and the number of squares applied, as well as determining the type of the coefficients. The kwarg input determines the name of the matrix, in the graph.

References

1. N. J. Higham. "Functions of Matrices". SIAM, Philadelphia, PA, 2008. DOI: 10.1137/1.9780898717778

2. N. J. Higham. "The Scaling and Squaring Method for the Matrix Exponential Revisited". SIAM Journal on Matrix Analysis and Applications, 26(4):1179-1193, 2005. DOI: 10.1137/04061101X

3. "Julia's matrix exponential", at the time of conversion.

 (graph,crefs)=graph_exp_native_jl_degopt(A; input=:A)

Same as graph_exp_native_jl but with calls to graph_degopt for contruction of numerator and denominator polynomials.

(graph,cref)=graph_sastre_exp(k,method=:auto)

Computes a polynomial evaluation approximating the exponential using k matrix multiplications following a method given in the reference. The schemes are embedded into degop-format, see graph_degopt.

The methods are (from the paper referenced below):

• :ps_degopt, Paterson–Stockmeyer method embedded into degopt-format.
• :y1s, given by equations (34)-(35)
• :z1ps, given by equations (34)-(35) and (52)
• :h2m, given by equations (34)-(35) and (69)

Not all combinations of k and method are implemented. Available ones are:

• k<3, method=:ps_degopt
• k=3, method=:y1s, as per Table 4 in the reference
• k=4, method=:y1s
• k=6, method=:h2m, as per Table 11 in the reference
• k=8, method=:z1ps, as per Table 7 in the reference

The default option method=:auto will choose method according to the value of k, as prescribed above.

Reference

1. J. Sastre. "Efficient evaluation of matrix polynomials". Linear Algebra and its Applications, 539:229-250, 2018. DOI: 10.1016/j.laa.2017.11.010

(graph,cref)=graph_sid_exp(k;T=Float64)

Computes a polynomial evaluation approximating the exponential using k matrix multiplications following the procedure in the reference. The coefficients are directly copied from the paper.

The evaluation is embedded in the degopt-format, see graph_degopt, using the function graph_sastre_yks_degopt. Moreover, for k<=3 it uses the Paterson–Stockmeyer method.

Reference

1. J. Sastre, J. Ibáñez, E. Defez. "Boosting the computation of the matrix exponential". Applied Mathematics of Computation, 340:206-220, 2019. DOI: 10.1016/j.amc.2018.08.017

(graph,cref)=graph_bbc_exp(k;T=Float64)

Computes a polynomial evaluation approximating the exponential using k matrix multiplications following the procedure in the reference. The coefficients are directly copied from the paper.

The evaluation is embedded in the degopt-format, see graph_degopt, and for k< 3, the evaluation is using the Paterson–Stockmeyer method.

Reference

1. P. Bader, S. Blanes, and F. Casas. "Computing the matrix exponential with an optimized Taylor polynomial approximation". Mathematics, 7(12), 2019. DOI: 10.3390/math7121174

(graph,cref)=graph_bbcs_cheb_exp(k;T=Complex{BigFloat})

Computes a polynomial evaluation approximating the exponential for skew-Hermitian matrices using k matrix multiplications following the procedure in the reference.

For k > 5 it resorts to scaling-and-squaring of the k = 5 graph. The graph is in the degopt format, see graph_degopt.

Reference

1. P. Bader, S. Blanes, F. Casas, M. Seydaoglu. "An efficient algorithm to compute the exponential of skew-Hermitian matrices for the time integration of the Schrödinger equation". arXiv:2103.10132 [math.NA], 2021.

    (graph,cref)=graph_denman_beavers(k;T=ComplexF64,cref_mode=0)

Creates the graph corresponding to the Denman–Beavers iteration for the matrix square root, using k iterations. The kwarg cref_mode specifies if references to all X and Y (cref_mode=0) should be stored or just Y (cref_mode=-1) or just X (cref_mode=1).

Reference

1. E. Denman and A. Beavers. "The matrix sign function and computations in systems". Applied Mathematics and Computation, 2(1), 1976. DOI: 10.1016/0096-3003(76)90020-5

 (graph,crefs)=graph_newton_schulz(
    k,
    T = ComplexF64;
    input = :A,
    Z = :Z,
    Q = :Q,
    V = :V,
)

Does k iterations of the Newton–Schulz iteration for approximating the inverse of a matrix A (name given by input), i.e., the recursion

V_k+1 = V_k*(2*I - A*V_k),

with V_0=A. The recursion is implemented using the graph operations

Z_i=A*V_i
Q_i=2*I-Z_i
V_{i+1}=V_i*Q_i,

and the kwargs Z, Q, and V determines the naming of the corresponding intermediate steps.

References

1. G. Schulz. "Iterative Berechnung der reziproken Matrix". Zeitschrift für Angewandte Mathematik und Mechanik, 13:57–59, 1933. DOI: 10.1002/zamm.19330130111

2. N. J. Higham. "Functions of Matrices". SIAM, Philadelphia, PA, 2008. DOI: 10.1137/1.9780898717778

 (graph,crefs)=graph_newton_schulz_degopt(k, T=ComplexF64; input=:A)

Does k iterations of the Newton–Schulz iteration for approximating the inverse, using the recursion

Z_i=A*V_i
V_{i+1}=V_i*(2*I-Z_i).

The function makes a call to graph_degopt, resulting in more degrees of freedom in crefs. See also graph_newton_schulz.