Speed up calculation in Physics simulation in Matlab

0
votes

I am working on a MR-physic simulation written in Matlab which simulates bloch's equations on an defined object. The magnetisation in the object is updated every time-step with the following functions. 

function Mt = evolveMtrans(gamma, delta_B, G, T2, Mt0, delta_t)
    % this function calculates precession and relaxation of the
    % transversal component, Mt, of M
    delta_phi = gamma*(delta_B + G)*delta_t;        
    Mt = Mt0 .* exp(-delta_t*1./T2 - 1i*delta_phi);
end

This function is a very small part of the entire code but is called upon up to 250.000 times and thus slows down the code and the performance of the entire simulation. I have thought about how I can speed up the calculation but haven't come up with a good solution. There is one line that is VERY time consuming and stands for approximately 50% - 60% of the overall simulation time. This is the line,

Mt = Mt0 .* exp(-delta_t*1./T2 - 1i*delta_phi);

where

Mt0 = 512x512 matrix

delta_t = a scalar

T2 = 512x512 matrix

delta_phi = 512x512 matrix

I would be very grateful for any suggestion to speed up this calculation.

More info below,

The function evovleMtrans is called every timestep during the simulation. The parameters that are used for calling the function are,

gamma = a constant. (gyramagnetic constant)

delta_B = the magnetic field value

G = gradientstrength

T2 = a 512x512 matrix with T2-values for the object

Mstart.r = a 512x512 matrix with the values M.r had the last timestep

delta_t = a scalar with the difference in time since the last calculated M.r

The only parameters of these that changed during the simulation are, G, Mstart.r and delta_t. The rest do not change their values during the simulation.

The part below is the part in the main code that calls the function. 

        % update phase and relaxation to calcTime

        delta_t = calcTime - Mstart_t;
        delta_B = (d-d0)*B0;
        G = Sq.Gx*Sq.xGxref + Sq.Gz*Sq.zGzref;

        % Precession around B0 (z-axis) and B1 (+-x-axis or +-y-axis)
        % is defined clock-wise in a right hand system x, y, z and
        % x', y', z (see the Bloch equation, Bloch 1946 and Levitt
        % 1997). The x-axis has angle zero and the y-axis has angle 90.
        % For flipping/precession around B1 in the xy-plane, z-axis has
        % angle zero.
        % For testing of precession direction:
        % delta_phi = gamma*((ones(size(d)))*1e-6*B0)*delta_t;
        M.r = evolveMtrans(gamma, delta_B, G, T2, Mstart.r, delta_t);
        M.l = evolveMlong(T1, M0.l, Mstart.l, delta_t);
2
matlabperformancesimulation
Are some of your variables like for example the relaxation time T2 or gamma actually constant during the evolution? Could you eloborate on why you are calling the function so often? If you are calling it for a very fine grid of delta_t it might be better to use a differential equation solver than using the explicit solution - Jannick
All T2 values and the gamma variable are constants. When I run the simulation I choose how many sample points I want i k-space. So when I for example choose 512x512 sample points I have to update the magnetisation at least 512*512 times with a small delta_t in between. That is why I call the function so many times. - Jens
I think it would be helpful if you could post some more code/information. You are aware that evolveMtrans includes the full analytical solution and is true for arbitrarily large delta_t? - Jannick
Thanks! I will post some more code/ information above. Yes, I am aware that it is true for larger delta_t. I use small delta_t to be able to sample the signal with small timesteps. - Jens

2 Answers

2
votes

This is not a surprise.

That "single line" is a matrix equation. It's really 1,024 simultaneous equations.

Per Jannick, that first term means element-wise division, so "delta_t/T[i,j]". Multiplying a matrix by a scalar is O(N^2). Matrix addition is O(N^2). Evaluating exponential of a matrix will be O(N^2).

I'm not sure if I saw a complex argument in there as well. Does that mean complex matricies with real and imaginary entries? Does your equation simplify to real and imaginary parts? That means twice the number of computations.

Your best hope is to exploit symmetry as much as possible. If all your matricies are symmetric, you cut your calculations roughly in half.

Use parallelization if you can.

Algorithm choice can make a big difference, too. If you're using explicit Euler integration, you may have time step limitations due to stability concerns. Is that why you have 250,000 steps? Maybe a larger time step is possible with a more stable integration schema. Think about a higher order adaptive scheme with error correction, like 5th order Runge Kutta.

0
votes

There are several possibilities to improve the speed of the code but all that I see come with a caveat.

Numerical ode integration

The first possibility would be to change your analytical solution by numerical differential equation solver. This has several advantages

  1. The analytical solution includes the complex exponential function, which is costly to calculate, while the differential equation contains only multiplication and addition. (d/dt u = -a u => u=exp(-at))
  2. There are plenty of built-in solvers for matlab available and they are typically pretty fast (e.g. ode45). The built-ins however all use a variable step size. This improves speed and accuracy but would be a problem if you really need a fixed equally spaced grid of time points. Here are unofficial fixed step solvers.

As a start you could also try to use just an euler step by replacing

M.r = evolveMtrans(gamma, delta_B, G, T2, Mstart.r, delta_t);

by

delta_phi = gamma*(delta_B + G)*t_step;
M.r += M.r .* (1-t_step*1./T2 - 1i*delta_phi);

You can then further improve that by precalculating all constant values, e.g. one_over_T1=1/T1, moving delta_phi out of the loop.

Caveat: You are bound to a minimum step size or the accuracy suffers. Therefore this is only a good idea if you time-spacing is quite fine.

Less points in time

You should carfully analyze whether you really need so many points in time. It seems somewhat puzzling to me that you need so many points. As you know the full analytical solution you can freely choose how to sample the time and maybe use this to your advantage.

Going fortran

This might seem like a grand step but in my experience basic (simple loops, matrix operations etc.) matlab code can be relatively easily translated to fortran line-by-line. This would be especially helpful in addition to my first point. If you still want to use the full analytical solution probably there is not much to gain here because exp is already pretty fast in matlab.