Implement iterative steady-state solver

[david-pl]

This is my attempt to revive #252.

@PhilipVinc @Z-Denis I'd appreciate your input/review on this. Large parts of this are just copy-pasted code from https://github.com/Z-Denis/SteadyState.jl, I hope that's okay.

Some other comments:

• This supports lazy operators as well as ones with a matrix representation.
• I always cache Jdagger and Hnhdagger, but since QO.jl v0.8 adjoints are lazy (i.e. dagger(op) just creates an operator with an <:Adjoint as data) so it shouldn't matter.
• GPU support should also be there, but I didn't test it.

Are there any other features of SteadyState.jl that are held back with this implementation?

[codecov[bot]]

Codecov Report

Merging #308 (5f2bd3f) into master (2751c4c) will decrease coverage by 0.00%. The diff coverage is 98.27%.

@@            Coverage Diff             @@
##           master     #308      +/-   ##
==========================================
- Coverage   97.98%   97.98%   -0.01%     
==========================================
  Files          15       16       +1     
  Lines        1288     1337      +49     
==========================================
+ Hits         1262     1310      +48     
- Misses         26       27       +1     

Impacted Files	Coverage Δ
src/steadystate.jl	91.17% <ø> (ø)
src/steadystate_iterative.jl	97.61% <97.61%> (ø)
src/master.jl	97.57% <100.00%> (+0.10%)	:arrow_up:

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 2751c4c...5f2bd3f. Read the comment docs.

[PhilipVinc]

Hi @david-pl , Thanks for doing this! I think it will be a great addition!

Im working mostly in Python nowadays (unfortunately) but I’ll be glad to review this.

Im busy those days so give me a few days to get it.

[david-pl]

@PhilipVinc Thanks for the quick response. Sure thing, take your time. I appreciate the effort!

[Z-Denis]

Hi @david-pl,

Thank you very much for taking the time to add this. I am a bit busy too, but I would be more than glad to review it.

[davidschlegel]

I am happy to see that iterative solvers are being added to the package. I am using my own implementation of bicgstab, which is also using Reverse-Cuthill-McKee reordering, which can make the problem much less hard to solve. For questions, feel free to ask.

[PhilipVinc]

Overall it's very nice! The recent changes to QuantumOptics make everything easier, indeed...

The only thing that is missing is proper attribution. Let me check my notes. While the method is trivial, the first person who I know talked about this trick is Prof V. Savona (EPFL), so it would be correct to put a brief comment about it in the source.

[PhilipVinc]

@davidschlegel Can you explain what is the Reverse-Cuthill-McKee reordering you are talking about and what is your algorithm? Maybe it's not relevant to this PR, but i'm still interested.

This PR is about solving the system of equations L\rho = 0 and Tr[\rho] = 1.

[davidschlegel]

@PhilipVinc , yes correct, in the function definition of f(x,y), this is done correctly. What I am talking about is that in principle, you could scale the trace condition on both sides of the equation, which is mathematically equivalent, but not numerically. So rather than taking something like

y[end] = tr(rho)
x[end] = 1

one could for instance take

weight = mean(abs.(L.data.nzval))
y[end] = weight*tr(rho)
x[end] = weight

Because you don't compute L directly, but by a linear operation, I don't know how this can be achieved here. The only point of this is that the trace condition should be enforced on the same scale as the solution of rho.

[davidschlegel]

Regarding Reverse-Cuthill-McKee reordering I simply use the following with the help of SymRCM

using SymRCM

    #Reordering matrix
    # reverse Cuthill-McKee ordering to make matrix diagonally dominant
    p = symrcm(L.data)
    ip = similar(p) # inverse permutation
    ip[p] = 1:length(p)
    sol = iterativemethod(L.data[p, p], x[p],...)[ip]

Again, because you don't explicitly construct the Liouvillian, I am not sure whether or how this can be done. Another thing that could also improve things which is well known for bicgstabl is preconditioning of the problem.

[PhilipVinc]

Ah, I see! This is interesting. As L=H\otimes 1 + 1 otimes H + L otimes L^dag, even if it is lazy-computed, I suppose that

mean(abs.(L.data.nzval))

is going to be upper bounded by something like

mean(abs.(H.data.nzval)) + mean(abs.(L.data.nzval).^2))

so one could try to do something like that. Though i'd leave this for a follow up PR and as an optional keyword.

We used extensively the current algorithm and it works very well.

As for the Ordering, that is interesting. In general not building the full liouvillian gives you a huge speedup so it would be interesting to see whever it is possible to precondition the linear operator.

It should be possible to do it, no?

[davidschlegel]

Yes, I think this would actually work! I think, typically it works fine by just setting the unit trace condition as it is. I think there is not much research on how it can be improved by setting the weight.

Regarding preconditioning, I don't know. Something which is known to give a substantial speedup is incomplete LU decomposition. I am not sure how this can be done implicitly on a linear operator approach in a smaller Hilbert space.

[PhilipVinc]

How does this particular preconditioner work?

[davidschlegel]

Interesting, if bicstabl accepts linear operators (defined by their action) as a precondition, then one could for example use something from LinearOperators, LLDL, LimitedLDLFactorizations or MUMPS.

[david-pl]

Seems good to me.

@PhilipVinc Cool!

The only thing that is missing is proper attribution. Let me check my notes.

What about this though?

Sorry for the delay, I was trying to test it on a GPU.

@Z-Denis Thanks! Unfortunately, I have the same issue (the only GPU I have available is quite old). I noticed a mistake I made though. This x0 = zeros(eltype(rho0), M^2+1) doesn't conserve container type. It should be x0 = similar(rho0.data, M^2+1) so everything really runs on the GPU when working with CuArrays. I'll fix it.

[davidschlegel]

Alright, I implemented a method using RCM-reordering and iLU preconditioning. It uses the full Liouvillian however to do so. It can be nicely integrated into the existing method via multiple dispatch. Should I however wait until this gets into the master?

[david-pl]

@davidschlegel Sorry for the delay. I'll merge this today so you can open a new PR for your method.