GRAPE algorithm formulation for a single spin 1/2Optimization of point-to-point pulses

A shaped pulse can be seen as a sequence of N short pulses of length \(\Delta t\) with piecewise constant rf-amplitudes, where the \(j^}\) pulse normally consists of two controls \(\omega _}(j)\) and \(\omega _}(j)\) that represent the x and y-components of the shaped pulse. While the offset frequency \(\omega _}\) at a certain position in the spectrum relative to the irradiation frequency constitutes the free evolution or drift Hamiltonian

$$\begin }_0 = \omega _} \ I_z = 2 \pi \ \nu _} \ I_z, \end$$

(1)

the pulses constitute the control Hamiltonian at the \(j^}\) pulse

$$\begin& }_}(j) = \omega _}(j) \ I_x + \omega _}(j) \ I_y \\ & \quad \qquad= 2 \pi \ \nu _}(j) \. \end$$

(2)

For a given sequence of N pulses with duration \(\Delta t\) at a specific offset \(\nu _}\) the propagator is given by

$$\begin U_j = \exp \}_0 + }_}(j) \} \Delta t \} \end$$

(3)

and the propagation of an initial spin density operator \(\rho _0\) can be written as

$$\begin \rho _N = U_N \cdots U_j \cdots U_1 \ \rho _0 \ U^_1 \cdots U^_j \cdots U^_N . \end$$

(4)

The goal of a point-to-point optimization is to find values of the controls that minimize the differences between the desired final state of the spin (\(\mathrm _}\)) and the obtained final density operator with the current control amplitudes (\(\mathrm _N\)), which is identical to maximizing their overlap according to (Khaneja et al. 2005)

$$\begin \Phi _}= Re \left\langle }}\left| \right. }}\right\rangle \end$$

(5)

In order to maximize the cost function \(\Phi _}\), we have to minimize its gradient with respect to every control at every timestep, which for the \(j^\) time point is resulting in

$$\begin & \Gamma _}(j)=\begin\dfrac}} \\ \dfrac}} \end \nonumber \\ & \quad = \begin \left\langle \bigg | \rho _ U^_j}\right\rangle + \left\langle \bigg | \dfrac_j}}\right\rangle \\ \left\langle \bigg | \rho _ U^_j}\right\rangle + \left\langle \bigg | \dfrac_j}}\right\rangle \end\nonumber \\ & \quad \approx \begin\left\langle \left| \right. }_(j), \rho _j]}\right\rangle \\ \left\langle \left| \right. }_(j), \rho _j]}\right\rangle \end \end$$

(6)

with \(\lambda _j = U^_ \cdots U^_ \ \lambda _} U_N \cdots U_\), \(\rho _j = U_j \cdots U_1 \ \rho _0 \ U^_1 \cdots U^_j\) and \(}_(j)\) and \(}_(j)\) being the x and y components of the control Hamiltonian. Thus, for each iteration i of the optimization, the controls \(\omega ^}_}(j)\), \(k\in \\), are guaranteed to increase the cost function \(\Phi _}\) for an infinitesimal \(\epsilon\) according to

$$\begin \omega _}^(j) \rightarrow \omega _}^(j)+\epsilon \frac}}}(j)}, \end$$

(7)

until convergence to a (local) optimum is reached. Please note that the partial derivatives may also be calculated using the rotation angles

$$\begin \theta _x(j) = \omega _x(j) \, \Delta t \, \, ; \, \, \theta _y(j) = \omega _y(j) \, \Delta t. \end$$

(8)

In this case, the update can be rewritten into

$$\begin \theta _}^(j) \rightarrow \theta _}^(j)+\epsilon \frac}}}(j)}, \end$$

(9)

using the relation \(d\omega / d\theta = 1/\Delta t\). Instead of the GRAPE-approximation given at the end of Eq. 6, it is advantageous to use the computationally more costly exact gradients for the pulse sequence update (de Fouquieres et al. 2011). The critical point in these exact gradients is the calculation of \(\partial U_j / \partial \omega _x(j)\) or \(\partial U_j / \partial \theta _x(j)\), respectively, and its complex conjugate, for which several solutions have been proposed (Goodwin and Kuprov 2015; Levante et al. 1996), but elegant analytical solutions might be of particular interest for specific problems. For a single spin we will derive such very compact analytical forms for various optimization types in the following.

In order to design broadband PP pulses with \(B_1\) inhomogeneity compensation, it is further necessary to calculate the cost function and the gradient for \(n_}\) offsets linearly distributed over the desired bandwdith \(\Delta \nu\) and \(n_}\) different rf-amplitudes \(\nu _}\) covering the desired \(B_1\)-compensated range. The global quality factor and the global gradient can be calculated as averages according to

$$\begin \overline}}=\frac}n_}}\displaystyle \sum _^}}\sum _^}} \Phi _}(\nu ^i_},\nu ^l_}) \end$$

(10)

and

$$\begin \overline}(j)}=\frac}n_}}\sum _^}}\sum _^}} \Gamma _}(j,\nu ^i_},\nu ^l_}) . \end$$

(11)

Optimization of universal rotation pulses

If not the transformation of a given state to a final state, but rather the rotation in Hilbert space itself is the target, the propagator itself can be optimized (de Fouquieres et al. 2011; Khaneja et al. 2005; Kobzar et al. 2012; Luy et al. 2005; Skinner et al. 2012). Considering the total propagator U(T) at time point \(T = N \Delta t\) given by

$$\begin U(T) = U_N \cdots U_j \cdots U_1, \end$$

(12)

a cost function with respect to the desired propagator \(U_}\) can be formulated as

$$\begin \Phi _} = Re\left\langle }}\left| \right. \right\rangle , \end$$

(13)

and the corresponding gradients to the \(j^\) time step for x and y components of the pulse train s are approximated by (Khaneja et al. 2005)

$$\begin & \Gamma _}(j)=\begin\dfrac}} \\ \dfrac}} \end = \beginRe \left\langle \left| \right. X_}\right\rangle \\ Re \left\langle \left| \right. X_}\right\rangle \end\nonumber \\ & \quad \approx \beginRe \left\langle \left| \right. }_(j) \, X_j}\right\rangle \\ Re \left\langle \left| \right. }_(j)} \, X_j}\right\rangle \end \end$$

(14)

with \(P_j = U^_ \cdots U^_ \ U_}\), \(X_j = U_j \cdots U_1\) and previously used notations for the control Hamiltonian components \(}_(j)\) and \(}_(j)\).

Analytical exact point-to-point gradients using x, y, and optional z controls

The usual equations for evolution of a spin density matrix have been used in the previous section. However, it might be advantageous to use the Liouville superoperator formalism for a single spin \(\frac\) in Cartesian coordinate representation. The spin density operator can then be represented by a four vector \(\rho = (\rho _}, \rho _x, \rho _y, \rho _z)^T\), where the contribution of the identity matrix may be neglected. If furthermore relaxation is neglected, the spin density at time point j is fully represented by the vector

$$\begin \rho _}= \begin \rho _} \\ \rho _} \\ \rho _} \end \end$$

(15)

and propagation is achieved by the Liouville superoperators \(R_j\) in the reduced Cartesian representation

$$\begin \rho _}=R_}...R_}\rho _0 \end$$

(16)

and the co-state is represented accordingly by

$$\begin \lambda _j=R_}^...R_}^\lambda _}. \end$$

(17)

Please note, that the spin density as well as the co-state vector only have to be multiplied single-sided by the corresponding Liouville superoperators. Since we are in the Cartesian component basis set, we furthermore see that the superoperator is represented by a simple rotation matrix, which can be written as

$$\begin R_}= \begin \cos (\theta )+n_x^2(1-\cos (\theta )) & -n_z\sin (\theta )+n_xn_y(1-\cos (\theta )) & n_y\sin (\theta )+n_xn_z(1-\cos (\theta )) \\ n_z\sin (\theta )+n_xn_y(1-\cos (\theta )) & \cos (\theta )+n_y^2(1-\cos (\theta )) & -n_x\sin (\theta )+n_yn_z(1-\cos (\theta )) \\ -n_y\sin (\theta )+n_xn_z(1-\cos (\theta )) & n_x\sin (\theta )+n_yn_z(1-\cos (\theta )) & \cos (\theta )+n_z^2(1-\cos (\theta )), \end \end$$

(18)

with the overall rotation angle \(\theta\) and the normalized components \(n_x\), \(n_y\) and \(n_z\) of the rotation axis of timestep j. It can also be noted more compact using the Rodrigues formula using its individual elements \(R_\):

$$\begin R_= \cos ^\bigg (\dfrac\bigg )+(2n_h^2-1)\sin ^\bigg (\dfrac\bigg ) & \text \ddot}\text ~h=k \\ 2n_hn_k\sin ^\bigg (\dfrac\bigg )-\epsilon _~n_l\sin (\theta ) & \text \ddot}\text ~h\ne k, \end\right. } \end$$

(19)

where \(\epsilon _\) is the Levi-Civita symbol.

In each step j, \(\theta\), \(n_x\), \(n_y\) and \(n_z\) can be calculated from the control values \(\omega _x\), \(\omega _y\) and \(\omega _z\) as follows:

$$\begin \theta&=\Delta t\sqrt})^2} = \sqrt\end$$

(20)

$$\begin n_x&=\dfrac})^2}}=\dfrac} = \dfrac\end$$

(21)

$$\begin n_y&= \dfrac})^2}} =\dfrac} = \dfrac \end$$

(22)

$$\begin n_z&= \dfrac}}})^2}} =\dfrac}}} = \dfrac. \end$$

(23)

with the effective flip angles around the Cartesian axes \(\theta _x = \omega _x \Delta t\), \(\theta _y = \omega _y \Delta t\), \(\theta _z = (\omega _z + \omega _}) \Delta t\), where the rotation around the z-axis now consists of the offset term introduced above and the introduction of a potential z-control that cannot be applied directly on a spectrometer, but that can be used to allow direct implementation of pulse sequence bound offset changes in optimizations (Stockmann and Wald 2018; Vinding et al. 2017, 2021). With these equations in hand, it is straightforward to reformulate the cost function for point-to-point pulses in Cartesian Liouvielle representation as

$$\begin & \Phi _}=\left\langle }}\left| \right. }}\right\rangle =\left\langle }}\left| \right. }...R_}\rho _0}\right\rangle =\Big \langle R_}^...R_}^\lambda _}\Big | R_}...R_}\rho _0\Big \rangle \nonumber \\ & \quad =\Big \langle \underbrace}^...R_}^\lambda _}}_}}\Big | \underbrace}...R_}\rho _0}_}}\Big \rangle . \end$$

(24)

Accordingly, the gradients for \(j^\) time point for the single spin is

$$\begin \Gamma _}(j)=\begin\dfrac}} \\ \dfrac}} \\ \dfrac}} \end = \begin\left\langle \left| \right. \rho _}\right\rangle \\ \left\langle \left| \right. \rho _}\right\rangle \\ \left\langle \left| \right. \rho _}\right\rangle \end, \end$$

(25)

which leaves an analytical derivation of the derivative of the rotation matrix with respect to the Cartesian controls to obtain an overall exact gradient. With the formulae given in this section, the derivatives of the individual matrix components can be straightforwardly achieved and corresponding results are shown in Table 1 for the x andy derivatives and in Table 2 for the z-control.

Table 1 Derivatives of rotation matrix components with respect to \(\theta _\) and \(\theta _\) for point-to-point optimizationsTable 2 Derivatives of rotation matrix components with respect to \(\theta _\) for point-to-point optimizationsAnalytical exact point-to-point gradients using amplitude, phase, and optional z controls

The normalized components of the rotation axes \(n_x\), \(n_y\) and \(n_z\) can also be written in polar coordinates, resulting in the pulse phase \(\alpha\) and rf-amplitude \(\omega _}\) as the input parameters for pulse shapes. The corresponding notation is

$$\begin n_x= & \dfrac}=\dfrac}\Delta t}})^2+\theta _z^2}}, ~n_y =\dfrac}\nonumber \\ & =\dfrac}\Delta t}})^2+\theta _z^2}}, ~n_z =\dfrac=\dfrac})^2+\theta _z^2}} \end$$

(26)

with

$$\begin \theta _&=\omega _}\Delta t \nonumber \\ \theta _z&=(\omega _z+\omega _})\Delta t \nonumber \\ \theta&=\sqrt^2 + \theta _z^2} \end$$

(27)

The derivatives for the rotational matrix components with respect to \(\alpha\), \(\omega _}\) and \(\omega _z\) can then be derived the same way as in the previous section. The gradients for the \(j^\) time point for the single spin are

$$\begin \Gamma _}(j)=\begin\dfrac}} \\ \dfrac}}}(j)} \\ \dfrac}} \end = \begin\left\langle \left| \right. \rho _}\right\rangle \\ \left\langle \left| \right. }(j)} \rho _}\right\rangle \\ \left\langle \left| \right. \rho _}\right\rangle \end, \end$$

(28)

which again leaves an analytical derivation of the derivative of the rotation matrix with respect to the Cartesian controls to obtain an overall exact gradient. The calculation can be done by hand or a symbolic mathematics program and the result for the rotation matrix components are listed in Tables 3 and 4, respectively. The amplitude and phase representation is particularly useful in cases of restricted rf-amplitudes or even constant rf-amplitudes (Skinner et al. 2006), where in the latter case only derivatives to the phase \(\alpha\) need to be considered for a minimum set of parameters.

Table 3 Derivatives of rotation matrix components with respect to phase \(\alpha\) and rf-amplitude \(\theta _\) for point-to-point optimizationsTable 4 Derivatives of rotation matrix components with respect to z rotations in polar coordinates of the xy-plane for point-to-point optimizationsAnalytical exact universal rotation gradients

We have seen that propagation can be expressed in terms of simple rotations in the case of a single spin \(\frac\). If effective rotations themselves need to be optimized, it is best to express them with minimum storage and computation time. As such, Cayley (1997) or, equivalently, quaternions (Emsley and Bodenhausen 1990; Hamilton 1843; Rodrigues 1840; Slad et al. 2022) can be used to express rotations. Although the minimum set of numbers to represent a rotation is the vector spanned by \((\theta _x, \theta _y, \theta _z)\), it is better to use a four-dimensional vector

$$\begin Q_j = \begin A_j \\ B_j\\ C_j \\ D_j \end \end$$

(29)

for the rotation at time step j, where the four components are defined by

$$\begin A_}&= n_x(j) \ \sin (\frac),\end$$

(30)

$$\begin B_}&= n_y(j) \ \sin (\frac),\end$$

(31)

$$\begin C_}&= n_z(j) \ \sin (\frac),\end$$

(32)

$$\begin D_}&= \cos (\frac) \end$$

(33)

using the definitions of Eqs. (20)-(22). Quaternions have the advantage that a series of rotations can be directly evaluated via the product

$$\begin \cdot _= \begin + & - & + & +A_2 \\ + & + & -A_2 & + \\ - & +A_2 & + & + \\ -A_2 & - & - & + \end\begin _ \\ _ \\ _ \\ _ \end \end$$

(34)

which, with its 16 simple multiplications and sums, is computationally more efficient than the construction of a rotation matrix involving at least one sine/cosine calculation.

Using the formalism for universal rotation pulses, the cost function can be written directly as

$$\begin \Phi _}=Re\left\langle }}\left| \right. \right\rangle =Re\left\langle }}\left| \right. \right\rangle \end$$

(35)

and the gradients are derived as

$$\begin \Gamma _}(j)=\begin\dfrac}} \\ \dfrac}} \\ \dfrac}} \end = \beginRe \left\langle }}\left| \right. \cdot \dfrac \cdot Q_ \cdots Q_1 }\right\rangle \\ Re \left\langle }}\left| \right. \cdot \dfrac \cdot Q_ \cdots Q_1 }\right\rangle \\ Re \left\langle }}\left| \right. \cdot \dfrac \cdot Q_ \cdots Q_1 }\right\rangle \end, \end$$

(36)

where again the overall gradient calculation is reduced to simple quaternion propagation plus the calculation of the rotation derivative - this time in quaternion notation. With Eqs. (30)-(33) the involved partial derivatives of the individual quaternion components can be calculated and corresponding results are summarized in Table 5. Obviously, the components of interest can also be derived with respect to polar coordinates \(\alpha\) and \(\theta _\) in the xy-plane, for which the terms are given in the right column of the same Table. It should be noted that the derivatives with respect to \(\alpha\) are particularly simple, which guarantees fastest gradient calculation times for the case of constant amplitude pulses.

Table 5 Exact gradients of the quaternion elements A, B, C, D with respect to cartesian coordinates (\(\theta _x\), \(\theta _y\) and \(\theta _z\)) and polar coordinates (\(\alpha\),\(\theta _\) and \(\theta _z\)). Auxiliary variable is \(n_ = \frac}\), other variables as defined in the main textLimited rf-amplitudes as holonomic constraints

Optimizing efficient pulses in magnetic resonance usually implies optimizations with boundary conditions. One of the fundamental boundaries concern large rf-amplitudes that can cause experimental problems in the form of arching, violated duty cycles, and exceedance of allowed energy depositions. As such, they pose clear restrictions to applicable rf-amplitudes, rf-power, or rf-energy, respectively. Simple implementations with hard cut-off limits have been proposed early on in optimal control optimizations (Gershenzon et al. 2008, 2007; Kobzar et al. 2004, 2008; Skinner et al. 2004). However, mathematically more sound are so-called holonomic constraints that can be included in optimizations as Lagrange multipliers. Holonomic constraints essentially require continuously differentiable functions for restrictions and are usually implemented by auxiliary variables. A multitude of periodic functions have been proposed as auxiliary functions (Goodwin 2017), but we will concentrate here on a widely used trigonometric function, the hyperbolic tangent \(\tanh\), which can be incorporated as a reduced rf-amplitude in polar coordinate systems according to

$$\begin \theta _^}(j) = \theta ^}(j) \ \cdot \ \tanh \left( \frac(j)}}(j)} \right) \end$$

(37)

where the reduced amplitude \(\theta _^}(j)\) of the \(j^\) digit continuously converges only to the maximum value \(\theta ^}(j)\). In an optimization, \(\theta _(j)\) instead can now be changed as an auxiliary variable without any restriction, while \(\theta _^}(j)\) will be used for the actual values of rf-amplitudes in the pulse shape.

For a simple amplitude restriction to a j-dependent maximum value,

$$\begin \theta ^}(j) = \textrm (j) = 2 \pi \ \Delta t (j) \ \nu _}^}(j) \end$$

(38)

with maximum rf-amplitude \(\nu _}^}(j)\) can be imposed. Usually a constant restriction for all time steps is chosen, but in special cases, physics may imply a j-dependent upper limit function (Skinner et al. 2012). In an actual optimization, we restrict ourselves to the polar case where the rf-amplitude is inherently included in the variable \(\theta _\). However, using the auxiliary approach, the variables \(n_x\), \(n_y\), \(n_z\), \(\theta _\), and \(\theta\) from equations (18,19,26,27,29 - 33) used in the rotation matrices as well as in the quaternion formalism, have to be replaced by their reduced variants according to

$$\begin \theta _^}(j)= & \theta ^}(j) \ \cdot \ \tanh \left( \frac(j)}}(j)} \right) \end$$

(39)

$$\begin \theta ^}(j)= & \sqrt^}(j))^2 + \theta _z(j)^2} \end$$

(40)

$$\begin n_x^}(j)= & \frac^}(j)}}(j)} \end$$

(41)

$$\begin n_y^}(j)= & \frac^}(j)}}(j)} \end$$

(42)

$$\begin n_x^}(j)= & \frac(j)}}(j)} \end$$

(43)

The resulting gradients \(\Gamma _}\) and \(\Gamma _}\) are then essentially of the same form as derived above. Only for the derivation to the auxiliary variable itself the \(\tanh\) is not a constant and care has to be taken. Taking the derivation of the \(j^\) element of an UR optimization as an example, we can derive

$$\begin \frac}} (j)} = \underbrace}}^} (j)}}_} \, \cdot \, \underbrace^} (j)} (j)}}_} \end$$

(44)

With the auxiliary derivative term the actual control \(\theta _^} (j)\) (i.e. the rf-amplitudes of the shaped pulse) can be used, for the optimization, instead, the variable \(\theta _\) has to be used, making it necessary to include the terms

$$\begin \frac^} (j)} (j)} = \text ^2 \left( \frac(j)}}(j)} \right) = 1 - \tanh ^2 \left( \frac(j)}}(j)} \right) \end$$

(45)

into the equations of the gradient in question. Please note that the second solution is derived from the general relation \(\text ^2(x) + \tanh ^2(x) = 1\) and results in a computationally friendly term, as \(\tanh \left( \frac(j)}}(j)} \right)\) has to be calculated anyway.

If overall rf-power restrictions have to be imposed, the maximum rotation angle can also be defined via the reduced angle at time point j according to

$$\begin \theta _^}(j) = \theta _(j) \ \cdot \ \sqrt}}}}} \ \cdot \ \tanh \left( \sqrt}}}}} \right) \end$$

(46)

with the squareroots of the maximum allowed average rf-power \(\overline}}\) and actual average rf-power expressed by

$$\begin \overline = \sum _^N \frac(j))^2} \end$$

(47)

where piecewise constant rf-amplitudes and uniform time steps \(\Delta t\) have been assumed. Again, the derivative of the reduced rotation angle has to be calculated for the gradient, but this time the reduction involves rotation angles at all time steps and the derivative in e.g. an UR optimization with rf-power restriction will involve

$$\begin \frac}} (j)} = \sum _^N \ \frac}}^} (k)} \, \cdot \, \frac^} (k)} (j)} \end$$

(48)

with the reduced rotation angle derivatives defined by

$$\begin \frac^}(j)} (j)}= & \left( 1 - \frac(j))^2}^N (\theta _(j))^2} \right) \cdot \sqrt}}}}} \tanh \left( \sqrt}}}}} \right) \nonumber \\ & + \frac(j))^2}^N (\theta _(j))^2} \left( 1 - \tanh ^2 \left( \sqrt}}}}} \right) \right) \ \ \ \ \forall \ \ \ k=j \;\end$$

(49)

$$\begin \frac^}(k)} (j)}= & \frac(j))^2}^N (\theta _(j))^2} \cdot \sqrt}}}}} \tanh \left( \sqrt}}}}} \right) \nonumber \\ & + \frac(j))^2}^N (\theta _(j))^2} \left( 1 - \tanh ^2 \left( \sqrt}}}}} \right) \right) \ \ \ \ \forall \ \ \ k \ne j \ . \end$$

(50)

Finally, also the overall rf-energy of a pulse can be restricted, using, for example, the energy defined in Hz

$$\begin \frac = \overline \, t_p, \end$$

(51)

or the more convenient expression in terms of rotation angles

$$\begin E_ = \sum _^N (\theta _(i))^2 = (2 \pi )^2 \ \Delta t \ \frac \end$$

(52)

when equations (46 - 50) are used with substituting \(\overline}}\) with the maximum allowed energy \(E_^} = (2 \pi )^2 \Delta t E^} / h\) and \(\overline\) by \(E_\).