Faster prephasors with linear programming

[paquiteau]

The problem

In #276 and #276 we added support for designing gradient waveform to connect two point in k-spac effiently, or so I though. While working on new trajectories, I found some MATLAB implementation use to stitches spokes together, and that was using Linear Programming.

TLDR: I decided to try it and compare with the current implementation:

• Blue is the trajectory to connect, red is using Linear Programming, green is current method using get_gradient_amplitudes_to_travel_for_set_time
• The Linear Programming solution finds shorter path (22 points vs >500) between spokes
• Computing the connecting path is slower (but can be parallellize between spokes)
• Unlike the existing implementation, this does not enforce the use of trapezes for the gradient waveforms.

A prototype and more insights can be found in: lip.ipynb

Longer explanations

Math Formulation of Linear Programming

• Variables: The gradient values at each time point for the x, y, and z axes.
• LP Constraints:
  • Integral: The sum of gradient values must equal a target k-space displacement ($\Delta k$).
  • Slew Rate: The change between consecutive points must not exceed a maximum rate ($s_{max}$).
  • Amplitude: The values must be within a maximum amplitude ($g_{max}$).
  • Boundary Conditions: The start and end values are fixed.

All constraints are linear, allowing the use of a standard LP solver. The constraints are formalized as matrices and vectors: see the Scipy doc for how this is done.

However, A big drawback of Linear Programming is that the length of the solution vector must be known, and if the prescribe size is to small, there is no solution available.

To circumvent that, I implemented a binary search between a lower and upper bound on the number of points required to create the connection, and the goal is to find the smallest solution that is acceptable.

More Gotcha's:

• The method is slow (4s for 256 spokes, vs 0.18 with the existing solution), but can be parallelize (e.g. using joblib) to perform multiple connections.
• differents connections path may have different length (the current implementation enforce the same lenght for all connection to compute), which could create some problems for echo spacing for a trajectory that passes multiple time in the center during a single (connected) shot.
• Getting the number of points right is tricky, if it fails outside the currently provided bounds, it fails hard. Getting tighter bounds also help the binary search to go faster.
• I think some edge cases (cartesian displacement for instance) could be done indepently, and faster (instead of solving of gx,gy,gz, we only have to solve the simpler problem for gx instead, the other waveform being zeroed.)
• LLM were use in the design of this, and I am not an expert in Linear Programming.

Next steps

• I would like some feedback on this methodology, and help on checking that it does work as we expect.
• Open a PR with
• Write some more lengthy documentation (that could take the form of a blog post)

References

• Linear Programming: https://en.m.wikipedia.org/wiki/Linear_programming
• Usage in MRI: I think I have seen a paper talking about it, but can't find it again....

[paquiteau]

Alright, after some tinkering, I made it twice as fast as before (reducing the problem to a 1D case, and only optimizing the length for the dimension with the largest area to cover), and clean the code a bit. I also added joblib support for multi spokes computations. Looks fast enough for me :)
lip.ipynb

[chaithyagr]

The optimizer is clearly underperforming. In my indivitual case by case solutions, I got very fast changes. Ill get back to this soon, but its lower priority for now for me.