Added pfaffian(::AbstractMatrix{<:Complex}) and isskewhermitian (#137)

So far I added only computation of pfaffians for complex skew symmetric
matrices, but if we want to take them seriously we should probably also
add a type `SkewSymmetric` and use that for pfaffian computations
instead of `SkewSymmetric` since that is what pfaffians are actually
defined for. However, this is out of the scope of this PR

This fixes #136

Author: jlbosse

Project.toml

@@ -1,7 +1,7 @@
 name = "SkewLinearAlgebra"
 uuid = "5c889d49-8c60-4500-9d10-5d3a22e2f4b9"
 authors = ["smataigne <[email protected]> and contributors"]
-version = "0.1.1"
+version = "0.1.2"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

README.md

@@ -10,6 +10,6 @@ The `SkewLinearAlgebra` package provides specialized matrix types, optimized met
 
 In particular, it defines new `SkewHermitian` and `SkewHermTridiagonal` matrix types supporting optimized eigenvalue/eigenvector,
 Hessenberg factorization, and matrix exponential/trigonometric functions.  It also provides functions to compute the
-[Pfaffian](https://en.wikipedia.org/wiki/Pfaffian) of real skew-symmetric matrices, along with a Cholesky-like factorization.
+[Pfaffian](https://en.wikipedia.org/wiki/Pfaffian) of skew-symmetric matrices, along with a Cholesky-like factorization for real skew-symmetric matrices.
 
 See the [Documentation](https://julialinearalgebra.github.io/SkewLinearAlgebra.jl/dev/) for details.

docs/src/pfaffian.md

@@ -1,10 +1,10 @@
 # Pfaffian calculations
 
-A real skew-symmetrix matrix ``A = -A^T`` has a special property: its determinant is the square of a polynomial function of the
+A skew-symmetrix matrix ``A = -A^T`` has a special property: its determinant is the square of a polynomial function of the
 matrix entries, called the [Pfaffian](https://en.wikipedia.org/wiki/Pfaffian).   That is, ``\mathrm{det}(A) = \mathrm{Pf}(A)^2``, but
 knowing the Pfaffian itself (and its sign, which is lost in the determinant) is useful for a number of applications.
 
-We provide a function `pfaffian(A)` to compute the Pfaffian of a real skew-symmetric matrix `A`.
+We provide a function `pfaffian(A)` to compute the Pfaffian of a skew-symmetric matrix `A`.
 ```jl
 julia> using SkewLinearAlgebra, LinearAlgebra
 
@@ -49,7 +49,7 @@ julia> pfaffian!(BigInt[0 2 -7; -2 0 -8; 7 8 0])
 ```
 (Note that the Pfaffian is *always zero* for any *odd* size skew-symmetric matrix.)
 
-Since the computation of the pfaffian can easily overflow/underflow the maximum/minimum representable floating-point value, we also provide a function `logabspfaffian` (along with an in-place variant `logabspfaffian!`) that returns a tuple `(logpf, sign)` such
+Since the computation of the pfaffian can easily overflow/underflow the maximum/minimum representable floating-point value, we also provide a function `logabspfaffian` for real skew-symmetric matrices (along with an in-place variant `logabspfaffian!`) that returns a tuple `(logpf, sign)` such
 that the Pfaffian is `sign * exp(logpf)`.   (This is similar to the [`logabsdet`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.logabsdet) function in Julia's `LinearAlgebra` library to compute the log of the determinant.)
 
 ```jl
@@ -81,4 +81,4 @@ SkewLinearAlgebra.pfaffian
 SkewLinearAlgebra.pfaffian!
 SkewLinearAlgebra.logabspfaffian
 SkewLinearAlgebra.logabspfaffian!
-```
+```

src/SkewLinearAlgebra.jl

@@ -15,6 +15,7 @@ export
     JMatrix,
     #functions
     isskewhermitian,
+    isskewsymmetric,
     skewhermitian,
     skewhermitian!,
     pfaffian,

src/pfaffian.jl

@@ -51,7 +51,7 @@ end
 
 function exactpfaffian!(A::AbstractMatrix)
     LinearAlgebra.require_one_based_indexing(A)
-    isskewhermitian(A) || throw(ArgumentError("Pfaffian requires a skew-Hermitian matrix"))
+    isskewsymmetric(A) || throw(ArgumentError("Pfaffian requires a skew-symmetric matrix"))
     return _exactpfaffian!(A)
 end
 exactpfaffian!(A::SkewHermitian) = _exactpfaffian!(A.data)
@@ -87,28 +87,77 @@ function _pfaffian!(A::SkewHermTridiagonal{<:Real})
     return pf
 end
 
-pfaffian!(A::Union{SkewHermitian{<:Real},SkewHermTridiagonal{<:Real}})= _pfaffian!(A)
-pfaffian(A::Union{SkewHermitian{<:Real},SkewHermTridiagonal{<:Real}})= pfaffian!(copyeigtype(A))
+# Using the Parley-Reid algorithm from https://arxiv.org/abs/1102.3440 as implented
+# in https://github.com/KskAdch/TopologicalNumbers.jl/blob/42bda634d47aecb67cfcd495d97e0723b0e73a4f/src/pfaffian.jl#L332
+function _pfaffian!(A::AbstractMatrix{<:Complex})
+    n = size(A, 1)
+    isodd(n) && return zero(eltype(A))
+    tau = Array{eltype(A)}(undef, n - 2)
+    pf = one(eltype(A))
+    @inbounds for k in 1:2:n-1
+        tauk = @view tau[k:end]
+
+        # Pivot if neccessary
+        @views kp = k + argmax(abs.(A[k+1:end, k]))
+        if kp != k + 1
+            @inbounds @simd for l in k:n
+                A[k+1,l], A[kp,l] = A[kp,l], A[k+1,l]
+            end
+            @inbounds @simd for l in k:n
+                A[l,k+1], A[l,kp] = A[l,kp], A[l,k+1]
+            end
+            pf *= -1
+        end
+
+        # Apply Gauss transformation and update `pf`
+        @inbounds if A[k+1,k] != zero(eltype(A))
+            @inbounds @views tauk .= A[k,k+2:end] ./ A[k,k+1]
+            pf *= @inbounds A[k,k+1]
+            if k + 2 <= n
+                @inbounds for l1 in eachindex(tauk)
+                    @simd for l2 in eachindex(tauk)
+                        @fastmath A[k+1+l2,k+1+l1] +=
+                            tauk[l2] * A[k+1+l1,k+1] -
+                            tauk[l1] * A[k+1+l2,k+1]
+                    end
+                end
+            end
+        else
+            return zero(eltype(A)) # Pfaffian is zero if there is a zero on the super/subdiagonal
+        end
+    end
+
+    return pf
+end
+
+pfaffian!(A::Union{SkewHermitian{<:Real},SkewHermTridiagonal{<:Real}}) = _pfaffian!(A)
 
 """
     pfaffian(A)
 
-Returns the pfaffian of `A` where a is a real skew-Hermitian matrix.
-If `A` is not of type `SkewHermitian{<:Real}`, then `isskewhermitian(A)`
-is checked to ensure that `A == -A'`
+Returns the pfaffian of `A` where a is a  skew-symmetric matrix.
+If `A` is not of type `SkewHermitian{<:Real}`, then `isskewsymmetric(A)`
+is checked to ensure that `A == -transpose(A)`
 """
-pfaffian(A::AbstractMatrix{<:Real}) = pfaffian!(copyeigtype(A))
+pfaffian(A::AbstractMatrix) = pfaffian!(copyeigtype(A))
 
 """
     pfaffian!(A)
 
 Similar to [`pfaffian`](@ref), but overwrites `A` in-place with intermediate calculations.
 """
 function pfaffian!(A::AbstractMatrix{<:Real})
-    isskewhermitian(A) || throw(ArgumentError("Pfaffian requires a skew-Hermitian matrix"))
+    isskewhermitian(A) || throw(ArgumentError("Pfaffian requires a skew-symmetric matrix"))
     return _pfaffian!(SkewHermitian(A))
 end
 
+function pfaffian!(A::AbstractMatrix{<:Complex})
+    LinearAlgebra.require_one_based_indexing(A)
+    isskewsymmetric(A) || throw(ArgumentError("Pfaffian requires a skew-symmetric matrix"))
+    return _pfaffian!(A)
+end
+
+
 function _logabspfaffian!(A::SkewHermitian{<:Real})
     n = size(A, 1)
     isodd(n) && return convert(eltype(A), -Inf), zero(eltype(A))

src/skewhermitian.jl

@@ -82,7 +82,7 @@ Base.size(A::SkewHermitian) = size(A.data)
 Returns whether `A` is skew-Hermitian, i.e. whether `A == -A'`.
 """
 function isskewhermitian(A::AbstractMatrix{<:Number})
-    axes(A,1) == axes(A,2) || throw(ArgumentError("axes $(axes(A,1)) and $(axex(A,2)) do not match"))
+    axes(A,1) == axes(A,2) || throw(ArgumentError("axes $(axes(A,1)) and $(axes(A,2)) do not match"))
     @inbounds for i in axes(A,1)
         for j = firstindex(A, 1):i
             A[i,j] == -A[j,i]' || return false
@@ -93,6 +93,23 @@ end
 isskewhermitian(A::SkewHermitian) = true
 isskewhermitian(a::Number) = a == -a'
 
+"""
+    isskewsymmetric(A)
+
+Returns whether `A` is skew-symmetric, i.e. whether `A == -transpose(A)`.
+"""
+function isskewsymmetric(A::AbstractMatrix{<:Number})
+    axes(A,1) == axes(A,2) || throw(ArgumentError("axes $(axes(A,1)) and $(axes(A,2)) do not match"))
+    @inbounds for i in axes(A,1)
+        for j = firstindex(A, 1):i
+            A[i,j] == -A[j,i] || return false
+        end
+    end
+    return true
+end
+isskewsymmetric(A::SkewHermitian{<:Real}) = true
+isskewsymmetric(a::Number) = iszero(a)
+
 """
     skewhermitian!(A)
 

test/runtests.jl

@@ -386,6 +386,7 @@ end
 end
 
 @testset "pfaffian.jl" begin
+    # real skew-hermitian matrices
     for n in [1, 2, 3, 10, 11]
         A = skewhermitian(rand(-10:10,n,n) * 2)
         Abig = BigInt.(A.data)
@@ -401,6 +402,20 @@ end
         logpf, sign = logabspfaffian(S)
         @test pfaffian(S) ≈ sign * exp(logpf) ≈ sign * sqrt(det(Matrix(S)))
     end
+
+    # complex skew-symmetric matrices
+    # n=11 is a bit jinxed because of https://github.com/JuliaLang/julia/issues/54287
+    for n in [1, 2, 3, 10, 20]
+        A = rand((-10:10) .+ 1im * (-10:10)', n, n)
+        A = (A .- transpose(A)) ./ 2
+        Abig = Complex{Rational{BigInt}}.(A)
+        @test pfaffian(A) ≈ SkewLinearAlgebra.exactpfaffian(Abig)
+        @test pfaffian(A)^2 ≈ det(A) atol=√eps(Float64) * max(1, abs(det(A)))
+        if VERSION ≥ v"1.7" # for exact det of BigInt matrices
+            @test SkewLinearAlgebra.exactpfaffian(Abig)^2 == det(Abig)
+        end
+    end
+
     # issue #49
     @test pfaffian(big.([0 14 7 -10 0 10 0 -11; -14 0 -10 7 13 -9 -12 -13; -7 10 0 -4 6 -17 -1 18; 10 -7 4 0 -2 -4 0 11; 0 -13 -6 2 0 -8 -18 17; -10 9 17 4 8 0 -8 12; 0 12 1 0 18 8 0 0; 11 13 -18 -11 -17 -12 0 0])) == -119000