Matrix exponential

This page illustrates using the Julia package BlochSim to compute matrix exponentials, i.e., expm, that are needed for Bloch simulations.

This page comes from a single Julia file: expm-bench.jl.

You can access the source code for such Julia documentation using the 'Edit on GitHub' link in the top right. You can view the corresponding notebook in nbviewer here: expm-bench.ipynb, or open it in binder here: expm-bench.ipynb.

Setup

First we add the Julia packages that are need for this demo. Change false to true in the following code block if you are using any of the following packages for the first time.

if false
    import Pkg
    Pkg.add([
        "BenchmarkTools"
        "BlochSim"
        "ForwardDiff"
        "InteractiveUtils"
        "LaTeXStrings"
        "ExponentialAction"
        "LinearAlgebra"
        "MIRTjim"
        "Plots"
        "Random"
    ])
end

Tell this Julia session to use the following packages for this example. Run Pkg.add() in the preceding code block first, if needed.

using BenchmarkTools: @benchmark
using BlochSim: expm_bloch3, matrix_bloch3
using ExponentialAction: expv
using LinearAlgebra: I
using MIRTjim: prompt
using Plots: gui, plot, plot!, default
using Random: seed!

seed!(0)
default(titlefontsize = 10, markerstrokecolor = :auto, label="", width = 1.5)

The following line is helpful when running this file as a script; this way it will prompt user to hit a key after each figure is displayed.

isinteractive() || prompt(:draw);

Benchmark exp*

• ExponentialAction.expv: general-purpose matrix exponential

• BlochSim.expm_bloch3 specific matrix exponential for 3×3 Bloch matrix

r1 = rand() * 3
r2 = rand() * 0.1
ϕ = rand() * 2π
w = (rand() - 0.5) * 5
t = rand() * 0.4
Ω = (rand() - 0.5) * 7 # RF amplitude
s, c = sincos(ϕ) # for speed and consistency
s *= Ω
c *= Ω
A = matrix_bloch3(r1, r2, w, s, c)

3×3 Matrix{Float64}:
 -0.00685458  -2.07015      2.04616
  2.07015     -0.00685458  -1.73962
 -2.04616      1.73962     -1.2171

timing test

x = [r1, r2, w, s, c]
f3(x) = expm_bloch3(x..., t) # using BlochSim
b3 = @benchmark f3($x) # 1.3 μs (22 allocations: 1.20 KiB)

BenchmarkTools.Trial: 10000 samples with 10 evaluations per sample.
 Range (min … max):  1.491 μs …  6.709 ms  ┊ GC (min … max):  0.00% … 99.87%
 Time  (median):     1.536 μs              ┊ GC (median):     0.00%
 Time  (mean ± σ):   2.268 μs ± 67.074 μs  ┊ GC (mean ± σ):  29.54% ±  1.00%

  ▆█▇▆▅▄▃▁                                                   ▂
  █████████▇▇▆▅▆▆▅▆▆▆▅▆▅▆▄▄▃▃▁▃▁▃▁▁▁▁▁▁▁▃▁▁▁▃▁▃▁▁▁▃▄▇███▇▆▆▄ █
  1.49 μs      Histogram: log(frequency) by time     2.76 μs <

 Memory estimate: 1.42 KiB, allocs estimate: 26.

vs general-purse exp

fv(x) = expv(t, matrix_bloch3(x...), I(3)) # ExponentialAction
bv = @benchmark fv($x) # 3.8 μs (147 allocations: 10.4 KiB)

BenchmarkTools.Trial: 10000 samples with 8 evaluations per sample.
 Range (min … max):  3.683 μs … 726.380 μs  ┊ GC (min … max):  0.00% … 98.69%
 Time  (median):     3.951 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   4.920 μs ±  23.983 μs  ┊ GC (mean ± σ):  17.39% ±  3.55%

  ▁▅▇██▇▅▃▁                  ▁▁                               ▂
  █████████▇▇▇▇▇▅▄▅▆▅▅▄▅▄▅▇█████▇▆▆▆▆▄▅▅▄▄▁▄▅▄▄▄▄▅▆▅▅▆▆▅▅▅▅▆▆ █
  3.68 μs      Histogram: log(frequency) by time      7.44 μs <

 Memory estimate: 10.14 KiB, allocs estimate: 137.

it is more fair to precompute matrices:

I3 = I(3)
fp(t) = expv(t, A, I3) # ExponentialAction
bp = @benchmark fp($t) # 3.4 μs (128 allocations: 9.9 KiB)

BenchmarkTools.Trial: 10000 samples with 8 evaluations per sample.
 Range (min … max):  3.480 μs … 719.776 μs  ┊ GC (min … max):  0.00% … 98.68%
 Time  (median):     3.753 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   4.674 μs ±  23.691 μs  ┊ GC (mean ± σ):  18.08% ±  3.55%

     ▁▄▅▇█▄                                                    
  ▁▃▆███████▄▃▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  3.48 μs         Histogram: frequency by time        5.88 μs <

 Memory estimate: 9.86 KiB, allocs estimate: 128.

This page was generated using Literate.jl.