The joint design of gradient and radiofrequency (RF) waveforms for small-tip 3D spatially tailored excitation poses a non-linear, non-convex, and constrained optimization problem that remains unsolved. Most existing approaches to tailored RF design either pre-define the gradients and only solve the least-squares RF design sub-problem, or restrict the gradients to certain classes such as echo-planar¹ or concentric shells² that permit low-order parametrization. Recently, a more general method for joint gradient and RF design was proposed³, in which the excitation k-space trajectory is expressed using a 2nd-order B-spline basis leading to the following constrained optimization problem:

$$ \arg\min_{\mathbf{c},\mathbf{b}} ||\mathbf{A}(\mathbf{c}) - \mathbf{b}||_2^2 + R(\mathbf{b}) \textrm{, subject to gradient hardware constraints,} \qquad \qquad \qquad (1) $$

where $$$\mathbf{c}$$$ is the B-spline coefficients, $$$\mathbf{A}$$$ is the small-tip system matrix⁴, and $$$\mathbf{b}$$$ is the RF waveform. The method in [3] alternates between updating $$$\mathbf{c}$$$ and $$$\mathbf{b}$$$, resulting in locally optimal gradients that are not restricted to the heuristic trajectories of previous approaches. However, the gradients obtained depend on the initial k-space trajectory that "seeds" the algorithm. For example, in [3] it was observed that a stack-of-spirals initial trajectory produced a relatively poor excitation (after local optimization) compared to the other initial trajectories tested. Here we propose an extension of [3] that reduces the dependence of the final excitation patterns on the initial (seed) k-space trajectory, possibly leading to designs that are closer to being globally optimal.

Our approach is based on the Iterated Local Search (ILS) method⁶, an iterative optimization strategy that employs a local optimization routine at each iteration. Our ILS search begins by finding a local optimum $$$k^*_0$$$ starting from some initial k-space trajectory (e.g., stack-of-spirals or concentric shells). We then attempt to sample "nearby" local minima by "perturbing" $$$k^*_0$$$ multiple times to produce a set of $$$N_{test}$$$ candidate seed trajectories $$$k_{1,j}$$$, $$$j=1...N_{test}$$$, for the next iteration (Fig. 1). For each $$$k_{1,j}$$$, Eq. (1) is descended (in parallel on a 32-core desktop computer) using the method in [3] to yield a local optimum $$$k^*_{1,j}$$$, and among all $$$k^*_{1,j}$$$, $$$j=1...N_{test}$$$, the trajectory producing the miminum cost function value is chosen as the value $$$k^*_1$$$ for the next iteration. This procedure is then repeated until convergence. We generate the candidate trajectories $$$k_{i+1,j}$$$, $$$j=1...N_{test}$$$, at iteration $$$i$$$ by expressing the locally optimal trajectory $$$k^*_i$$$ using a reduced number $$$N_{basis,j}$$$ ($$$<N_{basis}$$$) of B-spline basis functions, i.e., by projecting $$$k^*_i$$$ onto a set of $$$N_{test}$$$ lower-dimensional B-spline bases. For a given $$$j$$$ we do this projection by first obtaining an unconstrained least-squares fit $$$\mathbf{c}_j^u$$$ to $$$k^*_i$$$ using $$$N_{basis,j}$$$ 2nd-order B-spline basis functions, and then minimizing the following constrained cost function using Matlab's quadprog function:

$$\arg\min_{\mathbf{c}_j} ||\mathbf{c}_j - \mathbf{c_j^u}||_2^2 \textrm{, subject to gradient hardware constraints.}$$

Here we use $$$N_{basis}=50$$$ and $$$N_{basis,j} = \{ 10,12,14,...,50 \}$$$ (i.e., $$$N_{test}=20$$$, with the last value of 50 corresponding to the unperturbed trajectory). In essence, the $$$k_{i+1,j}$$$, $$$j=1...N_{test}$$$, correspond to smoothed versions of the current iterate $$$k^*_i$$$ with varying degree of smoothness. We evaluated three different k-space initializations: stack-of-spirals, SPINS⁵, and the extended-KT trajectory from [3]. The target pattern was a 3D cube. RF pulse duration was approximately equal (4ms) for all three k-space initializations. Simulations and experiments were done for single-coil transmission, however the principles introduced here would also apply to parallel transmission. We ran a total of five ILS iterations.

Figure 2 shows the design cost function value (Eq. (1)) versus computation time for both local and the proposed global (ILS) optimization, for the three different k-space initializations. Each curve shows the cost-vs-time for the best local minimum $$$k^*_i$$$ obtained at each iteration. A sharp spike along a curve indicates that the ILS algorithm identified a nearby local minimum with lower cost than the unperturbed trajectory. Stack-of-spirals produces significantly higher cost (poorer excitation accuracy) after local optimization compared to the other two initializations, however after the ILS search all three k-space intializations produce more similar final costs.

Figure 3 shows Bloch simulated excitation patterns obtained with both local and ILS optimization. All three k-space initializations produce similar excitation accuracy only after the proposed ILS optimization, as expected from Fig. 2. Figure 4 validates these excitation results experimentally, for the case of stack-of-spirals initialization [GE 3T scanner; uniform gel phantom; body coil RF transmission; 8-channel headcoil reception; matrix 240x240x48; FOV 24x24x24 cm³; TR=100ms; flip angle 7.5^o].

We have proposed an Iterated Local Search strategy for the joint design of gradients and RF waveforms for small-tip 3D tailored excitation, that reduces the influence of the choice of initial (seed) k-space trajectory on excitation accuracy.