pFq with Matrix argument

[MikaelSlevinsky]

This could be pretty useful. In Base, exp(z::Matrix{BigFloat}) doesn't exist but we could support it with some minor modifications of the rational algorithms

julia> using LinearAlgebra

julia> @inline errcheck(x::Matrix, y::Matrix, tol) = (norm(x-y) > max(norm(x), norm(y))*tol)
errcheck (generic function with 1 method)

julia> function pFqdrummond(::Tuple{}, ::Tuple{}, z::Matrix{T}; kmax::Int = 10_000) where T
           if norm(z) < eps(real(T))
               return Matrix{T}(I, size(z, 1), size(z, 2))
           end
           ζ = inv(z)
           Nlo = Matrix{T}(I, size(z, 1), size(z, 2))
           Dlo = Matrix{T}(I, size(z, 1), size(z, 2))
           Tlo = Nlo/Dlo
           Nhi = (2I - z)*Nlo + 2z
           Dhi = (2I - z)*Dlo
           Thi = Nhi/Dhi
           k = 1
           Nhi, Nlo = ((k+2)*I-z)*Nhi + k*z*Nlo + z*z, Nhi
           Dhi, Dlo = ((k+2)*I-z)*Dhi + k*z*Dlo, Dhi
           Thi, Tlo = Nhi/Dhi, Thi
           k += 1
           while k < kmax && errcheck(Tlo, Thi, 10eps(real(T)))
               Nhi, Nlo = ((k+2)*I-z)*Nhi + k*z*Nlo, Nhi
               Dhi, Dlo = ((k+2)*I-z)*Dhi + k*z*Dlo, Dhi
               Thi, Tlo = Nhi/Dhi, Thi
               k += 1
           end
           return Thi
       end
pFqdrummond (generic function with 1 method)

julia> Z = big.(rand(2, 2))
2×2 Matrix{BigFloat}:
 0.82811   0.0613331
 0.449884  0.073366

julia> pFqdrummond((), (), Z)
2×2 Matrix{BigFloat}:
 2.31398   0.0990111
 0.726256  1.09558

julia> exp(Float64.(Z))
2×2 Matrix{Float64}:
 2.31398   0.0990111
 0.726256  1.09558

julia> pFqdrummond((), (), Z) - exp(Float64.(Z))
2×2 Matrix{BigFloat}:
 6.85325e-16   3.73548e-17
 2.05869e-16  -7.37797e-17

julia> exp(Z)
ERROR: MethodError: no method matching exp(::Matrix{BigFloat})
Closest candidates are:
  exp(::Union{Float16, Float32, Float64}) at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/base/special/exp.jl:296
  exp(::StridedMatrix{var"#s859"} where var"#s859"<:Union{Float32, Float64, ComplexF32, ComplexF64}) at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/dense.jl:560
  exp(::StridedMatrix{var"#s859"} where var"#s859"<:Union{Integer, Complex{<:Integer}}) at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/dense.jl:561
  ...
Stacktrace:
 [1] top-level scope
   @ REPL[7]:1

julia> 

[dlfivefifty]

FYI Alan had some RMT work related to matrix hypergeom functions: https://math.mit.edu/~edelman/publications/computing_with_beta.pdf

[MikaelSlevinsky]

I'm confused by their equation (3.1). Given a diagonal matrix X = diag(x_1,...x_m), isn't exp(tr(X)) a scalar?

[dlfivefifty]

I think they are all scalar-valued

[MikaelSlevinsky]

Hmm, are they the determinants of the matrix-valued pFq by the standard definition? It seems to hold for ₀F₀^(β) = exp(tr(X)), and ₁F₀^(β) = |I-X|^{-a}.

If there's no obvious relationship, I'm not sure they're in the scope of this package.

[dlfivefifty]

I doubt it since one of them is a pfaffian

[aravindh-krishnamoorthy]

Hello @MikaelSlevinsky, in case you haven't considered these already, I'd be very much interested in understanding the trade-offs (computational complexity, accuracy, and domain) of the methods in New Solvable Matrix Integrals by AY Orlov; for a simple example, please see my blog post here.

A better (?) alternative is perhaps to implement versions based on Zonal polynomials. You're probably aware of these reference publications and reference software.

In case you're interested in adding matrix-valued functions with matrix arguments to this package, I could also help with evaluation and implementation a bit later (after one more implementation on my TODO list). I'd especially be interested in a comparison between Orlov's, Koev's, and other new methods.

[MikaelSlevinsky]

I guess I still haven't wrapped my head around why the hypergeometric functions of matrix argument (in the literature) do not agree with the definition through their Maclaurin series.

[aravindh-krishnamoorthy]

@MikaelSlevinsky I can tell you my own story... Perhaps guys at MathOverflow or Reddit Math can help with a broader historical context.

I first came across a good application of scalar-valued Gamma, Bessel, and hypergeometric functions of matrix arguments in [1] (please see Note 1) in the context of matrix-variate distributions. Turns out that several interesting integrals over matrices can be expressed in terms of these functions and evaluated efficiently using zonal polynomials, see Koev's links in my previous post. Alternatively, the analytical framework in Orlov's papers is useful if we want to work further on the results. I also found [2], which might give you a bit more context on the utility of these functions in "hot" topics such as ML :)

I hope this gives you a data point on why the scalar-valued definition is useful. As alluded to in [Preface, 3], there can be other definitions of hypergeometric functions of matrix arguments as well, which agree with the series definitions of the scalar equivalents.

Note 1: Please note that you can get a free "preview" pdf in the link given below. Unfortunately, Chapter 1.6 on hypergeometric functions is missing, but Chapter 1.5 on zonal polynomials is present.

[1]	Gupta, A.K. and Nagar, D.K., 2018. Matrix variate distributions. Chapman and Hall/CRC. Preview: https://www.taylorfrancis.com/books/mono/10.1201/9780203749289/matrix-variate-distributions-nagar-gupta
[2]	A. Karbalayghareh, X. Qian and E. R. Dougherty, "Optimal Bayesian Transfer Learning," in IEEE Transactions on Signal Processing, vol. 66, no. 14, pp. 3724-3739, 15 July15, 2018, doi: 10.1109/TSP.2018.2839583.
[3]	Higham, N.J., 2008. Functions of matrices: theory and computation. Society for Industrial and Applied Mathematics.