Direct Linear Solvers for Symmetric Positive Definite Toeplitz Matrices

[SebastianAment]

Overview

This PR adds implementations of the Durbin and Levinson algorithms that improve on the performance of StatsBase's implementations by a factor of two for larger matrices. The performance advantage starts to appear before n = 32 and grow with n. Experiments were run on a 2021 MacBook Pro with an M1 Pro.

Since this package specializes on Toeplitz matrices, I think the implementations are better incorporated here than at StatsBase. Further, this PR also adds Trench's algorithm for the efficient inversion of symmetric positive definite Toeplitz matrices. See below for performance benchmarks.

For a reference, see Golub and Van Loan, page 208.

Benchmarks

Durbin

n = 4096

durbin
  5.571 ms (2 allocations: 32.05 KiB)
durbin!
  5.568 ms (0 allocations: 0 bytes)
StatsBase.durbin
  11.043 ms (2 allocations: 32.05 KiB)
StatsBase.durbin!
  11.041 ms (0 allocations: 0 bytes)
canonical solve
  161.529 ms (8 allocations: 128.09 MiB)

n = 32

durbin
  716.241 ns (1 allocation: 336 bytes)
durbin!
  646.341 ns (0 allocations: 0 bytes)
StatsBase.durbin
  818.494 ns (1 allocation: 336 bytes)
StatsBase.durbin!
  794.708 ns (0 allocations: 0 bytes)
canonical solve
  54.375 μs (4 allocations: 9.11 KiB)

Levinson

n = 4096

levinson
  9.533 ms (4 allocations: 64.09 KiB)
levinson!
  9.531 ms (2 allocations: 32.05 KiB)
StatsBase.levinson
  23.107 ms (4 allocations: 64.09 KiB)
StatsBase.levinson!
  23.110 ms (0 allocations: 0 bytes)
canonical solve
  161.600 ms (6 allocations: 128.06 MiB)

n = 32

levinson
  905.273 ns (2 allocations: 672 bytes)
levinson!
  897.925 ns (1 allocation: 336 bytes)
StatsBase.levinson
  1.221 μs (2 allocations: 672 bytes)
StatsBase.levinson!
  1.175 μs (0 allocations: 0 bytes)
direct solve
  52.750 μs (3 allocations: 8.78 KiB)

Trench

n = 4096

trench
  18.469 ms (14 allocations: 128.06 MiB)
canonical inv
  499.677 ms (6 allocations: 130.03 MiB)

n = 32

trench
  1.279 μs (7 allocations: 8.94 KiB)
canonical inv
  59.542 μs (4 allocations: 24.58 KiB)

[coveralls]

Coverage increased (+0.2%) to 91.667% when pulling e2a1199dd0f2478a3c903d674f4a67c18bd9ffcd on SebastianAment:efficient-symmetric-solvers into 91c1c540c7e7cca4f4d87003ab7c5c3d5bb52e87 on JuliaMatrices:master.

[codecov[bot]]

Codecov Report

Base: 91.45% // Head: 91.66% // Increases project coverage by +0.21% :tada:

Coverage data is based on head (e2a1199) compared to base (91c1c54). Patch coverage: 90.65% of modified lines in pull request are covered.

Additional details and impacted files

@@            Coverage Diff             @@
##           master      #64      +/-   ##
==========================================
+ Coverage   91.45%   91.66%   +0.21%     
==========================================
  Files           2        3       +1     
  Lines         468      600     +132     
==========================================
+ Hits          428      550     +122     
- Misses         40       50      +10     

Impacted Files	Coverage Δ
src/ToeplitzMatrices.jl	90.37% <12.50%> (-1.19%)	:arrow_down:
src/directLinearSolvers.jl	96.96% <96.96%> (ø)
src/iterativeLinearSolvers.jl	91.66% <0.00%> (+0.65%)	:arrow_up:

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

:umbrella: View full report at Codecov.
:loudspeaker: Do you have feedback about the report comment? Let us know in this issue.

[SebastianAment]

Further notes:

• StatsBase.levinson! mutates b, which is unexpected (and should probably be called levinson!!). For this reason, the new levinson! implementation still allocates a single vector of length n.
• If the current maintainers want to incorporate the new implementations, I'd also change lines 687 to 701 to use the new levinson implementation.
• We could also add a SymmetricPositiveDefiniteToeplitz type, for which inv and \ automatically call trench and levinson, respectively.

[dlfivefifty]

• We could also add a SymmetricPositiveDefiniteToeplitz

Not a good idea: there’s no equivalent in LinearAlgebra. The analogue there is explicitly calling cholesky

[SebastianAment]

• We could also add a SymmetricPositiveDefiniteToeplitz

Not a good idea: there’s no equivalent in LinearAlgebra. The analogue there is explicitly calling cholesky

This is not central to this PR, just adding my thoughts:

• There are the Symmetric and Hermitian wrappers in LinearAlgebra which do impact how factorize behaves. There's no clear reason in my mind why one property should be dispatched on, while the other one shouldn't. In the case of Toeplitz matrices, checking definiteness is also an O(n^2) operation, same as checking symmetry for regular matrices.
• Signaling positive definiteness via types would make it easier to write generic code that should work efficiently for both definite and indefinite systems.
• It would be easier for non-experts (i.e. people who have not heard of Trench's or Levinson's algorithms) to take advantage of the efficient solver and inversion algorithms by just turning on a flag. This could be done by having a PositiveDefinite trait as a field in the SymmetricToeplitz structure.

[SebastianAment]

I found a few minor style issues.

I just committed your suggestions.

[ParadaCarleton]

Is anything holding this up? @dkarrasch @SebastianAment

[dkarrasch]

Out of curiosity, why not make the overloaded functions available in StatsBase.jl? Do you need special fields of structs that are owned by this package here?

[ParadaCarleton]

@dkarrasch if the only thing missing is a couple style changes, perhaps we should just apply the changes and merge, to avoid holding this up too long?

[dkarrasch]

Out of curiosity, why not make the overloaded functions available in StatsBase.jl? Do you need special fields of structs that are owned by this package here?

I still wonder... There will be durbin methods and others both in StatsBase.jl and here. How will people judge which one to use? The backend appears to be fairly generic.

[ParadaCarleton]

Out of curiosity, why not make the overloaded functions available in StatsBase.jl? Do you need special fields of structs that are owned by this package here?

I still wonder... There will be durbin methods and others both in StatsBase.jl and here. How will people judge which one to use? The backend appears to be fairly generic.

Presumably they will be, but that would require a PR in StatsBase.jl.

[ParadaCarleton]

@dlfivefifty could you merge this?