JuliaAttic / QuBase.jl

A foundational library for quantum mechanics in Julia
Other
43 stars 6 forks source link

Sparse support for QuArray #50

Open amitjamadagni opened 7 years ago

amitjamadagni commented 7 years ago

Currently the functionality is not supported and it would be great if we could have a work around as it would improve the efficiency in few cases. In sense for a many body system the sparse support would increase the efficiency and performance when used in conjunction to the solvers in QuDynamics. There have been some discussions in the past but any pointers on this would be great. Consider the following :

julia> v = statevec(1, FiniteBasis(2))
2-element QuVector in QuBase.FiniteBasis{QuBase.Orthonormal}:
...coefficients: Array{Float64,1}
[1.0,0.0]

julia> @time tensor(v,v)
  0.112194 seconds (104.74 k allocations: 4.530 MB)
4-element QuVector in QuBase.FiniteBasis{QuBase.Orthonormal}:
...coefficients: Array{Float64,1}
[1.0,0.0,0.0,0.0]

julia> v_s = SparseVector([1;0])
Sparse vector of length 2 with 1 Int64 nonzero entries:
  [1]  =  1

julia> function Base.kron(v1::SparseVector, v2::SparseVector)
           return kron(SparseMatrixCSC(v1), SparseMatrixCSC(v2))
       end

julia> @time kron(v_s, v_s)
  0.082897 seconds (42.55 k allocations: 1.709 MB)
4×1 sparse matrix with 1 Int64 nonzero entries:
    [1, 1]  =  1

julia> @time tensor(v,v)
  0.000153 seconds (26 allocations: 960 bytes)
4-element QuVector in QuBase.FiniteBasis{QuBase.Orthonormal}:
...coefficients: Array{Float64,1}
[1.0,0.0,0.0,0.0]

julia> @time kron(v_s, v_s)
  0.000010 seconds (20 allocations: 1.203 KB)
4×1 sparse matrix with 1 Int64 nonzero entries:
    [1, 1]  =  1

Though kron for SparseVectors is still not in base, the above has been constructed using the following issues : https://github.com/JuliaLang/julia/issues/19426 https://github.com/JuliaInv/jInv.jl/pull/26/files

It would be great if we could make QuArray entirely return a sparse construction instead of a dense version, or create an option if the coeffs construct should be sparse or dense with the default being set to a sparse construct. Something like :

julia> statevec(1, FiniteBasis(2))
2-element QuVector in QuBase.FiniteBasis{QuBase.Orthonormal}:
...coefficients: SparseVector
  [1]  =  1

julia> statevec(1, FiniteBasis(2), dense=true)
2-element QuVector in QuBase.FiniteBasis{QuBase.Orthonormal}:
...coefficients: Array{Float64,1}
[1.0,0.0]

One other option would be to just add sparse functionality on top of the current functionality, which I am not sure would be a good option. Any pointers on this would be useful.