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.233 μs …  10.531 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.293 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.319 μs ± 226.596 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ██▇▇▃▁                                                      
  ▆██████▆▅▄▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▂▁▂▂▂▂▂▂ ▃
  1.23 μs         Histogram: frequency by time        2.15 μ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.012 μs … 594.688 μs  ┊ GC (min … max):  0.00% … 98.72%
 Time  (median):     3.255 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   4.074 μs ±  20.722 μs  ┊ GC (mean ± σ):  18.86% ±  3.69%

       ▃▅█▇▆▄▁                                                 
  ▂▂▃▅▇███████▆▄▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▁▁▂▂▂▂▂▁▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
  3.01 μs         Histogram: frequency by time        4.73 μ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 9 evaluations per sample.
 Range (min … max):  2.848 μs … 518.707 μs  ┊ GC (min … max):  0.00% … 98.70%
 Time  (median):     3.063 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   3.832 μs ±  18.646 μs  ┊ GC (mean ± σ):  18.68% ±  3.81%

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

 Memory estimate: 9.86 KiB, allocs estimate: 128.

This page was generated using Literate.jl.