FFTOperators on sparse operators gives wrong answers

[AmitRotem]

While going for a FFT style momentum operation in the position basis i noticed the FFT fails on sparse matrices;

x_basis = PositionBasis(-1,1,10);
p_basis = MomentumBasis(x_basis);

# momentum operator in position basis
p_in_x_basis = transform(x_basis, p_basis) * momentum(p_basis) * transform(p_basis, x_basis);
p_in_x_basis ≈ momentum(x_basis) # false

The lhs is a diagonal matrix
Diagonal(p_in_x_basis.data) == (p_in_x_basis).data # true

The definition in the code is similar, but with a dense operator
transform(x_basis, p_basis) * dense(momentum(p_basis)) * transform(p_basis, x_basis) == momentum(x_basis) # true

My motivation was to get LazyProduct( transform(x_basis, p_basis), momentum(p_basis)^2, transform(p_basis, x_basis) )

Working with QuantumOpticsBase v0.5.4

[david-pl]

Hm, that is strange.

For what it's worth, your resolution of 10 in that range is quite low. So I'd never expect the ≈ operator to return true here.

That said, increasing the step-size doesn't work. I think there's really a bug in the multiplication of a transform operator with a sparse operator. The resulting operator is always still diagonal, which shouldn't be the case IMO. I'll have to investigate further. The reason no one noticed this so far is that you usually don't apply it to sparse operators, but rather states which are either Kets or dense operators.

Still, you should be fine with your application of LazyProduct( transform(x_basis, p_basis), momentum(p_basis)^2, transform(p_basis, x_basis) ) as that doesn't apply the transform to the momentum operator. Rather, it applies it to whatever you then multiply the LazyProduct instance with and then multiplies the momentum operator with the result. So long as that's a Ket or a (dense) operator you're fine.

[david-pl]

Okay, I found the reason. The problem is due to the in-place application of the planned fft!s. For example, in this line here

QuantumOpticsBase.jl/src/particle.jl

Line 591 in 12765d5

A.fft_r2! * reshape(data, n...,n...)

While the reshape is fine as it uses shared data, the application of a AbstractFFT.Plan to a sparse array falls back to the generic method here:

https://github.com/JuliaMath/AbstractFFTs.jl/blob/d64a878f3bda8b883a33f7bb05518ee8d8dfc707/src/definitions.jl#L224

And that appears to create a copy of the data, so we are applying the fft! to a copy that no longer shares the data.

Just as a test, I tried modifying the fft! application explicitly copying the result of the fft! back into the data:

data[:] .= (A.fft_r2! * reshape(data, n...,n...))[:]

With this, your test passes. However, if we just update the method like this, it might introduce performance regressions for dense arrays. We could potentially dispatch on dense arrays here and use the above for AbstractArray cases, as the fft! application copies stuff there.

[AmitRotem]

What are your thoughts about making the Lazy product the default?
Applying 2 dft and a diagonal matrix (either sparse or dense) is more efficient than applying the full matrix representation.

[david-pl]

While it doesn't sound completely unreasonable, I'd still vote against it. At no other point in the package does * return a LazyProduct. I'd prefer to keep semantics consistent so users don't get surprised here.

On a more general note: it never really makes sense to use a sparse array as output here as the result is always dense. So another approach would be to return a dense matrix here and keep the existing implementation. If we do that we have to make sure we don't accidentally convert other special array types (e.g. CuArray) to standard Matrix types.