[ANN] GraphMatFun.jl: computation graphs for matrix functions

Package Announcements
, ,
1

We are proud to announce the first release of the package GraphMatFun.jl. This package contains functionality to represent, manipulate and optimize algorithms for matrix functions using directed acyclic graphs. It also contains code-generation features for several languages (Julia, MATLAB, C/BLAS).

The Padé-based scaling and squaring algorithm is represented in this way:

Screenshot_20210810_100543
Screenshot_20210810_100543786×532 55.9 KB

The package has a corresponding manuscript and data files repository. In the data repository you can download optimized graphs and generated code to be used directly in computations. For example, the matrix exponential can in many cases be computed more efficiently than the native implementation:

julia> using GraphMatFun, BenchmarkTools, LinearAlgebra;
julia> # Download from https://github.com/matrixfunctions/GraphMatFunData/tree/main/data/generated/exp
julia> include("exp_sid_m5_opt_rho1_9.jl");
julia> A=randn(1000,1000); A=1.8*A/norm(A);
julia> E1=@btime exp($A);
  1.111 s (23 allocations: 76.31 MiB)
julia> E2=@btime exp_sid_m5_opt_rho1_9($A);
  332.667 ms (19 allocations: 68.67 MiB)
julia> norm(E1-E2)/norm(E2)
3.1592474046054336e-16
9 Likes
2

More about the example: The example illustrates something I believe we can expect in the future. The matrix exponential can be more efficiently computed with efficient polynomial evaluation schemes, instead of what has been state-of-the-art for very long: Padé based methods. Polynomial based methods only require matrix-matrix multiplications in contrast to Padé-based which require an inverse (or rather a left matrix division) which is more complicated to implement and parallelize than matrix-matrix multiplication. In MKL, left-division is more efficiently implemented than in OpenBLAS, which is why the difference in CPU-time in the example above will be smaller if you add using MKL. On the other hand, high order polynomial evaluations are formed from linear combinations of many matrices, which is perfect for Julia’s dot-fusion functionality. Differences are illustrated in the manuscript.