A Julia package for computing scalar functions of matrix variables and their Fréchet derivatives. The matrix functions computation (for arbitrary square matrices) is based on the Schur-Parlett algorithm (with improvements)a. The higher order Fréchet derivatives (for Hermitian matrices) are formulated similarly to the Daleckii-Krein theorem, where the divided differences are calculated accurately by the Opitz' formula. In particular, MatrixFuns supports the computation of discontinuous functions.
For smooth functions, such as exp:
julia> using MatrixFuns
julia> using LinearAlgebra
julia> A = [1 1 0; 0 1.1 1; 0 0 1.2]
3×3 Matrix{Float64}:
1.0 1.0 0.0
0.0 1.1 1.0
0.0 0.0 1.2
julia> mat_fun(exp, A) # computed by Schur-Parlett
3×3 Matrix{Float64}:
2.71828 2.85884 1.50334
0.0 3.00417 3.15951
0.0 0.0 3.32012
julia> X = 0.05 * Symmetric(A) # generate a symmetric matrix
3×3 Symmetric{Float64, Matrix{Float64}}:
0.05 0.05 0.0
0.05 0.055 0.05
0.0 0.05 0.06
julia> a = eigvals(X)
3-element Vector{Float64}:
-0.01588723439378913
0.055000000000000014
0.12588723439378913
julia> div_diff(exp, a) # returns the 2nd order divided difference
0.5284915575854794
julia> hs = map(x->x*X, [1, 2])
2-element Vector{Symmetric{Float64, Matrix{Float64}}}:
[0.05 0.05 0.0; 0.05 0.05500000000000001 0.05; 0.0 0.05 0.06]
[0.1 0.1 0.0; 0.1 0.11000000000000001 0.1; 0.0 0.1 0.12]
julia> mat_fun_frechet(exp, X, hs) # returns the 2nd order Fréchet derivative d^2exp(X)hs_1hs_2
3×3 Matrix{Float64}:
0.0110862 0.0119138 0.00588556
0.0119138 0.0181632 0.0130909
0.00588556 0.0130909 0.0135867For special functions, such as the error function erf:
julia> using SpecialFunctions
julia> mat_fun(erf, A)
3×3 Matrix{Float64}:
0.842701 0.375043 -0.369768
0.0 0.880205 0.301089
0.0 0.0 0.910314For singular functions, such as Fermi-Dirac functions with temperatures close to 0, user can set a smaller scale to reduce the spread of each block to avoid large Taylor expansion errors near the singularities, and can also customize color function to avoid the block's spectral range including singularities.
julia> μ = 1.15
julia> f(x) = 1/(1+exp(1000*(x-μ))); # Fermi-Dirac function with temperature equal to 1e-3
julia> color(x) = x < μ ? 1 : 2; # use `color` to avoid singularities.
julia> mat_fun(f, A; scale=0.01, color)
3×3 Matrix{Float64}:
1.0 0.0 -50.0
0.0 1.0 -10.0
0.0 0.0 1.92875e-22For piecewise functions consisting of continuous intervals, user can customize color function to first split the eigenvalues into different continuous intervals.
julia> x0 = 0.051 # discontinuous point
julia> g(x) = x < x0 ? sin(x) : (x+1)^2;
julia> color(x) = x < x0 ? 1 : 2; # specify a different color for each continuous interval
julia> mat_fun_frechet(g, X, hs; color)
3×3 Matrix{Float64}:
0.019713 0.0213782 0.00975085
0.0213782 0.0316017 0.0233283
0.00975085 0.0233283 0.0241837Specially, for constant functions, or piecewise functions consisting of constant intervals, user can further set sep=Inf to only split the eigenvalues by color, which will be more efficient.
julia> color(x) = Int(sign(x))
julia> div_diff(sign, -1.1, 0.02, -0.01, 0.3; color, sep=Inf) # returns the 3rd order divided difference sign[-1.1, 0.02, -0.01, 0.3]
-196.12683783190693For more details, please see the documentation.
julia> using Pkg
julia> Pkg.add("MatrixFuns")