LazyArrays.jl

LazyArrays.LazyArraysModule

LazyArrays.jl

DevTraviscodecov

Lazy arrays and linear algebra in Julia

This package supports lazy analogues of array operations like vcat, hcat, and multiplication. This helps with the implementation of matrix-free methods for iterative solvers.

The package has been designed with high-performance in mind, so should outperform the non-lazy analogues from Base for many operations like copyto! and broadcasting. Some operations will be inherently slower due to extra computation, like getindex. Please file an issue for any examples that are significantly slower than their the analogue in Base.

Lazy operations

To construct a lazy representation of a function call f(x,y,z...), use the command applied(f, x, y, z...). This will return an unmaterialized object typically of type Applied that represents the operation. To realize that object, call materialize, which will typically be equivalent to calling f(x,y,z...). A macro @~ is available as a shorthand:

julia> using LazyArrays, LinearAlgebra

julia> applied(exp, 1)
Applied(exp,1)

julia> materialize(applied(exp, 1))
2.718281828459045

julia> materialize(@~ exp(1))
2.718281828459045

julia> exp(1)
2.718281828459045

The benefit of lazy operations is that they can be materialized in-place, possible using simplifications. For example, it is possible to do BLAS-like Matrix-Vector operations of the form α*A*x + β*y as implemented in BLAS.gemv! using a lazy applied object:

julia> A = randn(5,5); b = randn(5); c = randn(5); d = similar(c);

julia> d .= @~ 2.0 * A * b + 3.0 * c # Calls gemv!
5-element Array{Float64,1}:
 -2.5366335879717514
 -5.305097174484744  
 -9.818431932350942  
  2.421562605495651  
  0.26792916096572983

julia> 2*(A*b) + 3c
5-element Array{Float64,1}:
 -2.5366335879717514
 -5.305097174484744  
 -9.818431932350942  
  2.421562605495651  
  0.26792916096572983

julia> function mymul(A, b, c, d) # need to put in function for benchmarking
       d .= @~ 2.0 * A * b + 3.0 * c
       end
mymul (generic function with 1 method)

julia> @btime mymul(A, b, c, d) # calls gemv!
  77.444 ns (0 allocations: 0 bytes)
5-element Array{Float64,1}:
 -2.5366335879717514
 -5.305097174484744  
 -9.818431932350942  
  2.421562605495651  
  0.26792916096572983

julia> @btime 2*(A*b) + 3c; # does not call gemv!
  241.659 ns (4 allocations: 512 bytes)

This also works for inverses, which lower to BLAS calls whenever possible:

julia> A = randn(5,5); b = randn(5); c = similar(b);

julia> c .= @~ A \ b
5-element Array{Float64,1}:
 -2.5366335879717514
 -5.305097174484744  
 -9.818431932350942  
  2.421562605495651  
  0.26792916096572983

Lazy arrays

Often we want lazy realizations of matrices, which are supported via ApplyArray. For example, the following creates a lazy matrix exponential:

julia> E = ApplyArray(exp, [1 2; 3 4])
2×2 ApplyArray{Float64,2,typeof(exp),Tuple{Array{Int64,2}}}:
  51.969   74.7366
 112.105  164.074 

A lazy matrix exponential is useful for, say, in-place matrix-exponetial*vector:

julia> b = Vector{Float64}(undef, 2); b .= @~ E*[4,4]
2-element Array{Float64,1}:
  506.8220830628333
 1104.7145995988594
 ```
 While this works, it is not actually optimised (yet). 

 Other options do have special implementations that make them fast. We
 now give some examples. 


### Concatenation

Lazy `vcat` and `hcat` allow for representing the concatenation of
vectors without actually allocating memory, and support a fast
`copyto!`  for allocation-free population of a vector.

julia julia> using BenchmarkTools

julia> A = ApplyArray(vcat,1:5,2:3) # allocation-free 7-element ApplyArray{Int64,1,typeof(vcat),Tuple{UnitRange{Int64},UnitRange{Int64}}}: 1 2 3 4 5 2 3

julia> Vector(A) == vcat(1:5, 2:3) true

julia> b = Array{Int}(undef, length(A)); @btime copyto!(b, A); 26.670 ns (0 allocations: 0 bytes)

julia> @btime vcat(1:5, 2:3); # takes twice as long due to memory creation 43.336 ns (1 allocation: 144 bytes)

Similar is the lazy analogue of `hcat`:

julia julia> A = ApplyArray(hcat, 1:3, randn(3,10)) 3×11 ApplyArray{Float64,2,typeof(hcat),Tuple{UnitRange{Int64},Array{Float64,2}}}: 1.0 1.16561 0.224871 -1.36416 -0.30675 0.103714 0.590141 0.982382 -1.50045 0.323747 -1.28173 2.0 1.04648 1.35506 -0.147157 0.995657 -0.616321 -0.128672 -0.671445 -0.563587 -0.268389 -1.71004 3.0 -0.433093 -0.325207 -1.38496 -0.391113 -0.0568739 -1.55796 -1.00747 0.473686 -1.2113 0.0119156

julia> Matrix(A) == hcat(A.args...) true

julia> B = Array{Float64}(undef, size(A)...); @btime copyto!(B, A); 109.625 ns (1 allocation: 32 bytes)

julia> @btime hcat(A.args...); # takes twice as long due to memory creation 274.620 ns (6 allocations: 560 bytes)




### Kronecker products

We can represent Kronecker products of arrays without constructing the full
array:

julia julia> A = randn(2,2); B = randn(3,3);

julia> K = ApplyArray(kron,A,B) 6×6 ApplyArray{Float64,2,typeof(kron),Tuple{Array{Float64,2},Array{Float64,2}}}: -1.08736 -0.19547 -0.132824 1.60531 0.288579 0.196093 0.353898 0.445557 -0.257776 -0.522472 -0.657791 0.380564 -0.723707 0.911737 -0.710378 1.06843 -1.34603 1.04876 1.40606 0.252761 0.171754 -0.403809 -0.0725908 -0.0493262 -0.457623 -0.576146 0.333329 0.131426 0.165464 -0.0957293 0.935821 -1.17896 0.918584 -0.26876 0.338588 -0.26381

julia> C = Matrix{Float64}(undef, 6, 6); @btime copyto!(C, K); 61.528 ns (0 allocations: 0 bytes)

julia> C == kron(A,B) true



## Broadcasting

Base includes a lazy broadcast object called `Broadcasting`, but this is
not a subtype of `AbstractArray`. Here we have `BroadcastArray` which replicates
the functionality of `Broadcasting` while supporting the array interface.

julia julia> A = randn(6,6);

julia> B = BroadcastArray(exp, A);

julia> Matrix(B) == exp.(A) true

julia> B = BroadcastArray(+, A, 2);

julia> B == A .+ 2 true

Such arrays can also be created using the macro `@~` which acts on ordinary 
broadcasting expressions combined with `LazyArray`:

julia julia> C = rand(1000)';

julia> D = LazyArray(@~ exp.(C))

julia> E = LazyArray(@~ @. 2 + log(C))

julia> @btime sum(LazyArray(@~ C .* C'); dims=1) # without @~, 1.438 ms (5 allocations: 7.64 MiB) 74.425 μs (7 allocations: 8.08 KiB) ```

LazyArrays.AppliedType
Applied(f, A...)

is a lazy version of f(A...) that can be manipulated or materialized in a non-standard manner.

LazyArrays.LazyArrayType
LazyArray(x::Applied) :: ApplyArray
LazyArray(x::Broadcasted) :: BroadcastArray

Wrap a lazy object that wraps a computation producing an array to an array.

LazyArrays.cacheMethod
cache(array::AbstractArray)

Caches the entries of an array.

LazyArrays.@~Macro
@~ expr

Macro for creating a Broadcasted or Applied object. Regular calls like f(args...) inside expr are replaced with applied(f, args...). Dotted-calls like f(args...) inside expr are replaced with broadcasted.(f, args...). Use LazyArray(@~ expr) if you need an array-based interface.

julia> @~ A .+ B ./ 2

julia> @~ @. A + B / 2

julia> @~ A * B + C