Abstract
In medical magnetic resonance imaging, parallel radio frequency excitation pulses have to respect a large number of specific absorption rate constraints. Geometrically, each of these constraints can be interpreted as a complex, centered ellipsoid. We propose to replace a collection of such constraints by the single constraint which corresponds to the associated maximum volume inscribed ellipsoid and implies all original constraints. We describe how to compute this ellipsoid via convex programming. Examples show that this reduction has very short computation times but cuts away parts of the feasible power domain.
1. Introduction
Parallel radio frequency (RF) excitation pulses for magnetic resonance imaging (MRI) deposit energy in the irradiated tissue. In the affected areas this creates an amount of heat which depends on the strength and duration of the superimposed RF fields and the conductivity and density of the irradiated tissue. Specific absorption rate (SAR) constraints are responsible for restricting this heating to a harmless level. Current International Electrotechnical Commission [1, 2] and Food and Drug Administration [3] regulations define a large number of local SAR constraints, i.e., one for each 10 gram volume equivalent. Medical RF excitation pulses not only have to produce the desired target magnetization with the scanner hardware capabilities, but also must make sure that every mandated SAR constraint is respected.
Physics dictates that the thermal energy deposited in a specific volume by oscillating electromagnetic fields is given by a time-homogeneous quadratic functional applied to the RF excitation waveforms; compare equation (4) in [4] or equation (1) in [5]. A SAR constraint sets an upper bound for a particular quadratic functional as in definition (2) below. The parameters for this representation can be obtained by measurements [6–8].
In the case of a single RF excitation coil the SAR constraints can easily be reduced to the single constraint which belongs to the most affected specific volume. However this is not possible in the case of parallel excitation with more than one coil, where the most affected specific volume can depend on the utilized coil voltages and phases. In this case the SAR constraints can present a serious additional burden for the design of RF waveforms, especially because of their large number. This obstacle might not only increase pulse design computation time, but also impact the quality of the calculated excitation pulse.
Special algorithms have been developed for coping with parallel RF pulse design problems in the presence of SAR constraints. In [4] the pulse design problem with SAR constraints is addressed via a quadratic small tip angle approximation [9] and a dimension-reduction method. However, only three SAR constraints are addressed in the numerical example. In [10], Sbrizzi et al. suggest a reduction of the number of SAR constraints through a selection via largest eigenvalues and clustering by a sufficiently large correlation of the associated eigenvectors. The remaining 12 to 36 SAR constraints in the pulse design examples are dealt with by Lagrange multipliers and a conjugate gradient algorithm in a small tip angle approximation. Guerin et al. [11] have implemented a primal-dual interior point algorithm for optimizing parallel RF sinc-pulses with 8 channels, up to 11 spokes and about 1300 SAR constraints. Several types of nonlinear algorithms also suited for large flip angles have been investigated by Hoyos-Idrobo et al. [12]. They were capable of handling RF pulse design problems with 8 parallel channels for 5 and 7 points and about 500 SAR constraints. That work has been extended to also include the optimization of the points in [13]. All mentioned algorithms have in common that they initially reduce the number of SAR constraints, in order to arrive at a number of constraints which is manageable by the utilized pulse design algorithm. The dominant technique for handling the remaining SAR constraints is Lagrange multiplier, which are updated during the RF pulse optimization. A notable exception is the algorithm in [14] which minimizes the maximum relative SAR energy subject to restrictions for the mean square error of the excited magnetization profile in several slices without a reduction of the SAR constraints.
A popular SAR reduction technique was introduced by Eichfelder and Gebhardt in [15]. The SAR matrices are sorted by eigenvalues and -vectors. Families of comparable SAR matrices are replaced by dominating virtual observation points VOPs. Formally these VOPs are identical to SAR constraints, but they are not necessarily connected to a physical heating of a specific volume of the irradiated tissue. In this work we will call them virtual SAR constraints. Lee et al. [5] have further refined this reduction technique. These reductions typically yield several hundreds of virtual SAR constraints, depending on the choice of an overestimation parameter. This overestimation parameter controls the maximum RF power reduction for the (dominating) virtual SAR constraints compared to the (dominated) original ones.
In this work we present a novel pure and extreme reduction technique. Geometrically, the quadratic functional defining a SAR constraint can be interpreted as a complex, centered ellipsoid with a dimension corresponding to the number of RF excitation coils. We propose to replace a collection of SAR constraints by the single virtual SAR constraint, which is associated with the maximum volume inscribed ellipsoid (MVIE). This single MVIE SAR constraint implies all original SAR constraints. The maximum volume characteristic makes sure that maximum room remains for the RF pulse design after this reduction. This intuitive idea has the appealing mathematical property that the parameters of the MVIE SAR constraint can be obtained as the solution of a convex optimization problem. We take advantage of this fact to develop methods for a quick numerical calculation of the parameters of the MVIE SAR constraint. The developed “implicit active constraints selection” heuristic in the iterations of an interior point algorithm shows particularly short numerical run times.
Besides these algorithmic aspects, this work focuses on the mathematics, which underlies the MVIE SAR constraints reduction. To this end we strengthen relevant theoretical results for MVIEs from the real-valued case to the complex-valued one, as required in our application; see the Appendix. Noticeable is the result that blowing up the complex MVIE ellipsoid by a factor corresponding to the square root of its dimension yields an ellipsoid, which contains the intersection of all original ellipsoids, which, to the best of our knowledge, has not been proved for the complex case before.
We present only few tentative numerical results. The two main examples consist of collections of 448480 and 111939 SAR constraints for 2 and 8 RF excitation coils, respectively. These examples show not only the feasibility of our method, but also its speed: it only takes seconds to reduce the example SAR constraint collections to the MVIE SAR constraint on a standard personal computer.
As drawback of our method we have to mention the missing ability to control how much power the MVIE SAR constraint is cutting away from the feasible RF space of the original SAR constraints. Such a reduction can impact the achieved flip angle or lead to an increased inhomogeneity of the magnetization (compare Figure 7 in [5]). The MVIE has contact points with the border of the intersection of the original SAR ellipsoids (see Theorem A.2 in the Appendix and Figure 1). However, away from these contact points, the MVIE does not allow as much RF power as the original SAR constraints. One can show that this power reduction is upper bounded by the number of excitation coils (see display (12) and Theorem A.3 in the Appendix). Figure 1 gives the visual impression that the practical average and maximum power reductions might be much smaller than this theoretical bound. However in our examples, a sampling method for the estimation of the practically incurred power reduction yields disappointingly large values; compare Table 2. Nevertheless there might exist specific pulse design applications, where the advantage of the reduction of the large number of SAR constraints to a single one outweighs the reduction of the feasible RF power in the associated RF pulse design problems. We demonstrate this possible advantage with an example of [12], in which the reduction of their 491 local and global SAR constraints to the single MVIE SAR constraint can substantially reduce the RF pulse optimization time without a perceivable deterioration of the obtained pulse quality.
An overview of this work is as follows: we formally introduce SAR constraints in Section 2.1 and present the maximum volume inscribed ellipsoid reduction technique in Section 2.2. The numerical implementation of the logdet objective function is discussed in Section 3.1 and of the positive semidefiniteness constraints are discussed in Section 3.2. Section 3.3 shows how the large number of constraints in the optimization problem can be dealt with. In Section 4.1 we present several MVIE SAR calculations and in Section 4.2 we exemplify the influence of the SAR reduction on a pulse optimization problem. We conclude the work in Section 5. In the Appendix we extend the theory of maximum volume inscribed ellipsoids from the real to the complex case.
2. Theory
2.1. SAR Constraints
We let be the number of RF excitation coils. We let be the set of -dimensional, complex Hermitian matrices, that is, matrices which satisfy for every , where denotes the conjugate of a complex number . For a column vector we let be the corresponding row vector with for every . We use the notation as a shorthand for the complex scalar product of two vectors . A matrix is positive semidefinite, if for every vector (since is Hermitian, the value of the sum is real). The matrix is positive definite, if additionally “” holds in this display whenever . We denote the set of positive semidefinite (resp., positive definite) matrices in with (resp., ).
We let be the duration of the excitation pulse. A complex -dimensional vector-valued function is called RF waveform. It specifies the amplitude and phase for each time point and RF coil of the excitation pulse.
We consider SAR constraints. Each SAR constraint is characterized by a positive definite Hermitian matrix . An RF waveform satisfies SAR constraint ifHere we have scaled the matrix in such a way that the upper bound has the value one. An RF waveform will be called SAR-compliant if it satisfies (2) for every . More background on SAR constraints can be found in [4, 5, 15].
The RF pulse design problem is to find a SAR-compliant RF waveform, which accurately produces the desired magnetization and also respects the physical capabilities of the scanner hardware; see [12, 16] and the example of Section 4.2.
If value and speed of an RF waveform are scaled by a factor (i.e., is replaced by and by ) the SAR power on the left hand side of (2) becomes times its previous value. If the gradient trajectory during this RF pulse is analogously scaled, the excited magnetization profile remains unchanged (under the condition that B0-inhomogeneity and T1- and T2-relaxations are negligible). In this way one can reduce SAR power at the expense of an increased RF pulse duration.
2.2. Maximum Volume Inscribed Ellipsoid
For a positive definite Hermitian matrix we define an associated -dimensional, complex, centered ellipsoid by For two matrices the ellipsoid is contained in if and only if is positive semidefinite, i.e.,For later use we also note the equivalence
If an RF waveform and a matrix satisfyand for an , then satisfies SAR constraint . This follows from equivalence (4). In particular, if the RF waveform and the matrix satisfy (6) andthen is SAR-compliant. Condition (7) says that the ellipsoid associated with must be contained in each ellipsoid which is associated with one of the SAR matrices for .
It is well known that there exists a unique MVIE for each closed compact convex set with a nonempty interior [17–19]. In the Appendix we tersely state and justify the necessary underlying theoretical results because the complex set-up which is needed in our work seems not to be covered by the literature. The intersection of the SAR ellipsoids is such a set. Since it is symmetric, its MVIE is also symmetric and can be represented as for a matrix . The main idea of this work is to calculate this MVIE SAR matrix and use it to replace the SAR constraints (2) in the RF pulse design problem by the single virtual SAR constraint (6). The “maximum volume” property makes sure that there is maximum possible room for the RF pulse design. In Figure 1 this idea is illustrated for three two-dimensional real ellipsoids.
For a matrix the volume of the ellipsoid is given by the expression (compare, e.g., [19] equation ) where is the volume of the -dimensional complex unit ball.
Finding the maximum volume ellipsoid can be casted as a convex optimization problem: Since the negative of the concatenated functions and is convex on the convex set [18], the following mathematical program is a convex optimization problem:According to Theorem A.1 in the Appendix below, this problem possesses a unique solution which we denote . Its inverse has property (7). This follows from condition (11) and the equivalence (5). Moreover the ellipsoid of this matrix has the maximum volume among all such matrices, since the objective function (9) is a strictly increasing function of the volume of the ellipsoid which is associated with .
It is interesting to note that when the maximum volume inscribed ellipsoid is blown up by the factor then it encompasses the intersection of all original SAR ellipsoids,This statement is proved in Theorem A.3 in the Appendix below. By performing an RF pulse design with the single SAR constraint corresponding to and comparing it to the one with , one can obtain an upper bound for the loss of the magnetization quality owing to the passage from the original SAR constraints (2) to the single MVIE SAR constraint (6).
3. Implementation
3.1. -Function for Positive Definite Matrices
We unify the objective function (9) and condition (10) into the concave function by setting Our numerical implementation of the calculation of firstly attempts a Cholesky decomposition of the matrix . When this decomposition fails the matrix is not positive definite and the value of is . Otherwise the Cholesky decomposition yields a lower triangular matrix with has strictly positive real diagonal elements and satisfies . Thus the value of can be calculated by The lower triangular matrix can be inverted by forward substitutions. Thus the inverse of can be obtained as . This inverse is important for the calculation of the first- and second-order derivatives of [18, 20, 21].
The calculation of for a matrix requires less than multiply, add, or subtract operations, plus square roots, inversions, and logarithms. The calculation of the first-order derivatives of requires less than additional multiply, add, or subtract operations and inversions. The calculation of the second-order derivatives of takes about additional multiply, add, or subtract operations.
We note that the numerical calculation of the inverses for appearing in condition (11) can also be done via a Cholesky decomposition, inversion of the resulting lower triangular matrix, and product of this triangular matrix with its transpose; compare [22].
3.2. Handling of Positive Semidefiniteness Constraints
The idea for dealing with convex positive definiteness or semidefiniteness constraints or as in (10) and (11) is to replace them by the stronger convex constraint where is a small number. The -function frequently acts as a barrier function for positive semidefinite constraints [20, 21].
By modifying constraints (11) in this way we obtain the nonlinear convex optimization problem with individual SAR constraint approximationsIf for every , this problem has again a unique solution. The inclusion,shows that by settingthis problem is a strengthening of problem (9). Its solution is denoted by and we set .
The converse inclusion,shows that by setting in (17) we obtain a relaxation of problem (9). Its solution is denoted by . Thus we get which implies with and converse inequalities for the volumes of the associated ellipsoids. In this way, we can assess the fraction of the maximum volume ellipsoid which is lost by the strengthening from to and adjust the parameter . We foundto be a suitable value in the examples of Section 4. As usual denotes the sum of the diagonal elements of a square matrix.
3.3. Reduction of Slack Variables
Interior point algorithms [18], which are well suited for solving convex optimization problems in general, introduce a slack variable for each of the constraints. Thus, an augmented linear system with variables must be solved in each iteration. This system contains first- and second-order derivatives at the current iterate. Since can be much larger than the solution of these linear problems can dominate the overall solution time. In order to reduce computation time we therefore propose two variants of the nonlinear program (16).
In the first variant we replace the SAR constraints (11) by a single SAR constraint approximation:This is a convex optimization problem with variables and one constraint. The step computation thus becomes much faster, but, as a drawback of this variant, we observed more step-backs of the solver, that is, step size reductions in order to ensure that every matrix remains positive definite. Clearly, condition (17) implies (26) which implies (11) for . If is the solution of this problem we therefore have Again, converse inequalities can be obtained for the volumes of the associated ellipsoids.
The second variant is what we call an external active constraint selection. Here we fix a number much smaller than , but in the order of the number of variables . For instance, in the examples we use The nonlinear solver is told that there are only constraints. For a given iterate during the optimization process, we calculate the values for every . The constraints are sorted by these values and the solver gets values and derivatives of the first constraints in ascending order.
The constraints which are closest to being infeasible thus become the active constraints, but the solver is not aware that the active constraints might leave, enter, or exchange positions. The active constraints are selected externally. If the solver converges to a feasible solution, this solution satisfies all original constraints.
Since, after the evaluation and sorting steps, constraints are discarded, first- and second-order derivatives do not have to be calculated for them which reduces computation time. Moreover, in the step computations the augmented linear system has only rather than variables. In the examples of Section 4.1, this algorithm variant calculates a numerical approximation of the MVIE SAR parameters with the shortest computation time. This heuristic might be useful in other applications with a large number of constraints and a small number of variables, but there is no guarantee that it converges.
4. Examples
The following numerical examples show that MVIE SAR parameters can be calculated quickly and that an overall approach consisting of the calculation of the MVIE SAR constraint and its subsequent utilization in an RF pulse optimization can reduce the overall computation time without affecting the excitation quality. However, these examples are not comprehensive enough to make statements to which MRI applications these observations can be carried over.
The calculation times which we report in this section have been obtained with an Intel(R) Core(TM) i7-4800MQ CPU at 2.70 GHz. We are using the interior point solver IpOpt [23] as workhorse for solving the considered nonlinear optimization problems.
In order to assess how much RF power can be lost by reducing to (or its numerical approximations) we define the power loss in direction for as From these fractions we deduce the corresponding maximum, average, and minimum power losses by settingwhere denotes the uniform probability distribution on the complex unit ball . From the properties of we know Since there are no general closed form expressions for these quantities we are going to estimate them by taking the maximum (resp., average, resp., minimum) of fractions of a large number of directions sampled from .
4.1. Reduction of SAR Matrices
Here we consider three example sets of SAR matrices:(1)A Siemens internal example with RF coils and SAR constraints.(2)An example of Nicolas Boulant with RF coils and SAR constraints.(3)The example from Hoyos-Idrobo et al. [12] with RF coils and SAR constraints.
In Section 4.2, the rather small example (3) is further utilized to investigate how a pulse design optimization behaves when the original SAR constraints are replaced by the MVIE SAR constraint.
For each of these examples we calculate a virtual SAR constraint with each of the three different methods of Section 3:(1)By solving the convex optimization problem (16) with the individual SAR constraints approximation (17).(2)By solving variant (25) with the single SAR constraint approximation (26).(3)As the first method, but with the external active constraint selection for active constraints as explained in Section 3.3.
Each of these methods delivers an approximation of the MVIE SAR parameters which imply the original SAR constraints. This has been checked by performing a Cholesky decomposition of for every . The utilized value of in these approximations is given in (24). The optimization is stopped after at most 100 iterations. These and further characteristics of the calculations are gathered in Table 1.
In examples (1) and (2) with a large number of SAR matrices the solution of the problem with individual SAR constraints utilizes most of the computation time for the symbolic factorization of the Hessian matrix. This part of the computation is negligible in the other proposed algorithms. The external active constraint selection method largely reduces the number of back steps in the first example but increases this number in the other examples. As intended it spends little time for the calculation of the Hessians and achieves the smallest overall computation time.
Up to numerical inaccuracies all methods deliver the same maximum volume inscribed ellipsoid. In Table 2 we present the power losses obtained by uniformly distributed random samples from -dimension complex unit vectors. Clearly these numbers bring the major drawback of the MVIE to light: in certain directions, it cuts away a substantial part of the available RF power. Also, the power reduction in average is disappointingly large. Finally, it is astonishing that the observed minimum power fractions of the examples with 8 excitation coils remain about 10% away from their theoretical value 1.
4.2. Excitation Pulse Design
We consider the inversion pulse design problems of [12] with either 5 or 7 points [24], a fixed total duration of 1.39 and 4.03 milliseconds, and a target flip angle of 30 and 180 degrees, respectively. There are RF excitation coils. Each coil must satisfy a maximum voltage restriction during each point and additionally an average power bound. The 491 SAR constraints are those considered in the third example of the previous section. These constraints can be reduced to the single MVIE SAR constraint in less than one second, as shown in the last column of Table 1. In this section we investigate the effect of this reduction on the performance of an RF pulse optimization.
For this optimization we follow the approach described in [16]: the interior point solver IpOpt is used to find the RF amplitudes and phases during the points which minimize the magnitude least squares distance between the achieved and desired magnetization. Both, first- and second-order derivatives are utilized in this process. Two different approaches are available for the calculation of the achieved magnetization:(1)The arbitrary flip angle optimization (AFA) based on a solution of the Bloch equation without relaxation terms.(2)The small tip angle approximation (STA) amounting to a linearization of the solution of the Bloch equation.
We are applying these two approaches to each of the two pulse design problems with either the full set of 491 SAR constraints or the single MVIE SAR constraint obtained in Section 4.1. In Table 3 the observed computation times, number of iterations, and achieved root mean square errors are reported. The middle column shows the results for the original SAR constraints and the right-hand column those for the MVIE SAR. However, a certain reduction of the available RF pulse power is incurred by our method.
The use of the MVIE SAR reduces the computation time by a factor of about 1/2 to 1/3 in the STA and AFA optimization. This reduction is larger than the time it takes to calculate the MVIE SAR parameters. Only in the 30-degree excitation example with 5 points the root mean square error of the magnetization is slightly reduced from about 1.6 degrees to 1.9 degrees. In the 180-degree inversion example the excitation pulse quality is not affected, since the SAR constraints are not the limiting constraints, but the efficient power constraints.
This RF pulse design example shows that even in the case where the SAR matrices have been already been reduced to about 500 they put a burden on the pulse optimization which dominates computation time. This burden is largely reduced by utilizing the MVIE SAR constraint.
5. Conclusions
Parallel RF excitation pulses not only have to form an accurate magnetization, but also must make sure that a large number of SAR constraints are respected. In this work we suggested to reduce this additional burden by replacing a collection of SAR constraints by a single virtual MVIE SAR constraint. This MVIE SAR constraint is characterized by the properties that it implies each of the original SAR constraints and that its associated ellipsoid has the maximum volume. The parameters of the MVIE SAR constraint can be obtained by solving a convex optimization problem. Taking advantage of the properties of the -function, which appears in the objective function of the problem and can also be used as a barrier function for the positive semidefiniteness constraints therein and techniques for reducing the number of slack variables, one can utilize interior point algorithms to obtain solutions for realistic problem sizes within seconds. Our tentative examples show the feasibility of the approach. Also, the RF pulse optimization can get much faster when the tighter MVIE SAR is respected rather than the large number of original SAR constraints. In general, this tightening can lead to a certain degradation of the magnetization quality, except for special cases as in the example of Section 4.2.
Appendix
Complex Maximum Volume Inscribed Ellipsoids
The theory of maximum volume inscribed ellipsoids has been developed for ellipsoids in Euclidian spaces [17]. In this appendix we shortly justify the analogous results for the space which were used in the main part of this work. Throughout this section we let be a closed, compact, convex, and symmetric subset of having a nonempty interior. Symmetry means that for every and with we have , where . We let be the border of and the border of . For two vectors we let be the matrix defined by for .
Theorem A.1. There exists a unique matrix such that is the MVIE of , i.e., maximizes the volume subject to the condition .
Proof. As in [19] Section 12.4, one shows that the set of matrices satisfying is compact and convex. The strict concavity of the -function on implies the existence and the uniqueness of a matrix in this set which maximizes the value of . Since the -function is strictly increasing with the volume of the ellipsoid , the inverse has the requested properties.
Theorem A.2. If the matrix is such that is the MVIE of , there exists a number , multipliers , and vectors such that
Proof. As preparation we define the mapping which transforms a real matrix to a complex Hermitian matrix by setting for where denotes the imaginary unit. This means that the diagonal and strictly lower triangular values of the real matrix specify the real parts of the entries of the Hermitian matrix , and the strictly upper triangular values of the imaginary parts of . The mapping allows us to parameterize with real parameters. The set is an open subset of .
Moreover, we define for a nonempty set the convex function by where denotes the real part of a complex number . This function has the property, that for nonempty convex sets we have (compare Corollary 12.2.1 in [25]) Since is homogeneous (i.e., for all and ) this equivalence remains valid if its right-hand side is restricted to vectors with . Moreover, one can show that for We can thus cast the problem to find the matrix such that is the maximum volume inscribed ellipsoid for as Except for the correspondence which transfers the MVIE problem from the space to this reformulation is the same as in display of [19]. As in the real case one can apply Theorem I of [17] for a given solution at this point: it gives the existence of a number , coefficients , and directions with for such that for and for Setting this implies for Because the matrices and are Hermitian, this is equivalent toThe proof now proceeds as the one of Theorem 12.12 of [19]: taking the trace on both sides of (A.10), we get . By rescaling the -values we can therefore assume without loss of generality . Defining and for and utilizing (A.10) we thus obtain Moreover, and yield for every . Substituting thus completes the proof.
Theorem A.3. The MVIE matrix of Theorem A.1 satisfies
Proof. According to Theorem A.2 there exist , and such that (A.1) is satisfied. Since and are symmetric we additionally have for every and with . Moreover, for every , with and the Cauchy-Schwarz inequality gives Hence is a tangent hyperplane for at the point . Since , is convex and this implies that both and lie in the closed affine half-space . If we define the set by we therefore conclude .
We let be the positive definite square root of . Using (A.1) twice we obtain for Hence which completes the proof.
Data Availability
The data sets are not publicly available.
Conflicts of Interest
The author declares no conflicts of interest.