LazyInverses.jl

LazyInverses provides a lazy wrapper for a matrix inverse, akin to Adjoint in Julia Base. See the README for example use cases.

A current highlight is the package's implementation of a type of energetic inner product (ternary dot product) for which a 7-fold increase in performance is observed compared to a naïve implementation.

Installation

Simply type ] follwed by add LazyInverses in the REPL.

Basic Usage

The package exports two types Inverse and PseudoInverse, and their corresponding "smart" constructors inverse and pinverse, which return the lazy wrappers, unless the input is a number or a 1 x 1 matrix, in which case the inverse is returned directly. This first example highlights the lazy behavior of inverse, and contrasts it to inv:

using LazyInverses
using LinearAlgebra

n = 1024
x = randn(n)
A = randn(n, n)
@time B = inv(A)
  0.058023 seconds (5 allocations: 8.508 MiB)
@time C = inverse(A)
  0.000003 seconds (1 allocation: 16 bytes)
A \ x ≈ C * x
  true
A * x ≈ C \ x
  true

Use Cases

Mahalanobis Norms and Energetic Inner Products

In a number of important models, including any that rely on computations with a multivariate Normal distribution, it is necessary to compute a Mahalanobis norm with an inverse of a positive semi-definite matrix. In order to allow for a particularly efficient implementation of this operation, LazyInverses.jl extends the ternary dot product like so:

dot(::AbstractVector, ::Inverse{<:Number, <:Union{<:Cholesky, <:CholeskyPivoted}}, ::AbstractVector)

Benchmarking the package's implementation against a naïve implementation, we observe an up to 7-fold increase in performance.

using LazyInverses
using LinearAlgebra
using BenchmarkTools

n = 1024
x = randn(n)
y = randn(n)
A = randn(n, n)
A = A'A
C = cholesky(A)
invC = inverse(C)

println("ternary dot-product multiplication")
@btime dot($x, $invC, $x)
  121.794 μs (3 allocations: 8.16 KiB)
@btime dot($x, $invC, $y)
  172.736 μs (15 allocations: 17.06 KiB)
@btime dot($x, $C \ $x)
  855.008 μs (1 allocation: 8.12 KiB)
@btime dot($x, $C \ $y)
  850.760 μs (1 allocation: 8.12 KiB)

Further, we observe a speed-up of up to a factor of 2 for a ternary matrix multiplication, where the middle matrix is an inverse Cholesky factorization.

k = n
X = randn(k, n)
Y = randn(n, k)
println("ternary multiplication")
@btime *($X, $invC, $X')
  32.755 ms (5 allocations: 16.00 MiB)
@btime *($X, $invC, $Y)
  54.706 ms (4 allocations: 16.00 MiB)
@btime *($X, $C \ $X')
  59.221 ms (4 allocations: 16.00 MiB)
@btime *($X, $C \ $Y)
  56.177 ms (4 allocations: 16.00 MiB)

Notably, the implementation takes advantage of threaded parallelism starting at a vector dimension of 1024, but calls a single-threaded implementation for smaller vectors to minimize constants and maximize performance.

WoodburyIdentity.jl

Coming soon.

KroneckerProducts.jl

Coming soon.

Zygote.jl

Coming soon.

Notes

The experiments were run on a 2017 MacBook Pro with an i7 dual core and 16 GB of RAM.