Profiling Howto

To choose an appropriate propagation method and parameters for a given problem, it is essential to benchmark and profile the propagation.

Consider the following simple example:

using QuantumPropagators: hamiltonian
using QuantumControlTestUtils.RandomObjects: random_dynamic_generator, random_state_vector

tlist = collect(range(0, step=1.0, length=101));
N = 200;  # size of Hilbert space
H = hamiltonian(random_dynamic_generator(N, tlist)...);
Ψ₀ = random_state_vector(N);

BenchmarkTools

The first line of defense is the use of BenchmarkTools. The @benchmark macro allows to generate statistics on how long a call to propagate takes.

Chebychev propagation

For example, we can time the propagation with the Chebychev method:

using BenchmarkTools
using QuantumPropagators
using QuantumPropagators: Cheby

@benchmark propagate($Ψ₀, $H, $tlist; method=Cheby, check=false) samples=10

BenchmarkTools.Trial: 10 samples with 1 evaluation per sample.
 Range (min … max):  20.378 ms …  20.925 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     20.630 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   20.610 ms ± 149.889 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▁        ▁    ▁   ▁        █ ▁  ▁  ▁                       ▁  
  █▁▁▁▁▁▁▁▁█▁▁▁▁█▁▁▁█▁▁▁▁▁▁▁▁█▁█▁▁█▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  20.4 ms         Histogram: frequency by time         20.9 ms <

 Memory estimate: 1023.66 KiB, allocs estimate: 3263.

Newton propagation

Or, the same propagation with the Newton method:

using QuantumPropagators: Newton

@benchmark propagate($Ψ₀, $H, $tlist; method=Newton, check=false) samples=10

BenchmarkTools.Trial: 10 samples with 1 evaluation per sample.
 Range (min … max):  80.171 ms … 84.078 ms  ┊ GC (min … max): 0.00% … 3.16%
 Time  (median):     82.892 ms              ┊ GC (median):    3.08%
 Time  (mean ± σ):   82.400 ms ±  1.441 ms  ┊ GC (mean ± σ):  2.22% ± 1.52%

  █         ▁                            █  ▁▁   ▁      ▁   ▁  
  █▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁██▁▁▁█▁▁▁▁▁▁█▁▁▁█ ▁
  80.2 ms         Histogram: frequency by time        84.1 ms <

 Memory estimate: 13.09 MiB, allocs estimate: 60295.

The result in this case illustrates the significant advantage of the Chebychev method for systems of moderate to small size and unitary dynamics.

When using custom data structures for the dynamical generators or states, @benchmark should also be used to optimize lower-level operations as much as possible, e.g. the application of the Hamiltonian to the state.

TimerOutputs

A lot more insight into the internals of a propagate call can be obtained by collecting timing data. This functionality is integrated in QuantumPropagators and uses the TimerOutputs package internally.

Enabling the collection of timing data

To enable collecting internal timing data, call QuantumPropagators.enable_timings:

QuantumPropagators.enable_timings()

The status of the data collection can be verified with QuantumPropagators.timings_enabled.

Chebychev propagation

Since the call to QuantumPropagators.enable_timings invalidates existing compiled code, and to avoid the compilation overhead showing up in the timing data, we call propagate once to ensure compilation:

propagate(Ψ₀, H, tlist; method=Cheby);

In any subsequent propagation, we could access the timing data in a callback to propagate:

function show_timing_data(propagator, args...)
    if propagator.t == tlist[end]
        show(propagator.timing_data, compact=true)
    end
end

propagate(Ψ₀, H, tlist; method=:cheby, callback=show_timing_data);

─────────────────────────────────────────────────────────────────────
                                        Time           Allocations
                                   ──────────────    ───────────────
        Tot / % measured:           829ms / 2.3%     22.5MiB / 0.5%
 ────────────────────────────────  ──────────────    ───────────────
 Section                   ncalls    time    %tot      alloc    %tot
─────────────────────────────────────────────────────────────────────
 prop_step!                   100  19.4ms  100.0%     127KiB  100.0%
 └─ matrix-vector product   1.30k  18.2ms   93.8%          ∅       ∅
─────────────────────────────────────────────────────────────────────

See the TimerOutputs documentation for details on how to print the timing_data.

Alternatively, without a callback:

propagator = init_prop(Ψ₀, H, tlist; method=Cheby)
for step ∈ 1:(length(tlist)-1)
    prop_step!(propagator)
end
show(propagator.timing_data, compact=true)

─────────────────────────────────────────────────────────────────────
                                        Time           Allocations
                                   ──────────────    ───────────────
        Tot / % measured:          40.7ms / 46.9%    476KiB / 26.6%
 ────────────────────────────────  ──────────────    ───────────────
 Section                   ncalls    time    %tot      alloc    %tot
─────────────────────────────────────────────────────────────────────
 prop_step!                   100  19.1ms  100.0%     127KiB  100.0%
 └─ matrix-vector product   1.30k  18.0ms   94.2%          ∅       ∅
─────────────────────────────────────────────────────────────────────

The reported runtimes here are less important than the number of function calls and the runtime percentages. In this case, the timing data shows that the propagation is dominated by the matrix-vector products (applying the Hamiltonian to the state), as it should. The percentage would go to 100% for larger Hilbert spaces.

Newton propagation

For the Newton method:

propagate(Ψ₀, H, tlist; method=Newton);   # recompilation
propagate(Ψ₀, H, tlist; method=Newton, callback=show_timing_data);

─────────────────────────────────────────────────────────────────────────────
                                                Time           Allocations
                                           ──────────────    ───────────────
            Tot / % measured:              157ms / 65.7%     14.9MiB / 87.3%
 ────────────────────────────────────────  ──────────────    ───────────────
 Section                           ncalls    time    %tot      alloc    %tot
─────────────────────────────────────────────────────────────────────────────
 prop_step!                           100   103ms  100.0%    13.0MiB  100.0%
 ├─ arnoldi!                          200  30.3ms   29.3%       224B    0.0%
 │  └─ matrix-vector product        2.00k  28.0ms   27.2%          ∅       ∅
 ├─ diagonalize_hessenberg_matrix     200  27.8ms   26.9%    9.38MiB   71.9%
 ├─ evaluate polynomial               200  21.7ms   21.0%    2.47MiB   19.0%
 ├─ get Leja points                   200  21.6ms   21.0%          ∅       ∅
 └─ get Newton coeffs                 200   207μs    0.2%          ∅       ∅
─────────────────────────────────────────────────────────────────────────────

We see here that the Newton propagation requires more matrix-vector products (2000 compared to 1200 for Chebychev), partly because the Newton propagator is "chunked" to m_max applications in each "restart" (10 by default, with 2 restarts required to reach machine precision in this case). Moreover, there is significant overhead beyond just matrix-vector multiplication, which will disappear only for significantly larger Hilbert spaces.

Disabling the collection of timing data

There there is a small overhead associated with collecting the timing data, it should not be enabled "in production". To QuantumPropagators.disable_timings function undoes the previous QuantumPropagators.enable_timings:

QuantumPropagators.disable_timings()

This again will trigger recompilation of any method that was collecting timing data, removing the associated overhead.