"""A package for computing Pfaffians"""
import numpy as np
import scipy.linalg as la
import scipy.sparse as sp
import math, cmath
def householder_real(x):
"""(v, tau, alpha) = householder_real(x)
Compute a Householder transformation such that
(1-tau v v^T) x = alpha e_1
where x and v a real vectors, tau is 0 or 2, and
alpha a real number (e_1 is the first unit vector)
"""
assert x.shape[0]>0
sigma=np.dot(x[1:],x[1:])
if sigma==0:
return (np.zeros(x.shape[0]), 0, x[0])
else:
norm_x=math.sqrt(x[0]**2+sigma)
v=x.copy();
if x[0]<=0:
v[0]-=norm_x
alpha=+norm_x
else:
v[0]+=norm_x
alpha=-norm_x
v/=np.linalg.norm(v)
return (v, 2, alpha)
def householder_complex(x):
"""(v, tau, alpha) = householder_real(x)
Compute a Householder transformation such that
(1-tau v v^T) x = alpha e_1
where x and v a complex vectors, tau is 0 or 2, and
alpha a complex number (e_1 is the first unit vector)
"""
assert x.shape[0]>0
sigma=np.dot(np.conj(x[1:]), x[1:])
if sigma==0:
return (np.zeros(x.shape[0]), 0, x[0])
else:
norm_x=cmath.sqrt(x[0].conjugate()*x[0]+sigma)
v=x.copy();
phase=cmath.exp(1j*math.atan2(x[0].imag, x[0].real))
v[0]+=phase*norm_x
v/=np.linalg.norm(v)
return (v, 2, -phase*norm_x)
def skew_tridiagonalize(A, overwrite_a=False, calc_q=True):
""" T, Q = skew_tridiagonalize(A, overwrite_a, calc_q=True)
or
T = skew_tridiagonalize(A, overwrite_a, calc_q=False)
Bring a real or complex skew-symmetric matrix (A=-A^T) into
tridiagonal form T (with zero diagonal) with a orthogonal
(real case) or unitary (complex case) matrix U such that
A = Q T Q^T
(Note that Q^T and *not* Q^dagger also in the complex case)
A is overwritten if overwrite_a=True (default: False), and
Q only calculated if calc_q=True (default: True)
"""
assert A.shape[0] == A.shape[1] > 0
assert abs((A+A.T).max())<1e-14
n = A.shape[0]
A = np.asarray(A)
if np.issubdtype(A.dtype, np.complexfloating):
householder = householder_complex
elif not np.issubdtype(A.dtype, np.number):
raise TypeError("pfaffian() can only work on numeric input")
else:
householder = householder_real
if not overwrite_a:
A = A.copy()
if calc_q:
Q = np.eye(A.shape[0], dtype=A.dtype)
for i in range(A.shape[0]-2):
v, tau, alpha = householder(A[i+1:,i])
A[i+1, i] = alpha
A[i, i+1] = -alpha
A[i+2:, i] = 0
A[i, i+2:] = 0
w = tau*np.dot(A[i+1:, i+1:], v.conj());
A[i+1:,i+1:]+=np.outer(v,w)-np.outer(w,v)
if calc_q:
y = tau*np.dot(Q[:, i+1:], v)
Q[:, i+1:]-=np.outer(y,v.conj())
if calc_q:
return (np.asmatrix(A), np.asmatrix(Q))
else:
return np.asmatrix(A)
def skew_LTL(A, overwrite_a=False, calc_L=True, calc_P=True):
""" T, L, P = skew_LTL(A, overwrite_a, calc_q=True)
Bring a real or complex skew-symmetric matrix (A=-A^T) into
tridiagonal form T (with zero diagonal) with a lower unit
triangular matrix L such that
P A P^T= L T L^T
A is overwritten if overwrite_a=True (default: False),
L and P only calculated if calc_L=True or calc_P=True,
respectively (default: True).
"""
assert A.shape[0] == A.shape[1] > 0
assert abs((A+A.T).max())<1e-14
n = A.shape[0]
A = np.asarray(A)
if not overwrite_a:
A = A.copy()
if calc_L:
L = np.eye(n, dtype=A.dtype)
if calc_P:
Pv = np.arange(n)
for k in range(n-2):
kp = k+1+np.abs(A[k+1:,k]).argmax()
if kp != k+1:
temp = A[k+1,k:].copy()
A[k+1,k:] = A[kp,k:]
A[kp,k:] = temp
temp = A[k:,k+1].copy()
A[k:,k+1] = A[k:,kp]
A[k:,kp] = temp
if calc_L:
temp = L[k+1,1:k+1].copy()
L[k+1,1:k+1] = L[kp,1:k+1]
L[kp,1:k+1] = temp
if calc_P:
temp = Pv[k+1]
Pv[k+1] = Pv[kp]
Pv[kp] = temp
if A[k+1,k] != 0.0:
tau = A[k+2:,k].copy()
tau /= A[k+1,k]
A[k+2:,k] = 0.0
A[k,k+2:] = 0.0
A[k+2:,k+2:] += np.outer(tau, A[k+2:,k+1])
A[k+2:,k+2:] -= np.outer(A[k+2:,k+1], tau)
if calc_L:
L[k+2:,k+1] = tau
if calc_P:
P = sp.csr_matrix( (np.ones(n), (np.arange(n), Pv)) )
if calc_L:
if calc_P:
return (np.asmatrix(A), np.asmatrix(L), P)
else:
return (np.asmatrix(A), np.asmatrix(L))
else:
if calc_P:
return (np.asmatrix(A), P)
else:
return np.asmatrix(A)
def pfaffian(A, overwrite_a=False, method='P'):
""" pfaffian(A, overwrite_a=False, method='P')
Compute the Pfaffian of a real or complex skew-symmetric
matrix A (A=-A^T). If overwrite_a=True, the matrix A
is overwritten in the process. This function uses
either the Parlett-Reid algorithm (method='P', default),
or the Householder tridiagonalization (method='H')
"""
assert A.shape[0] == A.shape[1] > 0
assert abs((A+A.T).max()) < 1e-14, abs((A+A.T).max())
assert method == 'P' or method == 'H'
if method == 'P':
return pfaffian_LTL(A, overwrite_a)
else:
return pfaffian_householder(A, overwrite_a)
def pfaffian_LTL(A, overwrite_a=False):
""" pfaffian_LTL(A, overwrite_a=False)
Compute the Pfaffian of a real or complex skew-symmetric
matrix A (A=-A^T). If overwrite_a=True, the matrix A
is overwritten in the process. This function uses
the Parlett-Reid algorithm.
"""
assert A.shape[0] == A.shape[1] > 0
assert abs((A+A.T).max())<1e-14
n = A.shape[0]
A = np.asarray(A)
if n%2==1:
return 0
if not overwrite_a:
A = A.copy()
pfaffian_val = 1.0
for k in range(0, n-1, 2):
kp = k+1+np.abs(A[k+1:,k]).argmax()
if kp != k+1:
temp = A[k+1,k:].copy()
A[k+1,k:] = A[kp,k:]
A[kp,k:] = temp
temp = A[k:,k+1].copy()
A[k:,k+1] = A[k:,kp]
A[k:,kp] = temp
pfaffian_val *= -1
if A[k+1,k] != 0.0:
tau = A[k,k+2:].copy()
tau /= A[k,k+1]
pfaffian_val *= A[k,k+1]
if k+2<n:
A[k+2:,k+2:] += np.outer(tau, A[k+2:,k+1])
A[k+2:,k+2:] -= np.outer(A[k+2:,k+1], tau)
else:
return 0.0
return pfaffian_val
def pfaffian_householder(A, overwrite_a=False):
""" pfaffian(A, overwrite_a=False)
Compute the Pfaffian of a real or complex skew-symmetric
matrix A (A=-A^T). If overwrite_a=True, the matrix A
is overwritten in the process. This function uses the
Householder tridiagonalization.
Note that the function pfaffian_schur() can also be used in the
real case. That function does not make use of the skew-symmetry
and is only slightly slower than pfaffian_householder().
"""
assert A.shape[0] == A.shape[1] > 0
assert abs((A+A.T).max())<1e-14
n = A.shape[0]
if n%2==1:
return 0
if np.issubdtype(A.dtype, np.complexfloating):
householder=householder_complex
elif not np.issubdtype(A.dtype, np.number):
raise TypeError("pfaffian() can only work on numeric input")
else:
householder=householder_real
A = np.asarray(A)
if not overwrite_a:
A = A.copy()
pfaffian_val = 1.
for i in range(A.shape[0]-2):
v, tau, alpha = householder(A[i+1:,i])
A[i+1, i] = alpha
A[i, i+1] = -alpha
A[i+2:, i] = 0
A[i, i+2:] = 0
w = tau*np.dot(A[i+1:, i+1:], v.conj());
A[i+1:,i+1:]+=np.outer(v,w)-np.outer(w,v)
if tau!=0:
pfaffian_val *= 1-tau
if i%2==0:
pfaffian_val *= -alpha
pfaffian_val *= A[n-2,n-1]
return pfaffian_val
def pfaffian_schur(A, overwrite_a=False):
"""Calculate Pfaffian of a real antisymmetric matrix using
the Schur decomposition. (Hessenberg would in principle be faster,
but scipy-0.8 messed up the performance for scipy.linalg.hessenberg()).
This function does not make use of the skew-symmetry of the matrix A,
but uses a LAPACK routine that is coded in FORTRAN and hence faster
than python. As a consequence, pfaffian_schur is only slightly slower
than pfaffian().
"""
assert np.issubdtype(A.dtype, np.number) and not np.issubdtype(A.dtype, np.complexfloating)
assert A.shape[0] == A.shape[1] > 0
assert abs(A + A.T).max() < 1e-14
if A.shape[0]%2 == 1:
return 0
(t, z) = la.schur(A, output='real', overwrite_a=overwrite_a)
l = np.diag(t, 1)
return np.prod(l[::2]) * la.det(z)