# Quick Reduction of Specific Absorption Rate Constraints in Parallel RF Excitation Pulse Design by Maximum Volume Inscribed Ellipsoids

Authors: Kurt Majewski
Journal: Concepts in Magnetic Resonance Part A
 (2019)

---

## 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.


---

## Body

## 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 kT points and about 500 SAR constraints. That work has been extended to also include the optimization of the kT 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 dominatingvirtual 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 themvirtual 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 anoverestimation 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.Figure 1
Maximum surface inscribed ellipse for three ellipses in light gray. Cut away difference between this ellipse and the intersection of the original ellipses in dark gray. Contact points are marked with circles and cut away corners with squares.An overview of this work is as follows: we formally introduce SAR constraints in Section2.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 letC∈N be the number of RF excitation coils. We let H be the set of C-dimensional, complex Hermitian matrices, that is, matrices M∈CC×C which satisfy Mi,j=Mj,i∗ for every i,j∈{1,…,C}, where α∗ denotes the conjugate of a complex number α∈C. For a column vector v∈CC we let v∗ be the corresponding row vector with (v∗)i:=(vi)∗ for every i∈{1,…,C}. We use the notation u∗v as a shorthand for the complex scalar product of two vectors u,v∈CC. A matrix M∈H ispositive semidefinite, if for every vector v∈CC(1)v∗Mv≔∑i,j=1Cvi∗Mi,jvj≥0(since M is Hermitian, the value of the sum is real). The matrix M ispositive definite, if additionally “>” holds in this display whenever v≠0. We denote the set of positive semidefinite (resp., positive definite) matrices in H with H+ (resp., H++).We letT>0 be the duration of the excitation pulse. A complex C-dimensional vector-valued function w:[0,T]→CC is calledRF waveform. It specifies the amplitude and phase for each time point and RF coil of the excitation pulse.We considerS∈N SAR constraints. Each SAR constraint s∈{1,…,S} is characterized by a positive definite Hermitian matrix A(s)∈H++. An RF waveform w satisfies SAR constraint s∈{1,…,S} if(2)∫0Twt∗Aswtdt≤1.Here we have scaled the matrix A(s) in such a way that the upper bound has the value one. An RF waveform w will be calledSAR-compliant if it satisfies (2) for every s∈{1,…,S}. 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 waveformw are scaled by a factor α>0 (i.e., w(t) is replaced by αw(αt) and T by T/α) 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 matrixH∈H++ we define an associated C-dimensional, complex, centered ellipsoid E(H)⊂CC by (3)EH≔v∈CC:v∗Hv≤1.For two matrices G,H∈H++ the ellipsoid E(G) is contained in E(H) if and only if G-H is positive semidefinite, i.e.,(4)EG⊂EH⇔G-H∈H+.For later use we also note the equivalence(5)G-H∈H++⇔H-1-G-1∈H++.If an RF waveformw and a matrix X∈H++ satisfy(6)∫0Twt∗Xwtdt≤1,and E(X)⊂E(A(s)) for an s∈{1,…,S}, then w satisfies SAR constraint s. This follows from equivalence (4). In particular, if the RF waveform w and the matrix X∈H++ satisfy (6) and(7)EX⊂⋂s=1SEAs,then w is SAR-compliant. Condition (7) says that the ellipsoid associated with X must be contained in each ellipsoid which is associated with one of the SAR matrices A(s) for s∈{1,…,S}.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 ⋂s=1SE(A(s)) of the SAR ellipsoids is such a set. Since it is symmetric, its MVIE is also symmetric and can be represented as E(X) for a matrix X∈H++. The main idea of this work is to calculate this MVIE SAR matrix X and use it to replace the S 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 matrixH∈H++ the volume of the ellipsoid E(H) is given by the expression (compare, e.g., [19] equation (12.5)) (8)VH≔VCdetH-1,where VC is the volume of the C-dimensional complex unit ball.Finding the maximum volume ellipsoid can be casted as a convex optimization problem: Since the negative of the concatenated functionslog and det is convex on the convex set H++ [18], the following mathematical program is a convex optimization problem:(9)maximizelog⁡det⁡Y(10)subjecttoY∈H++,(11)As-1-Y∈H+fors∈1,…,S.According to Theorem A.1 in the Appendix below, this problem possesses a unique solution which we denote YMVIE∈H++. Its inverse XMVIE≔YMVIE-1∈H++ has property (7). This follows from condition (11) and the equivalence (5). Moreover the ellipsoid E(XMVIE) of this matrix has the maximum volume among all such matrices, since the objective function (9) is a strictly increasing function of the volume V(Y-1) of the ellipsoid which is associated with Y-1.It is interesting to note that when the maximum volume inscribed ellipsoidE(XMVIE) is blown up by the factor C then it encompasses the intersection of all original SAR ellipsoids,(12)CEXMVIE=EXMVIEC⊃⋂s=1SEAs.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 XMVIE/C and comparing it to the one with XMVIE, 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).


## 2.1. SAR Constraints

We letC∈N be the number of RF excitation coils. We let H be the set of C-dimensional, complex Hermitian matrices, that is, matrices M∈CC×C which satisfy Mi,j=Mj,i∗ for every i,j∈{1,…,C}, where α∗ denotes the conjugate of a complex number α∈C. For a column vector v∈CC we let v∗ be the corresponding row vector with (v∗)i:=(vi)∗ for every i∈{1,…,C}. We use the notation u∗v as a shorthand for the complex scalar product of two vectors u,v∈CC. A matrix M∈H ispositive semidefinite, if for every vector v∈CC(1)v∗Mv≔∑i,j=1Cvi∗Mi,jvj≥0(since M is Hermitian, the value of the sum is real). The matrix M ispositive definite, if additionally “>” holds in this display whenever v≠0. We denote the set of positive semidefinite (resp., positive definite) matrices in H with H+ (resp., H++).We letT>0 be the duration of the excitation pulse. A complex C-dimensional vector-valued function w:[0,T]→CC is calledRF waveform. It specifies the amplitude and phase for each time point and RF coil of the excitation pulse.We considerS∈N SAR constraints. Each SAR constraint s∈{1,…,S} is characterized by a positive definite Hermitian matrix A(s)∈H++. An RF waveform w satisfies SAR constraint s∈{1,…,S} if(2)∫0Twt∗Aswtdt≤1.Here we have scaled the matrix A(s) in such a way that the upper bound has the value one. An RF waveform w will be calledSAR-compliant if it satisfies (2) for every s∈{1,…,S}. 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 waveformw are scaled by a factor α>0 (i.e., w(t) is replaced by αw(αt) and T by T/α) 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 matrixH∈H++ we define an associated C-dimensional, complex, centered ellipsoid E(H)⊂CC by (3)EH≔v∈CC:v∗Hv≤1.For two matrices G,H∈H++ the ellipsoid E(G) is contained in E(H) if and only if G-H is positive semidefinite, i.e.,(4)EG⊂EH⇔G-H∈H+.For later use we also note the equivalence(5)G-H∈H++⇔H-1-G-1∈H++.If an RF waveformw and a matrix X∈H++ satisfy(6)∫0Twt∗Xwtdt≤1,and E(X)⊂E(A(s)) for an s∈{1,…,S}, then w satisfies SAR constraint s. This follows from equivalence (4). In particular, if the RF waveform w and the matrix X∈H++ satisfy (6) and(7)EX⊂⋂s=1SEAs,then w is SAR-compliant. Condition (7) says that the ellipsoid associated with X must be contained in each ellipsoid which is associated with one of the SAR matrices A(s) for s∈{1,…,S}.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 ⋂s=1SE(A(s)) of the SAR ellipsoids is such a set. Since it is symmetric, its MVIE is also symmetric and can be represented as E(X) for a matrix X∈H++. The main idea of this work is to calculate this MVIE SAR matrix X and use it to replace the S 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 matrixH∈H++ the volume of the ellipsoid E(H) is given by the expression (compare, e.g., [19] equation (12.5)) (8)VH≔VCdetH-1,where VC is the volume of the C-dimensional complex unit ball.Finding the maximum volume ellipsoid can be casted as a convex optimization problem: Since the negative of the concatenated functionslog and det is convex on the convex set H++ [18], the following mathematical program is a convex optimization problem:(9)maximizelog⁡det⁡Y(10)subjecttoY∈H++,(11)As-1-Y∈H+fors∈1,…,S.According to Theorem A.1 in the Appendix below, this problem possesses a unique solution which we denote YMVIE∈H++. Its inverse XMVIE≔YMVIE-1∈H++ has property (7). This follows from condition (11) and the equivalence (5). Moreover the ellipsoid E(XMVIE) of this matrix has the maximum volume among all such matrices, since the objective function (9) is a strictly increasing function of the volume V(Y-1) of the ellipsoid which is associated with Y-1.It is interesting to note that when the maximum volume inscribed ellipsoidE(XMVIE) is blown up by the factor C then it encompasses the intersection of all original SAR ellipsoids,(12)CEXMVIE=EXMVIEC⊃⋂s=1SEAs.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 XMVIE/C and comparing it to the one with XMVIE, 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.logdet-Function for Positive Definite Matrices

We unify the objective function (9) and condition (10) into the concave function logdet:H→R∪{-∞} by setting (13)logdetH≔log⁡det⁡HifH∈H++,-∞otherwise.Our numerical implementation of the calculation of logdet(H) firstly attempts a Cholesky decomposition of the matrix H. When this decomposition fails the matrix H is not positive definite and the value of logdet(H) is -∞. Otherwise the Cholesky decomposition yields a lower triangular matrix L∈CC×C with has strictly positive real diagonal elements and satisfies LL∗=H. Thus the value of logdet(H) can be calculated by (14)logdetH=log⁡det⁡Ldet⁡L∗=2∑i=1Clog⁡Li,i.The lower triangular matrix L can be inverted by forward substitutions. Thus the inverse of H can be obtained as H-1=L-1(L-1)∗. This inverse is important for the calculation of the first- and second-order derivatives of logdet [18, 20, 21].The calculation oflogdet(H) for a matrix H∈H++ requires less than 4/3C3 multiply, add, or subtract operations, plus C square roots, inversions, and logarithms. The calculation of the C2 first-order derivatives of logdet(H) requires less than 8/3C3 additional multiply, add, or subtract operations and C inversions. The calculation of the C2(C2+1)/2 second-order derivatives of logdet(H) takes about 4C4 additional multiply, add, or subtract operations.We note that the numerical calculation of the inversesA(s)-1 for s∈{1,…,S} 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 constraintsH∈H+ or H∈H++ as in (10) and (11) is to replace them by the stronger convex constraint (15)logdetH≥Clog⁡ϵ,where ϵ>0 is a small number. The logdet-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 withindividual SAR constraint approximations(16)maximizelogdetY(17)subjecttologdetBs-Y≥Clog⁡ϵfors∈1,…,S.If B(s)∈H++ for every s∈{1,…,S}, this problem has again a unique solution. The inclusion,(18)H∈H:logdetH≥Clog⁡ϵ⊂H++,shows that by setting(19)Bs≔As-1fors∈1,…,Sthis problem is a strengthening of problem (9). Its solution is denoted by Yϵ and we set Xϵ≔Yϵ-1.The converse inclusion,(20)H+⊂H∈H:logdetH+ϵI≥Clog⁡ϵ,shows that by setting (21)Bs≔As-1+ϵIfors∈1,…,Sin (17) we obtain a relaxation of problem (9). Its solution is denoted by Yϵ′. Thus we get (22)logdetYϵ′≥logdetYMVIE≥logdetYϵ,which implies with Xϵ′≔(Yϵ′)-1(23)logdetXϵ′≤logdetXMVIE≤logdetXϵ,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 XMVIE to Xϵ and adjust the parameter ϵ. We found(24)ϵ≔1100mins=1StraceBsto be a suitable value in the examples of Section 4. As usual trace· 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 S constraints. Thus, an augmented linear system with C2+S variables must be solved in each iteration. This system contains first- and second-order derivatives at the current iterate. Since S can be much larger than C2 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 theS SAR constraints (11) by asingle SAR constraint approximation:(25)maximizelogdetY(26)subjectto∑s=1SlogdetBs-Y≥SClog⁡ϵ.This is a convex optimization problem with C2 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 B(s)-Y remains positive definite. Clearly, condition (17) implies (26) which implies (11) for B(s)≔A(s)-1. If Yϵ′′ is the solution of this problem we therefore have (27)logdetYMVIE≥logdetYϵ′′≥logdetYϵ.Again, converse inequalities can be obtained for the volumes of the associated ellipsoids.The second variant is what we call anexternal active constraint selection. Here we fix a number S′ much smaller than S, but in the order of the number of variables C2. For instance, in the examples we use (28)S′=min⁡4C2,S.The nonlinear solver is told that there are only S′≤S constraints. For a given iterate Y∈H++ during the optimization process, we calculate the values logdetB(s)-Y for every s∈{1,…,S}. The constraints are sorted by these values and the solver gets values and derivatives of the first S′ constraints in ascending order.TheS′ 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 S original constraints.Since, after the evaluation and sorting steps,S-S′ 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 C2+S′ rather than C2+S 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.


## 3.1.logdet-Function for Positive Definite Matrices

We unify the objective function (9) and condition (10) into the concave function logdet:H→R∪{-∞} by setting (13)logdetH≔log⁡det⁡HifH∈H++,-∞otherwise.Our numerical implementation of the calculation of logdet(H) firstly attempts a Cholesky decomposition of the matrix H. When this decomposition fails the matrix H is not positive definite and the value of logdet(H) is -∞. Otherwise the Cholesky decomposition yields a lower triangular matrix L∈CC×C with has strictly positive real diagonal elements and satisfies LL∗=H. Thus the value of logdet(H) can be calculated by (14)logdetH=log⁡det⁡Ldet⁡L∗=2∑i=1Clog⁡Li,i.The lower triangular matrix L can be inverted by forward substitutions. Thus the inverse of H can be obtained as H-1=L-1(L-1)∗. This inverse is important for the calculation of the first- and second-order derivatives of logdet [18, 20, 21].The calculation oflogdet(H) for a matrix H∈H++ requires less than 4/3C3 multiply, add, or subtract operations, plus C square roots, inversions, and logarithms. The calculation of the C2 first-order derivatives of logdet(H) requires less than 8/3C3 additional multiply, add, or subtract operations and C inversions. The calculation of the C2(C2+1)/2 second-order derivatives of logdet(H) takes about 4C4 additional multiply, add, or subtract operations.We note that the numerical calculation of the inversesA(s)-1 for s∈{1,…,S} 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 constraintsH∈H+ or H∈H++ as in (10) and (11) is to replace them by the stronger convex constraint (15)logdetH≥Clog⁡ϵ,where ϵ>0 is a small number. The logdet-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 withindividual SAR constraint approximations(16)maximizelogdetY(17)subjecttologdetBs-Y≥Clog⁡ϵfors∈1,…,S.If B(s)∈H++ for every s∈{1,…,S}, this problem has again a unique solution. The inclusion,(18)H∈H:logdetH≥Clog⁡ϵ⊂H++,shows that by setting(19)Bs≔As-1fors∈1,…,Sthis problem is a strengthening of problem (9). Its solution is denoted by Yϵ and we set Xϵ≔Yϵ-1.The converse inclusion,(20)H+⊂H∈H:logdetH+ϵI≥Clog⁡ϵ,shows that by setting (21)Bs≔As-1+ϵIfors∈1,…,Sin (17) we obtain a relaxation of problem (9). Its solution is denoted by Yϵ′. Thus we get (22)logdetYϵ′≥logdetYMVIE≥logdetYϵ,which implies with Xϵ′≔(Yϵ′)-1(23)logdetXϵ′≤logdetXMVIE≤logdetXϵ,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 XMVIE to Xϵ and adjust the parameter ϵ. We found(24)ϵ≔1100mins=1StraceBsto be a suitable value in the examples of Section 4. As usual trace· 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 S constraints. Thus, an augmented linear system with C2+S variables must be solved in each iteration. This system contains first- and second-order derivatives at the current iterate. Since S can be much larger than C2 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 theS SAR constraints (11) by asingle SAR constraint approximation:(25)maximizelogdetY(26)subjectto∑s=1SlogdetBs-Y≥SClog⁡ϵ.This is a convex optimization problem with C2 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 B(s)-Y remains positive definite. Clearly, condition (17) implies (26) which implies (11) for B(s)≔A(s)-1. If Yϵ′′ is the solution of this problem we therefore have (27)logdetYMVIE≥logdetYϵ′′≥logdetYϵ.Again, converse inequalities can be obtained for the volumes of the associated ellipsoids.The second variant is what we call anexternal active constraint selection. Here we fix a number S′ much smaller than S, but in the order of the number of variables C2. For instance, in the examples we use (28)S′=min⁡4C2,S.The nonlinear solver is told that there are only S′≤S constraints. For a given iterate Y∈H++ during the optimization process, we calculate the values logdetB(s)-Y for every s∈{1,…,S}. The constraints are sorted by these values and the solver gets values and derivatives of the first S′ constraints in ascending order.TheS′ 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 S original constraints.Since, after the evaluation and sorting steps,S-S′ 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 C2+S′ rather than C2+S 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⋂s=1SE(A(s)) to E(XMVIE) (or its numerical approximations) we define the power loss fv(X) in direction v∈{u∈CC,u∗u=1} for X∈H++ as (29)fvX≔mins=1Sv∗Asvv∗Xv.From these fractions we deduce the corresponding maximum, average, and minimum power losses by setting(30)fmaxX≔supv∈CC,v∗v=1fvX,(31)favgX≔∫v∈CC,v∗v=1fvXμdv,(32)fminX≔infv∈CC,v∗v=1fvX,where μ denotes the uniform probability distribution on the complex unit ball {u∈CC,u∗u=1}. From the properties of XMVIE we know (33)1=fminXMVIE≤favgXMVIE≤fmaxXMVIE≤C.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 f(v) fractions of a large number of directions v sampled from μ.


### 4.1. Reduction of SAR Matrices

Here we consider three example sets of SAR matrices:(1)
A Siemens internal example withC=2 RF coils and S=448480 SAR constraints.(2)
An example of Nicolas Boulant withC=8 RF coils and S=111939 SAR constraints.(3)
The example from Hoyos-Idrobo et al. [12] with C=8 RF coils and S=491 SAR constraints.In Section4.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 Section3:(1)
By solving the convex optimization problem (16) with theindividual SAR constraints approximation (17).(2)
By solving variant (25) with thesingle SAR constraint approximation (26).(3)
As the first method, but with theexternal active constraint selection for S′=4C2 active constraints as explained in Section 3.3.Each of these methods delivers an approximationX of the MVIE SAR parameters which imply the original SAR constraints. This has been checked by performing a Cholesky decomposition of A(s)-X for every s∈{1,…,S}. 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.Table 1
Numerical calculation of maximum volume inscribed ellipsoid for the three examples and the three proposed methods.
Example characteristics Example number ( 1 ) ( 2 ) ( 3 )  RF coilsC 2 8 8 SAR constraintsS 448480 111939 491  Method(1): individual SAR constraints  Number of iterations 33 100 22 Overall algorithm [s] 657.7 847.3 0.5 Symbolic factorization [s] 630.8 607.0 0.1 Update Hessian [s] 0.6 127.2 0.2 Compute direction [s] 22.2 97.0 0.1 Acceptable trial point [s] 2.4 0.04 0.0 Number of back steps 15 76 22  Method(2): single SAR constraint  Number of iterations 100 100 61 Overall algorithm [s] 27.9 129.9 0.9 Symbolic factorization [s] 0.1 0.1 0.1 Update Hessian [s] 0.7 115.5 0.5 Compute direction [s] 0.1 0.06 0.1 Acceptable trial point [s] 26.3 0.001 0.1 Number of back steps 1920 412 117  Method(3): external active constraint selection  Number of iterations 19 58 46 Overall algorithm [s] 2.9 4.4 0.6 Symbolic factorization [s] 0.1 0.1 0.1 Compute Hessian [s] 0.0 0.3 0.2 Compute direction [s] 0.0 0.2 0.1 Acceptable trial point [s] 0.0 3.6 0.1 Number of back steps 15 127 137Table 2
Power losses for the three examples.
Example number ( 1 ) ( 2 ) ( 3 ) RF coilsC 2 8 8  fmin(X) 1.00002 1.11045 1.10287 favg(X) 1.1802 2.56831 2.83011 fmax(X) 1.88863 5.18516 6.2459In 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 Table2 we present the power losses obtained by 107 uniformly distributed random samples from C-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 kT 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 C=8 RF excitation coils. Each coil must satisfy a maximum voltage restriction during each kT 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 kT 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)
Thearbitrary flip angle optimization (AFA) based on a solution of the Bloch equation without relaxation terms.(2)
Thesmall 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 Section4.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.Table 3
Optimization of AFA and STA objective functions for the excitation and inversion examples without and with proposed MVIE SAR reduction.
MVIE SAR reduction no yes  Number of SAR constraints 491 1 Average power constraints 8 8  30-degree excitation with 5 k T points  Number of RF constraints 40 40  Arbitrary flip angle optimization (AFA)  Computation time [s] 2.63 1.24 Iteration count 41 38 Root mean square error [deg] 1.63 1.87  Small tip angle approximation (STA)  Computation time [s] 1.87 0.89 Iteration count 38 41 Root mean square error [deg] 1.64 1.90  180-degree inversion with 7 k T points  Number of RF constraints 56 56  Arbitrary flip angle optimization (AFA)  Computation time [s] 8.52 2.87 Iteration count 84 55 Root mean square error [deg] 15.44 15.44  Small tip angle approximation (STA)  Computation time [s] 4.84 1.44 Iteration count 65 52 Root mean square error [deg] 10.38 10.38The 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 5kT 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.


## 4.1. Reduction of SAR Matrices

Here we consider three example sets of SAR matrices:(1)
A Siemens internal example withC=2 RF coils and S=448480 SAR constraints.(2)
An example of Nicolas Boulant withC=8 RF coils and S=111939 SAR constraints.(3)
The example from Hoyos-Idrobo et al. [12] with C=8 RF coils and S=491 SAR constraints.In Section4.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 Section3:(1)
By solving the convex optimization problem (16) with theindividual SAR constraints approximation (17).(2)
By solving variant (25) with thesingle SAR constraint approximation (26).(3)
As the first method, but with theexternal active constraint selection for S′=4C2 active constraints as explained in Section 3.3.Each of these methods delivers an approximationX of the MVIE SAR parameters which imply the original SAR constraints. This has been checked by performing a Cholesky decomposition of A(s)-X for every s∈{1,…,S}. 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.Table 1
Numerical calculation of maximum volume inscribed ellipsoid for the three examples and the three proposed methods.
Example characteristics Example number ( 1 ) ( 2 ) ( 3 )  RF coilsC 2 8 8 SAR constraintsS 448480 111939 491  Method(1): individual SAR constraints  Number of iterations 33 100 22 Overall algorithm [s] 657.7 847.3 0.5 Symbolic factorization [s] 630.8 607.0 0.1 Update Hessian [s] 0.6 127.2 0.2 Compute direction [s] 22.2 97.0 0.1 Acceptable trial point [s] 2.4 0.04 0.0 Number of back steps 15 76 22  Method(2): single SAR constraint  Number of iterations 100 100 61 Overall algorithm [s] 27.9 129.9 0.9 Symbolic factorization [s] 0.1 0.1 0.1 Update Hessian [s] 0.7 115.5 0.5 Compute direction [s] 0.1 0.06 0.1 Acceptable trial point [s] 26.3 0.001 0.1 Number of back steps 1920 412 117  Method(3): external active constraint selection  Number of iterations 19 58 46 Overall algorithm [s] 2.9 4.4 0.6 Symbolic factorization [s] 0.1 0.1 0.1 Compute Hessian [s] 0.0 0.3 0.2 Compute direction [s] 0.0 0.2 0.1 Acceptable trial point [s] 0.0 3.6 0.1 Number of back steps 15 127 137Table 2
Power losses for the three examples.
Example number ( 1 ) ( 2 ) ( 3 ) RF coilsC 2 8 8  fmin(X) 1.00002 1.11045 1.10287 favg(X) 1.1802 2.56831 2.83011 fmax(X) 1.88863 5.18516 6.2459In 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 Table2 we present the power losses obtained by 107 uniformly distributed random samples from C-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 kT 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 C=8 RF excitation coils. Each coil must satisfy a maximum voltage restriction during each kT 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 kT 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)
Thearbitrary flip angle optimization (AFA) based on a solution of the Bloch equation without relaxation terms.(2)
Thesmall 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 Section4.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.Table 3
Optimization of AFA and STA objective functions for the excitation and inversion examples without and with proposed MVIE SAR reduction.
MVIE SAR reduction no yes  Number of SAR constraints 491 1 Average power constraints 8 8  30-degree excitation with 5 k T points  Number of RF constraints 40 40  Arbitrary flip angle optimization (AFA)  Computation time [s] 2.63 1.24 Iteration count 41 38 Root mean square error [deg] 1.63 1.87  Small tip angle approximation (STA)  Computation time [s] 1.87 0.89 Iteration count 38 41 Root mean square error [deg] 1.64 1.90  180-degree inversion with 7 k T points  Number of RF constraints 56 56  Arbitrary flip angle optimization (AFA)  Computation time [s] 8.52 2.87 Iteration count 84 55 Root mean square error [deg] 15.44 15.44  Small tip angle approximation (STA)  Computation time [s] 4.84 1.44 Iteration count 65 52 Root mean square error [deg] 10.38 10.38The 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 5kT 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. TheC×C parameters of the MVIE SAR constraint can be obtained by solving a convex optimization problem. Taking advantage of the properties of the logdet-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.



---

*Source: 7369845-2019-06-03.xml*