pfaffian(A::SkewHermitian{<:Complex}) is missing

[jlbosse]

Currently there is no support for computing pfaffians of complex matrices. TopologicalNumbers.jl has an implementation of the algorithms from Efficient numerical computation of the pfaffian for dense and banded skew-symmetric matrices but is a fairly heavy package.

If there is interest I can open a PR adding support for pfaffians of complex matrices.

[stevengj]

The Pfaffian is defined for complex skew-antisymmetric matrices $A = -A^T$, but there is no conventional definition (as far as I can tell) for complex skew-Hermitian matrices $A = -A^*$.

This package doesn't have any type for complex skew-symmetric matrices $A = A^T$ right now. I suppose we could add one, but it's a fairly heavyweight change just to add one operation.

(We had a discussion offline about this with @alanedelman in 2022, and the conclusion was that there does not seem to be any useful natural definition of a Pfaffian for complex skew-Hermitian matrices.)

[stevengj]

Actually, I suspect that the SkewLinearAlgebra._exactpfaffian! function probably works as-is for complex skew-symmetric matrices (though I haven't tested it). So, if we didn't want to add a new type, all that would be needed would be to add a method like:

pfaffian(A::AbstractMatrix{T}) where {T<:Complex} =
    pfaffian!(copyto!(similar(A, typeof(zero(T)/one(T))), A))
function pfaffian!(A::AbstractMatrix{<:Complex})
    LinearAlgebra.require_one_based_indexing(A)
    isskewsymmetric(A) || throw(ArgumentError("Pfaffian requires a skew-symmetric matrix"))
    return _exactpfaffian!(A)
end

along with a function isskewsymmetric(A).

This seems like it would be a useful addition if you want to add it along with tests and documentation.

[jlbosse]

Ah sorry, yes I did indeed mean complex skew-symmetric, not skew-hermitian. Regarding using _exactpfaffian! instead of e.g. one of the algorithms described in Efficient numerical computation of the pfaffian for dense and banded skew-symmetric matrices: Is exact_pfaffian numerically stable enough and/or fast enough?

[stevengj]

As for speed, all algorithms for the Pfaffian of a dense $n \times n$ matrix are presumably $\Omega(n^3)$, and _exactpfaffian! is $\Theta(n^3)$ so it should be fine up to constant factors.

No idea about numerical stability for inexact arithmetic, but experimentally it seems fine:

julia> A = skewhermitian!(randn(100,100));

julia> p1, p2 = pfaffian(A), SkewLinearAlgebra._exactpfaffian!(copy(A.data))
(-7.982243401303055e30, -7.982243401308707e30)

julia> abs(p1 - p2) / abs(p1)
7.080738143649332e-13

julia> A = skewhermitian!(randn(300,300));

julia> p1, p2 = pfaffian(A), SkewLinearAlgebra._exactpfaffian!(copy(A.data))
(4.151962041555626e129, 4.151962041553932e129)

julia> abs(p1 - p2) / abs(p1)
4.0792054290158497e-13

though it does seem more susceptible to overflow than I would like:

julia> A = skewhermitian!(randn(500,500));

julia> p1, p2 = pfaffian(A), SkewLinearAlgebra._exactpfaffian!(copy(A.data))
(-1.5562908514995813e244, NaN)

but of course, once you get that big you will quickly overflow any algorithm and will have to switch to a logpfaffian method or similar.

[stevengj]

It would be great to benchmark it (for both performance and accuracy) against the algorithm in TopologicalNumbers.jl, though.

[smataigne]

Thank you @jlbosse for your remarks! And thank you @stevengj for your interesting answers!