PaperlayerView exact revision on arXiv

Source-linked AI summary

Efficient numerical computation of the Pfaffian for dense and banded skew-symmetric matrices

M. Wimmer

arXiv:1102.3440v2cond-mat.mes-hallcs.MSmath.NAphysics.comp-ph

TL;DR

The paper addresses efficient computation of Pfaffians and canonical forms for skew-symmetric matrices, where tridiagonal reduction is central. It develops structure-aware Gauss, Householder, and Givens algorithms for dense and banded matrices, with implementations across several languages. The methods benchmark favorably and reproduce equivalent topological-charge definitions for a disordered nanowire.

  • Problem

    Existing approaches do not comprehensively exploit skew-symmetry, level-3 operations, complex arithmetic, or banded sparsity for tridiagonalization.

  • Method

    The paper develops Gauss-based Parlett-Reid algorithms and unitary Householder/Givens transformations to reduce dense or banded skew-symmetric matrices to tridiagonal form.

  • Results

    The implementations are faster than competing approaches and numerically show equivalence between Hamiltonian- and scattering-matrix definitions of topological charge.

  • Takeaways & Limitations

    The methods enable Pfaffian-based computation for large disordered nanowire systems by exploiting their sparse banded structure.

Abstract

from arXiv · show

Computing the Pfaffian of a skew-symmetric matrix is a problem that arises in various fields of physics. Both computing the Pfaffian and a related problem, computing the canonical form of a skew-symmetric matrix under unitary congruence, can be solved easily once the skew-symmetric matrix has been reduced to skew-symmetric tridiagonal form. We develop efficient numerical methods for computing this tridiagonal form based on Gauss transformations, using a skew-symmetric, blocked form of the Parlett-Reid algorithm, or based on unitary transformations, using block Householder transformations and Givens rotations, that are applicable to dense and banded matrices, respectively. We also give a complete and fully optimized implementation of these algorithms in Fortran, and also provide Python, Matlab and Mathematica implementations for convenience. Finally, we apply these methods to compute the topological charge of a class D nanowire, and show numerically the equivalence of definitions based on the Hamiltonian and the scattering matrix.

I. INTRODUCTION

The paper motivates efficient Pfaffian computation by reducing skew-symmetric matrices to tridiagonal form, which also supports canonical-form computation. Existing direct evaluation is prohibitively expensive for larger matrices, while Gaussian and unitary approaches offer O(n^3) alternatives.

  • A skew-symmetric matrix satisfies A = −A^T, and its determinant is the square of its Pfaffian.
  • The Pfaffian has applications in topological charges, quantum Monte Carlo, Ising spin glasses, fractional quantum Hall states, quantum chaos, and lattice quantum field theory.
  • Direct evaluation of the permutation definition costs O(n!), making it impractical for larger matrices.
  • Reducing a skew-symmetric matrix to tridiagonal or partial tridiagonal form enables straightforward Pfaffian evaluation.
  • Skew-symmetric Gaussian elimination computes an LTLT factorization in O(n^3), while unitary transformations provide pivot-free stable computation for dense matrices and can exploit bandedness.

B. Tridiagonalization and the canonical form of skew-symmetric matrices

Tridiagonalization serves both efficient Pfaffian evaluation and canonical-form computation under unitary congruence. The canonical form connects to singular values and eigenvalue degeneracies relevant to physical applications.

  • Skew-symmetric tridiagonalization is the common computational route for both Pfaffian calculation and canonical-form construction.
  • Every skew-symmetric matrix admits a canonical form A = UΞU^T under unitary congruence, with 2 × 2 blocks and possible zero blocks.
  • This canonical form has been used to establish degeneracies of transmission and Andreev-reflection eigenvalues.
  • For even-dimensional skew-symmetric tridiagonal matrices, the canonical parameters are given by the nonzero singular values of an associated bidiagonal matrix.
  • In the real case, eigenvalues occur as ±iσ_j, while the complex matrix A* A has doubly degenerate eigenvalues σ^2.

C. Skew-symmetric tridiagonalization and existing approaches

The paper identifies gaps in existing tridiagonalization software: limited support for complex matrices, insufficient exploitation of skew-symmetry and sparsity, and a lack of efficient blocked implementations. It therefore targets real or complex dense and banded matrices.

  • Existing approaches connect Pfaffian and canonical-form computation to skew-symmetric tridiagonalization, but no comprehensive implementation addresses all cases.
  • Standard real decompositions do not exploit skew-symmetry for precision and speed and cannot handle complex skew-symmetric matrices.
  • A structure-exploiting implementation exists for real dense matrices, but its complex counterpart is unavailable.
  • Existing libraries and accompanying code lack efficiency, transformation-matrix access, or full exploitation of skew-symmetry.
  • Existing approaches generally omit level-3 block algorithms and the sparsity of banded matrices, despite their importance for memory access and practical structure.
  • The paper develops algorithms for tridiagonalizing real or complex skew-symmetric matrices while exploiting skew-symmetry and, when present, bandedness.

D. Outline

The paper outlines Gauss-, Householder-, and Givens-based reductions for dense and banded skew-symmetric matrices, emphasizing transformations that preserve structure and support efficient Pfaffian computation. The dense Parlett-Reid route uses pivoting and can be blocked for level-3 operations.

  • D. Outline: The paper considers dense and banded tridiagonalization using Gauss transformations, Householder reflections, and Givens rotations.
  • Gauss transformations: Gauss transformations eliminate selected row and column entries while permutations place a maximal nonzero pivot in the active position.
  • Gauss transformations: Applying transformations from both sides preserves skew-symmetry and progressively reduces leading rows and columns to tridiagonal form.
  • Blocked algorithms: Regrouping updates enables level-3 matrix operations with more favorable memory access patterns in blocked Parlett-Reid and Householder algorithms.
  • Cost and limitation: The full skew-symmetric LTLT decomposition costs 2n^3/3 flops, while partial tridiagonalization and Pfaffian computation require n^3/3 flops.
  • Cost and limitation: Although Aasen’s and Bunch-Kaufmann algorithms achieve n^3/3 flops for full decompositions, adapting them to the partial factorization needed for Pfaffians remains difficult.

C. Tridiagonalization of dense matrices with Householder reflections

Householder reflections reduce dense skew-symmetric matrices to tridiagonal form without pivoting, supporting both Pfaffian computation and canonical-form calculations. The method costs 4n^3/3 flops for complete tridiagonalization, but only 2n^3/3 for the partial form needed by the Pfaffian, while remaining twice as costly as Parlett-Reid.

  • Householder method: Householder transformations eliminate all entries below a chosen point in each column or row without pivoting, preserving the skew-symmetric tridiagonalization structure.The transformations are unitary and can be computed to machine precision; complex choices can produce a purely real tridiagonal matrix.
  • Householder method: Each step applies a block Householder transformation and updates the trailing matrix through a full matrix-vector product and outer-product operation.The vector w is formed from the trailing matrix and Householder vector, and the outer-product update costs 4(n −k)^2 flops per step.
  • Computational cost: Complete tridiagonalization requires 4n^3/3 flops, while Pfaffian computation requires only 2n^3/3 flops using a partial tridiagonal form.The reduced cost comes from skipping every other row/column elimination.
  • Transformation determinant: The determinant of the accumulated Householder transformation can be computed in O(n^2), making it negligible compared with tridiagonalization.For τ = 2/v†v, each Householder reflection has determinant −1; other choices require separate determinant evaluation.
  • Comparison: Householder tridiagonalization is twice as costly as Parlett-Reid for Pfaffian computation, but remains useful for obtaining the canonical form.Its connection to unitary-congruence canonical-form calculations gives it a role beyond Pfaffian evaluation.

D. Tridiagonalization of band matrices with Givens rotations

Dense algorithms do not preserve finite bandwidth: pivoting can expand the band unpredictably, while Householder outer-product updates create out-of-band entries. Givens rotations provide a selective band-preserving strategy by chasing limited fill-in toward the matrix boundary, with cost O(bn^2).

  • Motivation: Symmetric pivoting in Parlett-Reid can cause uncontrolled bandwidth growth, while Householder outer-product updates rapidly introduce values outside the band.These limitations motivate a separate unitary-transformation approach for banded skew-symmetric matrices.
  • Givens rotations: Givens rotations G_i,i+1 selectively modify adjacent rows and columns to reduce a banded skew-symmetric matrix to tridiagonal form.The rotation is applied from the left and its transpose from the right, preserving skew-symmetry.
  • Fill-in management: Each rotation introduces at most one additional out-of-band nonzero per affected row and column, which can be chased beyond the matrix boundary.Successive Givens transformations move the fill-in farther down the band.
  • Skew-symmetric update: The diagonal 2 × 2 blocks remain invariant under the paired left-right Givens transformations because of skew-symmetry.This is the principal update difference from the symmetric or Hermitian tridiagonalization routine.
  • Complexity: Banded tridiagonalization of an n × n matrix with bandwidth b scales as O(bn^2).The implementation also includes a routine for this banded reduction.

A. Fortran

The authors provide a Fortran implementation for dense and banded skew-symmetric tridiagonalization and Pfaffian computation, with LAPACK-compatible design and auxiliary interfaces for other languages. Storage schemes explicitly represent only the relevant triangular or banded entries, and benchmark tables compare runtimes across methods and architectures.

  • Storage: Dense matrices use ordinary two-dimensional storage with only one strictly triangular part set, while banded matrices use a K × N array for nonzero off-diagonals.The banded scheme explicitly stores the zero diagonal to match LAPACK’s symmetric and Hermitian band-storage design.
  • Storage: Upper-triangular band storage places A_i,j at AB(K + 1 + i −j, j), whereas lower-triangular storage places it at AB(1 + i −j, j).The index ranges differ according to whether the upper or lower triangle is stored.
  • Routine library: The Fortran library includes dense and banded tridiagonalization routines plus Pfaffian routines built on those reductions.Dense routines cover Parlett-Reid and Householder methods, while SKBTRD and SKBPFA address banded matrices.
  • Optimization: Blocked routines expose an internal block-size parameter that users can change to optimize performance for a specific architecture.Default block sizes match the corresponding LAPACK symmetric routines.
  • Dependencies: The implementation uses LAPACK and BLAS transformation routines and therefore requires both libraries to be linked.LAPACK and BLAS provide optimized Householder and Givens operations across common computer architectures.
  • Benchmarks: Benchmark Table II reports Pfaffian-computation times in seconds for multiple methods on two architectures, with banded cases indexed by strictly upper or lower bandwidth k.For banded matrices, the full bandwidth is 2k + 1.
  • Language interfaces: Python, Matlab, and Mathematica versions provide convenient standalone dense Householder tridiagonalization implementations, although they are slower than Fortran.The alternative-language versions are intended especially for users without Fortran access.

A. Benchmarks

Benchmarks show that blocked implementations improve Pfaffian and tridiagonalization performance, while banded performance depends on algorithm and matrix structure. The methods also enable large disordered-nanowire calculations, where Hamiltonian- and reflection-matrix-based topological charges agree closely.

  • Dense matrices: Blocked algorithms are faster than unblocked versions, with relative speed-ups reaching 60% depending on architecture.The Parlett-Reid speed-up exceeds that of Householder tridiagonalization because it contains more level-3 operations.
  • Banded matrices: For banded random matrices, Parlett-Reid performs well despite not explicitly exploiting bandedness.Its skew-symmetric outer-product update skips zeros in update vectors, whereas Householder suffers faster bandwidth growth for these matrices.
  • Banded matrices: The banded Givens approach uses less memory and is considerably faster than Householder tridiagonalization for canonical-form computation.For the tested matrix sizes, it remains slightly slower than Parlett-Reid; memory, rather than computation time, typically limits tractable sizes.
  • Comparison with existing approaches: The routines are typically about 10 times faster than other approaches, including real Hessenberg reduction in the real case.The full expected factor-of-10 speed-up is not always reached because of LAPACK optimization, while one comparison implementation has unfavorable memory access.
  • Nanowire application: The methods apply Pfaffian evaluation to large disordered nanowires by exploiting the banded skew-symmetric structure of the Majorana-basis Hamiltonian.For a disordered supercell, the matrix size scales as WL, and the work reports this as the first application of the Hamiltonian-based definition to such large disordered systems.
  • Nanowire application: Hamiltonian-Pfaffian and reflection-matrix definitions of the nanowire topological charge agree well for clean and disordered systems.Both predict a vanished topological phase at the largest disorder for 10 < EF/Eso < 25; small transition-point differences are attributed to finite-size geometry effects, while insulating regimes agree fully.
  • Nanowire application: The Hamiltonian-based calculation scales as W^3L^2, whereas the transport-based definition scales as W^3L and is computationally more favorable.The numerical equivalence of the two definitions supports using the transport formulation when computational scaling is decisive.

V. CONCLUSIONS

The paper reduces Pfaffian and canonical-form computation to skew-symmetric tridiagonal form, develops efficient tridiagonalization algorithms, and demonstrates their numerical utility for dense, banded, and nanowire problems.

  • Skew-symmetric tridiagonalization makes both Pfaffian computation and canonical-form computation straightforward.
  • Gauss-based blocked Parlett-Reid and unitary block Householder/Givens algorithms provide tridiagonalization methods for dense and banded skew-symmetric matrices.
  • A comprehensive numerical library implements the algorithms, with benchmark performance reported as superior to other approaches.
  • Applying the Pfaffian method to a disordered nanowire showed numerical equivalence between Hamiltonian- and scattering-matrix-based topological-charge definitions.

Appendix A: The computation of the canonical form of a skew-symmetric matrix

The appendices derive canonical forms and block tridiagonalization procedures for skew-symmetric matrices, emphasizing unitary reductions, accumulated updates, and implementation-aware computational efficiency.

  • Appendix A: The computation of the canonical form of a skew-symmetric matrix: A 2n × 2n skew-symmetric matrix admits a unitary congruence reduction to tridiagonal form, which can be reorganized using a permutation and a canonical block structure.
  • Appendix A: The computation of the canonical form of a skew-symmetric matrix: The canonical form follows from the singular value decomposition of the associated matrix J, with positive singular values determined by its rank.
  • Appendix A: The computation of the canonical form of a skew-symmetric matrix: In practice, real SVD routines suffice because complex matrices can first be reduced to real tridiagonal form by unitary transformations.
  • Appendix B: Block versions of the Parlett-Reid and Householder tridiagonalization algorithms: Accumulating Gauss or Householder transformations enables block updates with a higher level-3 fraction than individual level-2 operations.
  • Appendix B: Block versions of the Parlett-Reid and Householder tridiagonalization algorithms: Both block algorithms express successive reductions through structured transformations and update operations involving vectors generated from the current matrix.
  • Appendix B: Block versions of the Parlett-Reid and Householder tridiagonalization algorithms: Panel processing computes successive Householder reflections locally, accumulates them, and applies the combined update to the remaining matrix.
  • Appendix B: Block versions of the Parlett-Reid and Householder tridiagonalization algorithms: Parlett-Reid blocking is dominated by level-3 outer-product updates, whereas Householder blocking retains an inherently level-2 matrix-vector multiplication.

Appendix C: Upper versus lower triangle storage in the Fortran implementation

The Fortran implementation exploits skew-symmetry and column-major storage by operating on one triangle and favoring column operations, with upper-triangle variants reversing traversal.

  • The Fortran implementation works with only the lower or upper triangle and primarily uses column operations to match column-major memory layout.
  • Upper-triangle implementations begin tridiagonalization in the last column rather than the first.
  • Partial tridiagonalization changes the zero pattern, but the associated row-column permutation has determinant one and therefore preserves the Pfaffian.
  • The Parlett-Reid algorithm computes a UTU^T decomposition with U upper unit triangular.
Loading 1102.3440v2…