DSF for time-dependent Hamiltonian?

[leostenzel]

Bug Description

Hi,
I'm curious about trying the dynamic shifted Fock method for quantum annealing. But it doesn't seem to work with QobjEvos?
I'm not sure if I'm using it wrong, or how difficult it'd be to implement…

The error states something about an operator cache, I didn't dig any deeper:

ERROR: AssertionError: cache needs to be set up for operator of type
((MatrixOperator + ScalarOperator * MatrixOperator) ⊗ (MatrixOperator + ScalarOperator * MatrixOperator)). Set up cache by calling `cache_operator(L, u)`

The example below is just adapted from the QT.jl paper:

Code to Reproduce the Bug

using QuantumToolbox

tlist = range(0, 15, 1000)

f(p, t) = cos(t)

# Hamiltonian
function H(op_list, p)
    a1 = op_list[1]
    return a1' * a1 + a2' * a2 - (a2^2)' * a2^2 + a1 + a1' + a1' * a2^2 + QobjEvo(a1 * (a2^2)', f)
end

# Collapse operators
c_ops(op_list, p) = [op_list[1], a2, a2^2]

# Expectation operators
function e_ops(op_list, p)
    a1 = op_list[1]
    return [a1' * a1, a2' * a2, a1, a2]
end

N1 = 4
N2 = 40

a1 = tensor(destroy(N1), qeye(N2))
a2 = tensor(qeye(N1), destroy(N2))

ψ0 = kron(fock(N1, 0), fock(N2, 0))

mesolve(H([a1], ()), ψ0, tlist, c_ops([a1], ())) # Works

dsf_mesolve(H, ψ0, tlist, c_ops, [a1]) # Fails

Code Output

ERROR: AssertionError: cache needs to be set up for operator of type
((MatrixOperator(160 × 160) + ScalarOperator(0.0 + 0.0im) * MatrixOperator(160 × 160)) ⊗ (MatrixOperator(160 × 160) + ScalarOperator(0.0 + 0.0im) * MatrixOperator(160 × 160))). Set up cache by calling `cache_operator(L, u)`
Stacktrace:
  [1] mul!(w::SubArray{…}, L::SciMLOperators.TensorProductOperator{…}, v::SubArray{…})
    @ SciMLOperators ~/.julia/packages/SciMLOperators/gJ83E/src/tensor.jl:257
  [2] arnoldi_init!(A::SciMLOperators.TensorProductOperator{…}, b::Vector{…}, V::Matrix{…}, H::Matrix{…})
    @ QuantumToolbox ~/.julia/packages/QuantumToolbox/TVo9d/src/arnoldi.jl:25
  [3] arnoldi!(AS::ArnoldiSpace{…}, A::SciMLOperators.TensorProductOperator{…}, b::Vector{…})
    @ QuantumToolbox ~/.julia/packages/QuantumToolbox/TVo9d/src/arnoldi.jl:59
  [4] arnoldi(A::SciMLOperators.TensorProductOperator{Any, Tuple{…}, Nothing}, b::Vector{ComplexF64}, m::Int64)
    @ QuantumToolbox ~/.julia/packages/QuantumToolbox/TVo9d/src/arnoldi.jl:71
  [5] dsf_mesolveProblem(H::typeof(H_example), ψ0::QuantumObject{…}, tlist::StepRangeLen{…}, c_ops::typeof(c_ops), op_list::Vector{…}, α0_l::Vector{…}, dsf_params::@NamedTuple{}; e_ops::QuantumToolbox.var"#dsf_mesolve##0#dsf_mesolve##1", params::@NamedTuple{}, δα_list::Vector{…}, krylov_dim::Int64, kwargs::@Kwargs{})
    @ QuantumToolbox ~/.julia/packages/QuantumToolbox/TVo9d/src/time_evolution/time_evolution_dynamical.jl:364

Expected Behaviour

The call to dsf_mesolve should give similar results to the mesolve call.

Your Environment

Julia              Ver. 1.12.2
QuantumToolbox     Ver. 0.39.1
SciMLOperators     Ver. 1.13.0
LinearSolve        Ver. 3.48.0
OrdinaryDiffEqCore Ver. 1.36.0

System information:
====================================
OS       : macOS (arm64-apple-darwin24.0.0)
CPU      : 8 × Apple M1
Memory   : 8.0 GB
WORD_SIZE: 64
LIBM     : libopenlibm
LLVM     : libLLVM-18.1.7 (ORCJIT, apple-m1)
BLAS     : libopenblas64_.dylib (ilp64)

Additional Context

No response

[ytdHuang]

@albertomercurio,

This error also exists in v0.39.1, so I think there is some extra treatment in dsf_mesolve and requires extra operator cache. Could you have a look at this ?

I gave a look, but I don't understand the code structure of dsf_mesolve.