Prototype SciMLOperators-backed lazy operators and benchmark them

[Krastanov]

This issue is part of unitaryHACK26. You have to be registered to complete this issue.
Learn more about the PR submission process here and about unitaryHACK rules here!

A note about AI Slop: while we are open to collaboration with LLMs for unitaryHACK, fully AI-generated PRs are not acceptable. It is at the maintainers' discretion whether or not LLM-generated PRs are the right fit for the issues, and those that appear fully AI-generated may be immediately rejected. Read unitaryHACK's full AI Policy here.

Note from the maintainer: If you have not contributed before to open source projects, no AI agent use is permitted for this issue.

---

Context

QuantumOpticsBase.jl currently has custom lazy operator implementations:
LazySum, LazyProduct, and LazyTensor. These are important for avoiding
unnecessary materialization of full matrices, especially when applying structured
operators to states.

There is also a mature Julia package, SciMLOperators.jl, that provides a
matrix-free operator interface with lazy algebra, tensor-product operators,
cached mutating evaluation, and an AbstractMatrix-like API. This issue is to
prototype an alternative SciMLOperators-backed lazy operator path and benchmark
it against the current custom implementation.

Relevant code and docs:

• Current QuantumOpticsBase.jl lazy operators:
  • https://github.com/qojulia/QuantumOpticsBase.jl/blob/master/src/operators_lazysum.jl
  • https://github.com/qojulia/QuantumOpticsBase.jl/blob/master/src/operators_lazyproduct.jl
  • https://github.com/qojulia/QuantumOpticsBase.jl/blob/master/src/operators_lazytensor.jl
• Current lazy-operator documentation:
  • https://github.com/qojulia/QuantumOptics.jl/blob/master/docs/src/quantumobjects/operators.md#lazy-operators
• SciMLOperators.jl repository:
  • https://github.com/SciML/SciMLOperators.jl
• SciMLOperators.jl docs:
  • https://docs.sciml.ai/SciMLOperators/stable/
  • https://docs.sciml.ai/SciMLOperators/stable/interface/
  • https://docs.sciml.ai/SciMLOperators/stable/premade_operators/

Why this matters

The current custom lazy operators are tailored to QuantumOpticsBase.jl, but
they also duplicate functionality that exists in the wider SciML ecosystem.
SciMLOperators.jl supports matrix-free operators, lazy addition and
composition, tensor products, and cache_operator for allocation-reduced
in-place evaluation. A prototype would answer a concrete engineering question:

Can we reuse SciML's operator algebra without losing the performance and basis
semantics that QuantumOptics users rely on?

The expected result is not an immediate replacement. The expected result is a
tested prototype plus benchmark data that tells maintainers whether this path is
worth pursuing.

Minimal target example

The exact public API can be adjusted during implementation, but something like
the following should be possible:

using Test
using QuantumOpticsBase
using SciMLOperators

b0 = SpinBasis(1//2)
b = tensor(b0, b0, b0, b0)
sx = sigmax(b0)
sz = sigmaz(b0)

H_current = LazySum(
    LazyTensor(b, 1, sx),
    LazyTensor(b, 2, sz),
    LazyProduct(LazyTensor(b, 3, sx), LazyTensor(b, 4, sz)),
)

H_sciml = sciml_lazy_operator(H_current) # name is flexible

psi = Ket(b, randn(ComplexF64, length(b)))

@test dense(H_sciml) ≈ dense(H_current)
@test H_sciml * psi ≈ H_current * psi

Scope

In scope:

• Add an experimental SciMLOperators-backed operator representation for
  QuantumOpticsBase.jl operators.
• Preserve QuantumOpticsBase basis information (basis_l, basis_r) around
  the wrapped SciML operator.
• Provide conversion from common QuantumOpticsBase operators to SciML
  operators, for example dense/sparse operators to MatrixOperator.
• Support equivalents of the current lazy building blocks:
  • lazy sums via SciMLOperators lazy addition,
  • lazy products/composition via SciMLOperators lazy composition,
  • lazy tensor products via TensorProductOperator, including identity factors
    for missing subsystems in LazyTensor.
• Implement enough QuantumOpticsBase methods for correctness testing and
  benchmarking: at minimum dense, mul!/operator application to Ket, basic
  scalar multiplication, addition, multiplication, and dagger where feasible.
• Add a benchmark suite under QuantumOpticsBase.jl/benchmark/ using
  BenchmarkTools.jl.
• Benchmark the current custom lazy implementation against the
  SciMLOperators-backed prototype, including cached and uncached SciML
  evaluation where relevant.

Out of scope:

• Replacing LazySum, LazyProduct, or LazyTensor as the default
  implementation.
• Changing public APIs in QuantumOpticsBase.jl unless the maintainer agrees
  they are needed for the prototype.
• GPU support.
• Time-dependent operators, unless they are trivial to support through
  SciMLOperators' update interface.
• Full integration with QuantumOptics.jl time evolution solvers beyond a small
  benchmark smoke test.

Suggested implementation notes

• Prefer an optional dependency or extension-based design if practical, so
  SciMLOperators.jl does not become an unconditional dependency of
  QuantumOpticsBase.jl unless maintainers decide that is acceptable.
• Keep the prototype isolated, for example in an extension or a small
  SciMLOperators compatibility module.
• Avoid optimizing by materializing the full operator before applying it to a
  state. The point of the benchmark is matrix-free or lazy application.
• Compare against the existing lazy types as they are normally used, not against
  manually simplified dense matrices.
• Use cache_operator for a fair SciMLOperators in-place benchmark, since that
  is the documented path for reducing allocations.

Benchmark requirements

Add repeatable benchmarks covering at least these cases:

• LazySum: sum of local single-site terms on spin chains with increasing
  subsystem count, applied to a random Ket.
• LazyProduct: product/composition of several compatible operators, applied
  to a random Ket.
• LazyTensor: local operators embedded into a composite basis, including cases
  with missing identity subsystems.
• Mixed structured Hamiltonian: a small transverse-field or Heisenberg-like
  Hamiltonian built as a sum of tensor-product terms.
• Optional benchmark: one short time-evolution call using the existing
  custom lazy operator and the SciML-backed prototype, if the prototype supports
  the required mul! path.

Acceptance criteria

• A SciMLOperators-backed prototype exists and can be constructed from at least
  LazySum, LazyProduct, and LazyTensor examples.
• Correctness tests show that dense conversion and application to kets agree
  with the current custom implementation on representative examples.
• The benchmark suite is checked in under QuantumOpticsBase.jl/benchmark/.
• Benchmark results are summarized in the PR description, including whether
  SciMLOperators is faster, slower, or allocation-heavier in each tested case.

Difficulty

Advanced. This requires understanding the QuantumOpticsBase.jl operator and
basis interfaces, the current custom lazy implementations, and the
SciMLOperators.jl matrix-free operator model.

[swiatlowemnie333]

hi, i would like to take a look at this issue for unitary-hack if it is still free. setting up my local envirnoment right now to check how to map it to SciMLOperators
Zazul

[Krastanov]

Issues do not get reserved, you can work on it without asking :)

[Adithyaphani]

@Krastanov PR is up at qojulia/QuantumOpticsBase.jl#233 — linking here for tracking.

Two implementation findings worth noting:

QuantumOpticsBase uses column-major ordering (site 1 fastest-varying), while SciMLOperators' TensorProductOperator uses row-major (first argument slowest-varying). These conventions are opposite, so the site operator order needs to be reversed during conversion — otherwise the operator constructs fine but applies to the wrong tensor factor.
Any operator tree containing a TensorProductOperator requires cache_operator to be called before mul!, unlike MatrixOperator.

Both are documented in the PR. The benchmark covers edge/mid sites, dense/sparse factors, spin-1/2 and spin-1 local dimensions, and Heisenberg-like mixed Hamiltonians as requested.

[cmdr-chara]

@Krastanov I also have an implementation open at qojulia/QuantumOpticsBase.jl#230.

The current branch is back to one focused commit, with Buildkite and Codecov green. It keeps the SciMLOperators path isolated as an optional extension, preserves the QuantumOpticsBase basis wrapper, and includes conversions for LazySum, LazyProduct, and LazyTensor.

The benchmark suite compares the current lazy path against uncached and cached SciMLOperators-backed application. It covers sparse and dense tensor factors, spin-1/2 and spin-1 local dimensions, edge and non-edge tensor sites, and mixed Hamiltonian cases.

I noticed qojulia/QuantumOpticsBase.jl#233 is also open for this issue. If there is a specific benchmark or implementation detail you want compared between the two approaches, I can add that to qojulia/QuantumOpticsBase.jl#230.

[Adithyaphani]

@Krastanov PR : qojulia/QuantumOpticsBase.jl#233 , looking forward to your review and response on my PR, happy to make further changes if needed.