Pfaffian

[stevengj]

People are sometimes interested in computing this, and it's an interesting quantity defined for skew-Hermitian matrices. As for determinants, there are probably a large number of reasonable ways to compute it, but this algorithm is a possible starting point. (You could also use your Hessenberg factorization.)

As for the determinant, the Pfaffian will quickly overflow for large matrices, so you might want to compute its log instead.

[stevengj]

In particular, this formula from that paper makes it trivial to implement both the Pfaffian and the log-Pfaffian from your Hessenberg reduction, so it's worth doing that even if there are slightly more efficient ways to compute it: (noting that the determinant of your Q matrix in the Hessenberg factorization is just ±1 depending on whether the number reflectors is even/odd)

[stevengj]

This algorithm looks nice for exact Pfaffians of BigInt matrices: https://core.ac.uk/download/pdf/82502568.pdf

[stevengj]

In particular, here is an implementation of the exact Pfaffian of integer matrices using the algorithm from the paper above:

using SkewLinearAlgebra, LinearAlgebra
using LinearAlgebra: exactdiv

# in-place O(n³) algorithm to compute the exact Pfaffian of
# a skew-symmetric matrix over integers (or potentially any ring supporting exact division).
#
#     G. Galbiati & F. Maffioli, "On the computation of pfaffians,"
#     Discrete Appl. Math. 51, 269–275 (1994).
#     https://doi.org/10.1016/0166-218X(92)00034-J
function exactpfaffian!(A::AbstractMatrix{<:Integer})
    LinearAlgebra.require_one_based_indexing(A)
    isskewhermitian(A) || throw(ArgumentError("Pfaffian requires a skew-Hermitian matrix"))
    n = size(A,1)
    isodd(n) && return zero(eltype(A))
    c = one(eltype(A))
    signflip = false
    n = n ÷ 2
    while n > 1
        # find last k with A[2n-1,k] ≠ 0
        k = 2n
        while k > 0 && iszero(A[2n-1,k]); k -= 1; end
        iszero(k) && return zero(eltype(A))

        # swap rows/cols k and 2n
        if k != 2n
            for i = 1:2n
                A[i,k], A[i,2n] = A[i,2n], A[i,k]
                A[k,i], A[2n,i] = A[2n,i], A[k,i]
            end
            signflip = !signflip
        end

        # update, A, c, n
        for j = 1:2n-2, i = 1:j-1
            δ = A[2n-1,2n]*A[i,j] - A[i,2n-1]*A[j,2n] + A[j,2n-1]*A[i,2n]
            A[j,i] = -(A[i,j] = exactdiv(δ, c))
            # @assert A[i,j] * c == δ
        end
        c = A[2n-1,2n]
        n -= 1
    end
    return signflip ? -A[1,2] : A[1,2]
end
exactpfaffian(A::AbstractMatrix{<:Integer}) =
    exactpfaffian!(copyto!(similar(A, typeof(exactdiv(zero(eltype(A))^2, one(eltype(A))))), A))

You probably only want to use this for AbstractMatrix{<:BigInt} by default, since other types are likely to overflow etcetera, similar to how the LinearAlgebra package only uses the Bareiss algorithm by default for determinants of BigInt matrices (see https://github.com/JuliaLang/julia/pull/40868).

[stevengj]

And here's another paper on numerical algorithms: https://arxiv.org/abs/1012.5022

One of the algorithms they discuss is "tridiagonalization based on Householder transformations", i.e. by Hessenberg factorization as I suggested above. Probably I would just go with that approach since it is easiest.

Note also that they suggest that the Pfaffian of complex skew-Hermitian matrices is useful for some electronic-structure calculations, which adds some additional motivation to #15.

[smataigne]

Thank you very much for the documentation and the code above. I will add the functionality to the package.

[stevengj]

Closed by #47