``potentialoperator()`` could use better documentation

[rafaelcoldatoms]

While working to reproduce the results from Fitzek (2020), I tried to implement a cosine- shaped potential in the following way (this is a minimal-working example):

using QuantumOptics

b_position = PositionBasis(-1, 1, 100)
b_momentum = MomentumBasis(b_position)

x = position(b_position)
p = momentum(b_momentum)

H_kin = p^2/2
V = cos(x)
H = H_kin + V

This return an error:

MethodError: no method matching cos(::Operator{PositionBasis{-1, 1, Int64, Int64}, PositionBasis{-1, 1, Int64, Int64}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}})
Closest candidates are:
  cos(::T) where T<:Union{Float32, Float64} at special/trig.jl:98
  cos(::LinearAlgebra.UniformScaling) at /usr/share/julia/stdlib/v1.8/LinearAlgebra/src/uniformscaling.jl:173
  cos(::DualNumbers.Dual) at ~/.julia/packages/DualNumbers/5knFX/src/dual.jl:311
  ...

Stacktrace:
 [1] top-level scope
   @ In[6]:13
 [2] eval
   @ ./boot.jl:368 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1428

It looks like there's no cosine method for operators. If this is truly a missing feature and not a bug, I think it should be fixed ASAP.

[Krastanov]

You have a couple of options for this, and there are a few reasons it might be good this is not automated.

First, you can always manually write a Taylor expansion of the cosine. That is not a terrible idea, in particular because it preserves some sparsity in the operator representation, even if it is only approximate.

Second, you can try to extract the underlying matrix and perform the cosine call yourself. In general, it is implemented in Julia by eigenvector decomposition. After all cos(rand(N,N)) works as expected in Julia. But see what happens if we try to work directly with the array object stored in x.data:

julia> cos(x.data)
ERROR: MethodError: no method matching eigen!(::LinearAlgebra.Hermitian{ComplexF64, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}; sortby::Nothing)

The error happens because a Hermitian sparse matrix is used to represent the underlying data. Eigendecomposition of sparse matrices generally does not preserve sparsity so many of these functions are not implemented at all (reasons explained below). So the developers are left with two options now:

1. Automatically convert the sparse matrix to dense matrix and then perform the eigendecomposition, etc. But then you will have users complaining about "My sparse matrices became dense without my permission and my RAM consumption exploded".
2. Expect users to do it manually so that they are aware of the significant expense they are incurring. This is generally the attitude of the core Julia language developers, for instance see this discussion: JuliaLang/julia#32899 (comment)

Now, back to this particular issue. If you are fine with losing sparsity, you can do this:

cosx_dense = Operator(x.basis_l, x.basis_r, cos(dense(x).data))

Some of this should indeed be implemented as part of QuantumOptics. While the conversion to dense matrix should certainly be explicit, some helpers for the rest of the operation would be nice. I am sure the devs here would be very grateful if you submit a pull request with it -- probably a whole new set of operator function methods.

Also, this object is actually still quite sparse. It would be wonderful to implement special methods for operator functions that preserve sparsity, but doing that the correct way is a lot of work. The linked issue above provides a list of many implementations of the "proper" methods to do it, but someone needs to do the work to bridge them to QuantumOptics.jl. If you are interested in starting this work, you will find a lot of support here.

If this is truly a missing feature and not a bug, I think it should be fixed ASAP.

I believe you did not meant it that way, but comments like this are not taken very well in the open source community. The folks that have developed QuantumOptics have mostly done it in their free time and such demanding tone from the users of free software can be very discouraging.

[Krastanov]

And for special cases like this, where sparsity (surprisingly) remains preserved, you can use this very roundabout way to make a sparse cos operator

sparse(Operator(x.basis_l, x.basis_r, cos(dense(x).data)))

[ChristophHotter]

Thank you @Krastanov, I also think it should not be automated.

@rafaelcoldatoms you can use the function potentialoperator.

using QuantumOptics

b_position = PositionBasis(-1, 1, 100)
b_momentum = MomentumBasis(b_position)

x = position(b_position)
p = momentum(b_position)

H_kin = p^2/2
V_cos(x) = cos(x)
V = potentialoperator(b_position, V_cos)
H = H_kin + V

Note that you also need to define p on the position basis for the sum.

[Krastanov]

Some of this should indeed be implemented as part of QuantumOptics

It seems it is actually already implemented :D

[rafaelcoldatoms]

Hey, thanks for the help!
First of all:

I believe you did not meant it that way, but comments like this are not taken very well in the open source community.
The folks that have developed QuantumOptics have mostly done it in their free time and such demanding tone from the users of free software can be very discouraging.

I'm truly sorry if I came out a inconsiderate. I am grateful for the community and the devs that keep up such wonderful projects. It was a keyboard slip-up, nothing more.

Now, regarding your comments:

It seems it is actually already implemented :D

It here still a need for a pull request? It does look like the suggestion by @ChristophHotter works as intended, and I wouldn't want to burden the devs with unnecessary work.

Thanks again!

[Krastanov]

No worries, I have occasionally misphrased things in a similar manner and just wanted to warn fellow community members about this common source of misunderstanding.

If the potentialoperator function does what you need, feel free to close this issue, or to rename it to something like "potentialoperator could use better documentation" (if you feel it is underdocumented), or you could submit a PR with changes to the documentation if you think it could be better.

[rafaelcoldatoms]

The potentialoperator function does look like it works for me. Thanks!

I looked in the documentation, and found a single use of potentialoperator in total, in the Wavepacket in 2D example example (other then the API itself of course). I do think that it can appear in more places, especially in the Schrödinger time evolution documentation.

So I'll change the name of the issue, as you suggested. Maybe later, when I'll have time to learn how to properly do a PR (never tried one before) I'll try to fix it myself if.