Operators : Squeeze, displaceop. Addition : expm, full

[amitjamadagni]

An attempt at implementing expm, full (sparse to dense) as expm is not implemented for sparse matrix. With respect to the implementation of expm with respect to sparse matrix, I have been referring to the JuliaLang/LinearAlgebra.jl#174 which had hints at scipy implementation. The scipy implementation is based on the following algorithm (http://core.ac.uk/download/pdf/17664.pdf) as cited in the documentation. I will try a hand at implementing this if already not implemented, but that said as exponentiation is only available for dense matrix, I have implemented full function and added the expm on the top of that.
Also the * method using DualMatrix was giving some inconsistent results, an attempt has been made to get it work.
Any comments would be helpful. Thanks

[acroy]

Wouldn't it be easier to just say expm(full(coeffs(qarr)))? But maybe it doesn't matter.

[amitjamadagni]

Yeah !!! Even I was thinking after I pushed. Let me check on the efficiency ...

[amitjamadagni]

Seems like expm(full(coeffs(qarr))) is more efficient.

[acroy]

I think for sparse matrices one mostly wants to avoid calculating expm directly, since the latter is typically dense and large. This is the reason why empv et al are so important. For those one rather computes expm*v where v is some vector. That said, I think that expm(full(S)) is fine for small problems. If one is really serious about it, one needs something like expmv for QuArrays and functions that displace/squeeze/etc a vector. If we want to include such functions in QuBase, we can pull in Expokit.jl or ExpmV.jl for instance.

[jrevels]

@acroy At some point, we should probably include a definition like the following to dispatch on for sparse QuArrays:

typealias SparseQuArray{B,T,N,S<:AbstractSparseArray} QuArray{B,T,N,S}

If we want to include such functions in QuBase, we can pull in Expokit.jl or ExpmV.jl for instance.

+1, and I think it makes sense to include those functions (in the form of using the packages, of course).

[acroy]

I would suggest to wait with using Expokit.jl et al. We have to look into it for the propagators anyways.

+1 to SparseQuArray.

[amitjamadagni]

One build passes and the other build fails. The error says "does not pass tests in QuBase" for the second build. Any insights would be helpful.

[acroy]

I guess the Travis failure is not related to this PR. Might be a consequence of JuliaLang/julia#10380. @jrevels : Any ideas?

[amitjamadagni]

@acroy done with the edits.

[acroy]

Looks Ok, but I want to wait for @jrevels to comment on the Travis failure. Could you add a test for displaceop? If one applies the displacement op to the ground state, one gets a coherent state and you can test if the elements in the state vector are correct.

[amitjamadagni]

@acroy I get an

ERROR: InexactError()
 in setindex! at sparse/sparsematrix.jl:1253
 in generic_scale! at linalg/generic.jl:18

when I use the combination of complex type for alpha parameter in displaceop. So I have replaced it with scale. I will work on this and see if this can be made to work if possible with scale.

Edit : And also the build is failing, when I run the tests on my machine it passes

sv1 = [0.5403023058681398,0.8414709848078965]
sv2 = [0.5403023058681398 + 0.0im,0.0 + 0.8414709848078965im]
println(displaceop(2,1)*statevec(1,FiniteBasis(2)))
println(QuArray(sv1))
println(@assert displaceop(2,1)*statevec(1,FiniteBasis(2)) == QuArray(sv1))
println(@assert displaceop(2,1*im)*statevec(1,FiniteBasis(2)) == QuArray(sv2))
println(displaceop(2,1) * statevec(1,FiniteBasis(2)) == QuArray(sv1))

# output 
2-element QuVector in FiniteBasis{Orthonormal,1}:
...coefficients: Array{Float64,1}
[0.5403023058681398,0.8414709848078965]
2-element QuVector in FiniteBasis{Orthonormal,1}:
...coefficients: Array{Float64,1}
[0.5403023058681398,0.8414709848078965]
nothing
nothing
true

which is as expected. But there is an error in the first build saying the assertion fails.

[acroy]

I guess you have to use @test_approx_eq or similar. Why do you "hard code" the expected result instead of calculating exp(-1/2) or whatever for the coefficients?

[amitjamadagni]

@acroy we could have done something like

exp(-0.5)*(statevec(1,FiniteBasis(2))+statevec(2,FiniteBasis(2)))
2-element QuVector in FiniteBasis{Orthonormal,1}:
...coefficients: Array{Float64,1}
[0.6065306597126334,0.6065306597126334]

which is not equal to the required output. (Truncated version)

We could define coherent state as

coherent(n::Int,alpha::Number) = displaceop(n,alpha)*statevec(1,FiniteBasis(n))

# then test as

coherent(2,1) == displaceop(2,1)*statevec(1,FiniteBasis(2))

We could also have the expansion of coherent states i.e., |alpha> = exp(-mod(alpha)^2)/2)*(sigma[(alpha)^n/sqrt(n!)|n>)

[acroy]

A function coherentstatevec should probably either be implemented using the displace function (which will use expmv) or directly using the expansion in Fock states. But your version is probably Ok for now.

Anyways, the test probably fails because your basis is too small. Using QuDOS.coherentstatevec I get

julia> cs = QuDOS.coherentstatevec(10,1.)
QuStateVec{Array{Complex{Float64},1}}(Complex{Float64}[0.606531+0.0im,0.606531+0.0im,0.428882+0.0im,0.247615+0.0im,0.123808+0.0im,0.0553686+0.0im,0.022603+0.0im,0.00854887+0.0im,0.00299672+0.0im,0.00110007+0.0im],10,[10])

julia> cs[1]
0.6065306597114603 + 0.0im

julia> exp(-0.5)
0.6065306597126334

julia> Base.Test.@test_approx_eq_eps real(cs[1]) exp(-0.5) 1e-6

You might be able to use a similar test if you increase the basis size.

[amitjamadagni]

@acroy I have made an attempt at implementing coherentstatevec using expm as well the expanding the summation. I have added the tests for displaceop using the same. A review would be helpful.

[acroy]

I don't think we need both versions and for now I would prefer the second one. There might actually be a problem having both, because they only differ in the keywords and I am not sure how Julia handles this. You can check by calling

@which coherentstatevec(20,1,method = "direct")
@which coherentstatevec(20,1)

To see which versions are actually called.

[amitjamadagni]

I have tried this out

println(@which coherentstatevec(10,1,method="direct"))
println(coherentstatevec(10,1,method="direct"))
println(@which coherentstatevec(10,1))
println(coherentstatevec(10,1))

# gives out 

coherentstatevec(n::Int64,alpha::Number) at QuBase.jl/src/arrays/constructors.jl:28
10-element QuVector in FiniteBasis{Orthonormal,1}:
...coefficients: Array{Float64,1}
[0.6065306597126334,0.6065306597126334,0.4288819424803534,0.24761510494160166,0.12380755247080083,0.05536842069051653,0.02260406309258736,0.008543532794657104,0.0030205949871958465,0.001006864995731949]
coherentstatevec(n::Int64,alpha::Number) at QuBase.jl/src/arrays/constructors.jl:28
10-element QuVector in FiniteBasis{Orthonormal,1}:
...coefficients: Array{Float64,1}
[0.6065306597114603,0.6065306597355303,0.4288819421804789,0.24761510797604577,0.12380752738214919,0.05536859516814174,0.02260302507996636,0.008548870306885925,0.0029967233467164565,0.0011000726308558629]

The output is different though they point to the same method which is reflected by the @which

[acroy]

To be sure you can put a println statement in one of the functions ...

[amitjamadagni]

I gave in a println for the rand(n) and this is the output ... seems like the method is being called but is not reflected by @which

coherentstatevec(n::Int64,alpha::Number) at QuBase.jl/src/arrays/constructors.jl:28

# printing rand(n)

[0.019244684430714143,0.3719091096873184,0.9047165798455694,0.6866397480495858,0.5561245580875909,0.14391135614954487,0.6180840544733457,0.002417775152962154,0.9427717780056364,0.20641868199486169]

# returning the output 

10-element QuVector in FiniteBasis{Orthonormal,1}:
...coefficients: Array{Float64,1}
[0.6065306597126334,0.6065306597126334,0.4288819424803534,0.24761510494160166,0.12380755247080083,0.05536842069051653,0.02260406309258736,0.008543532794657104,0.0030205949871958465,0.001006864995731949]


coherentstatevec(n::Int64,alpha::Number) at QuBase.jl/src/arrays/constructors.jl:28
10-element QuVector in FiniteBasis{Orthonormal,1}:
...coefficients: Array{Float64,1}
[0.6065306597114603,0.6065306597355303,0.4288819421804789,0.24761510797604577,0.12380752738214919,0.05536859516814174,0.02260302507996636,0.008548870306885925,0.0029967233467164565,0.0011000726308558629]

[acroy]

I solved the issue with the Travis failure (see JuliaLang/julia#26). Please rebase and squash after addressing the final comments.

[acroy]

We shouldn't export this function. It probably will cause trouble. Let's just use it internally and if there is demand we can export later.

[amitjamadagni]

But we are using it one of the tests. So I guess we need to export it.

[acroy]

You can call it via QuBase.coherentstatevec_inf. (In Julia you can always access all methods and types defined in a module by putting Modulename. in front. Exported methods are in a way just a convenience and in that sense define a "public API".)

[amitjamadagni]

@acroy
I have tried these things out. Locally they pass for me, but the build-bot shows an error

println(coeffs(displaceop(2,1.0)'))
println(coeffs(squeezingop(lowerop(2), QuArray(eye(2)),2.0)))
println(@assert coeffs(squeezingop(lowerop(2), QuArray(eye(2)), 2.0)') == coeffs(displaceop(2,1)))
println(@assert squeezingop(lowerop(2),QuArray(eye(2)), 2.0)' == displaceop(2,1)) 

[0.5403023058681398 0.8414709848078965
 -0.8414709848078965 0.5403023058681397]
[0.5403023058681398 0.8414709848078965
 -0.8414709848078965 0.5403023058681397]
nothing
nothing

[acroy]

Does it make a difference if you take coeffs(displaceop(2,1.0)) in the last assert (note the 1.0)? Otherwise you have to resort to @test_approx_eq or @test_approx_eq_eps.

[amitjamadagni]

But locally for me it is passing, I guess it should pass by the build-bot as well, right ? I am missing something here ??

[acroy]

It's strange but it can happen. In the worst case you have to put in println statements like in your local test.

[acroy]

The == in your @test_approx_eq commit is superfluous ...

[acroy]

Hah, nice! Could you please squash your commits. I will have another look tomorrow, but I think everything looks good now.