Lazy arrays and linear algebra in Julia
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.
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
Appliedthat 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 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 + β*yas implemented in
BLAS.gemv!using a lazy applied object: ```julia 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(Ab) + 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(Ab) + 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 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
```
Often we want lazy realizations of matrices, which are supported via
ApplyArray. For example, the following creates a lazy matrix exponential:
julia 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.074A lazy matrix exponential is useful for, say, in-place matrix-exponetialvector: ```julia 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.
Lazy
vcatand
hcatallow 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
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) ```
We can represent Kronecker products of arrays without constructing the full array:
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.26381julia> C = Matrix{Float64}(undef, 6, 6); @btime copyto!(C, K); 61.528 ns (0 allocations: 0 bytes)
julia> C == kron(A,B) true
Base includes a lazy broadcast object called
Broadcasting, but this is not a subtype of
AbstractArray. Here we have
BroadcastArraywhich replicates the functionality of
Broadcastingwhile 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) ```