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.
Vcat is the lazy analogue of vcat. For lazy vectors like ranges, it
creating such a vector is allocation-free. copyto! allows for allocation-free
population of a vector.
julia> using LazyArrays, BenchmarkTools
julia> A = Vcat(1:5,2:3) # allocation-free
7-element Vcat{Int64,1,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)Similarly, Hcat is the lazy analogue of hcat.
julia> A = Hcat(1:3, randn(3,10))
3×11 Hcat{Float64,Tuple{UnitRange{Int64},Array{Float64,2}}}:
1.0 0.350927 0.339103 -1.03526 … 0.786593 0.0416694
2.0 -1.10206 1.52817 0.223099 0.851804 0.430933
3.0 -1.26467 -0.743712 -0.828781 -0.0637502 -0.066743
julia> Matrix(A) == hcat(A.arrays...)
true
julia> b = Array{Int}(undef, length(A)); @btime copyto!(b, A);
26.670 ns (0 allocations: 0 bytes)
julia> B = Array{Float64}(undef, size(A)...); @btime copyto!(B, A);
109.625 ns (1 allocation: 32 bytes)
julia> @btime hcat(A.arrays...); # takes twice as long due to memory creation
274.620 ns (6 allocations: 560 bytes)Base now 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> A = randn(6,6);
julia> B = BroadcastArray(exp, A);
julia> Matrix(B) == exp.(A)
true
julia> B = BroadcastArray(+, A, 2);
julia> B == A .+ 2
trueFollowing Base's lazy broadcasting, we introduce lazy multiplication. The type
Mul support multiplication of any two objects, not necessarily arrays.
(In the future there will be a MulArray a la BroadcastArray.)
Mul is designed to work along with broadcasting, and to lower to BLAS calls
whenever possible:
julia> A = randn(5,5); b = randn(5); c = randn(5); d = similar(c);
julia> d .= 2.0 .* Mul(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 .* Mul(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)We also have lazy inverses PInv(A), designed to work alongside Mul to
to lower to BLAS calls whenever possible:
julia> A = randn(5,5); b = randn(5); c = similar(b);
julia> c .= Mul(PInv(A), b)
5-element Array{Float64,1}:
-2.5366335879717514
-5.305097174484744
-9.818431932350942
2.421562605495651
0.26792916096572983
julia> c .= Ldiv(A, b) # shorthand for above
5-element Array{Float64,1}:
-2.5366335879717514
-5.305097174484744
-9.818431932350942
2.421562605495651
0.26792916096572983We can represent Kronecker products of arrays without constructing the full array.
julia> A = randn(2,2); B = randn(3,3);
julia> K = Kron(A,B)
6×6 Kron{Float64,2,Tuple{Array{Float64,2},Array{Float64,2}}}:
1.99255 -1.45132 0.864789 -0.785538 0.572163 -0.340932
-2.7016 0.360785 -1.78671 1.06507 -0.142235 0.70439
1.89938 -2.69996 0.200992 -0.748806 1.06443 -0.0792386
-1.84225 1.34184 -0.799557 -2.45355 1.7871 -1.06487
2.49782 -0.333571 1.65194 3.32665 -0.444258 2.20009
-1.75611 2.4963 -0.185831 -2.33883 3.32464 -0.247494
julia> C = Matrix{Float64}(undef, 6, 6); @btime copyto!(C, K);
61.528 ns (0 allocations: 0 bytes)
julia> C == kron(A,B)
true
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent