Efficient performance analysis and optimization of transient‐state sequences for multiparametric magnetic resonance imaging

Abstract In transient‐state multiparametric MRI sequences such as Magnetic Resonance Spin TomogrAphy in Time‐domain (MR‐STAT), MR fingerprinting, or hybrid‐state imaging, the flip angle pattern of the RF excitation varies over the sequence. This gives considerable freedom to choose an optimal pattern of flip angles. For pragmatic reasons, most optimization methodologies choose for a single‐voxel approach (i.e., without taking the spatial encoding scheme into account). Particularly in MR‐STAT, the context of spatial encoding is important. In the current study, we present a methodology, called BLock Analysis of a K‐space‐domain Jacobian (BLAKJac), which is sufficiently fast to optimize a sequence in the context of a predetermined phase‐encoding pattern. Based on MR‐STAT acquisitions and reconstructions, we show that sequences optimized using BLAKJac are more reliable in terms of actually achieved precision than conventional single‐voxel–optimized sequences. In addition, BLAKJac provides analytical tools that give insights into the performance of the sequence in a very limited computation time. Our experiments are based on MR‐STAT, but the theory is equally valid for other transient‐state multiparametric methods.

Examples of transient-state sequences are MR fingerprinting (MRF), 12,13 Magnetic Resonance Spin TomogrAphy in Time-domain (MR-STAT), 14,15 or transient/hybrid-state imaging. Transient-state sequences have the promise of more efficiently deducing several tissue properties from one single scan. In order to achieve that, these techniques abandon the steady-state approach and introduce time-varying sequence parameters (e.g., the radio-frequency [RF] flip angle). These techniques trade simplicity for efficiency, where "efficiency" is seen as obtaining more information (low noise or artifacts, high resolution) per unit timeor, conversely, faster acquisition without the loss of accuracy or precision.
When focusing on the inherent performance potential of a sequence, rather than on the deficiency of the reconstruction process, the Cramér-Rao lower bound (CRLB) is the most appropriate tool. The CRLB has been used for single-voxel optimization. 16,17,20,26 Similarly, Byanju et al. 27 focus exclusively on low spatial frequencies, bringing it close to single-voxel, while Leitão et al. 28 assume full sampling at every instant, also turning the optimization into a single-voxel approach. Other work [22][23][24] analyzes the efficiency in the context of the encoding pattern and focuses on the undersampling artefacts obtained by the specific reconstruction process; that work does not apply the CRLB. Also, Jordan et al. 30 take encoding into account but choose not to apply CRLB as a criterion. Up-to-date application of the CRLB in the context of a RF time-varying train combined with gradient encoding has not been presented before.
In our approach, we take the spatial encoding-pattern into accountand we focus on the theoretical CRLB, which we study by means of the Jacobian and its associated Fisher Information Matrix. This is a large matrix, of the order of n Â n; in non-Cartesian sequences, n is the number of samples in k-space; in Cartesian sequences, it is the number of different phase-encoding values (in our MR-STAT examples, n ¼ 224Þ. A full inversion of a matrix of that size incurs a prohibitive processing burden. In this paper, we present a solution that is both fast and useful. As a solution, we propose a methodology that factorizes the Jacobian, a methodology called BLock Analysis of a K-space-domain Jacobian (BLAKJac) (see Section 2). BLAKJac is fast relative to a full inversion because it achieves an n 2 factor in performance increase in the analysis of sequences. This allows for iterative sequence optimization in the context of encoding with a total computing time of the order of minutes. As one useful feature of BLAKJac, we show, on a Cartesian-encoded sequence reconstructed by MR-STAT, 31 that sequences optimized in the context of spatial encoding (2D) reliably predict the precision of T 1 -or T 2 -maps, as opposed to single-voxel (0D) optimized sequences. As another useful feature, BLAKJac allows for prospective analysis of spatial noise spectra in the obtainable tissue parameter maps. Given any sequence of time-varying flip angles and gradient-encoding patterns, it can predict and visualize which spatial frequencies will be most affected by noise, as shown in the simulations and experimental tests.

| THEORY
To evaluate the performance of a sequence, it is necessary to evaluate the inverse of the Fisher Information Matrix J H J, where J is the Jacobian of the associated signal model. In order to allow efficient inversion of the information matrix, we propose a factorization of the Jacobian based on the following idea (see also Figure 1): we focus on the spatial frequency domain and consider thus the Fourier transforms of properly scaled parameter maps (see Equation 2). These transformed maps are denoted as "auxiliary variables" with the symbol "f" in the sequel. Note that we do not need to actually calculate f as an intermediate step: this notion only serves to explain the factorization of the Jacobian.
Thus, we consider the composite signal model Þfor a given function h, where h ∘ f ¼ Γ. This leads to a splitting of the Jacobian into two factors: J Γ=u ¼ J Γ=f J f=u , where each of these two factors leads to a computationally efficient inversion of the information matrix.
F I G U R E 1 Proposed signal model decomposition and associated Jacobian factors. The arrows represent mappings, not processing steps. The introduction of the auxiliary variables f allows to factorize the Jacobian into two factors, J f=u and J Γ=f ; both can be efficiently evaluated (Here, Γ is the signal model, see also Equation (A1) in Appendix A.1; u is the set of unknowns [e.g., the T 1 and

| CRLB for spatially dependent parameter maps
Transient-state multiparametric MRI sequences are often used to estimate tissue properties like T 1 , T 2 , or proton density. [12][13][14][15][16]32 For reconstruction techniques that are unbiased, there is a theoretical limit on the uncertainty to the resulting parameter maps, namely, the CRLB: Here, "E½ " means "expectation value", u denotes the vector of unknowns, and b u is an unbiased estimator thereof. "H" denotes the Hermitian conjugate and I is the Fisher Information Matrix.
To evaluate the Fisher Information Matrix, we assume that we have a perfect forward MR signal model Γ : ð Þ that exactly models spin dynamics and sequence-design settings, leading to sampled data. We also assume that this data is only confounded by independent and identically distributed (i.i.d.) thermal noise with variance σ 2 . In that case, I u ð Þ reduces to 1 σ 2 J H Γ=u J Γ=u . Here, the notation J Γ=u is the Jacobian matrix ∂Γ u ð Þ ∂u: Note the similarity in formalism to, for example, the work by Zhao et al. 17 ; the fundamental difference is that Zhao et al. 17 consider the vector of unknowns to consist of just N p elements, with N p the number of properties (e.g., three), thereby defining u ¼ T 1 , T 2 , M 0 ½ T (with M 0 representing the proton density). In our work, u is a vector of N p N locations elements, which thus represents the whole spatial distribution of each parameter type.
In the sequel, the symbol r will denote a three-elements location vector (of which there are N locations [i.e., the number of voxels]); because we assume that the set of locations spans a Cartesian grid, we apply a slight notational abuse and use r also as an index to u (i.e., the full vector u consists of elements u p,r ). Similarly, the k-space vector k ¼ k x , k y , k z ð Þwill also be used as an index having N k distinct values. Here, we take For any flip-angle and gradient-encoding scheme, we can calculate J Γ=u , from which the expected noise statistics in the reconstructed parameters can be calculated by the inversion of J H Γ=u J Γ=u . Yet, in practice, this inversion is problematic. This Fisher Information Matrix has size N p N k Â N p N k ; the inversion of a full matrix of this size requires O N 3 k N 3 p operations. A factorization of the Jacobian, elaborated below, allows for a substantial speedup of the aforementioned inversion problem.

| Taylor expansion of the model Γ and factorization of the Jacobian
Appendix A elaborates on a Taylor expansion of the signal model Γ : ð Þ, which allows for a factorization of the Jacobian. That formalism is based on a first-order Taylor approximation of the Bloch model, that is, on The aforementioned factorization separates the (weighted) Fourier-transform component from the Jacobian, the "remaining" factor being close to a block-diagonal matrix. More precisely, we define auxiliary variables as where "⨀" is the Hadamard product and F Á ð Þ is shorthand for n-dimensional Fourier transform (i.e., F x ð Þ k ¼ P r x r e i2πkÁr ); as already mentioned, the vector k ¼ k x , k y , k z ð Þis also used as an index, similarly to vector r. While f 0 represents the Fourier transform of the proton-density map, the physical meaning of, for example, f 1 , can be seen as the Fourier transform of a weighted T1 map.
Using these auxiliary variables, we consider the signal model as Γ ¼ h f u ð Þ ð Þ (see Figure 1). We thereby split the Jacobian into two factors: • The factor J Γ=f is approximated as a block-diagonal. This diagonal structure is shown in Figure 2.
• The factor J f=u is a weighted Fourier transform matrix, with well-known noise properties, ; here, σ p is the Cramér-Rao lower limit to the standard deviation of the reconstructed noise in property p (e.g., the T 1 -map). The proportionality factor in the formula above includes the choice of u ref p , the noise level in the data (σ), as well as the proton density.
Note that J Γ=f is close to being a block-diagonal, while the original J Γ=u is not. This is a crucial aspect because it allows for the N 2 k efficiency factor when inverting the Fisher Information Matrix. The resulting performance gain is particularly beneficial if the inversion J H Γ=f J Γ=f À1 needs to be evaluated multiple times, for instance, during an iterative optimization process.

| Analysis and applications of BLAKJac
The proposed BLAKJac factorization allows for several useful applications, among which are noise spectrum analysis and sequence optimization.

| BLAKJac for noise spectrum analysis
BLAKJac can efficiently predict the noise spectrum in the reconstructed parameter maps u. Assuming an i.i.d. noise at the input, the noise covariance matrix at the auxiliary variables f is approximated by . By examining the resulting noise variance on a given property p and a measurement sequence, we can extract an N k Â N k submatrix, Given the block-diagonality of J Γ=f , this submatrix is diagonal.
Consequently, BLAKJac provides the information on the noise variance as a function of spatial frequency, as will be shown in the numerical and experimental tests ( Figures 10 and 11, details will follow).

| BLAKJac for sequence optimization
Because the Jacobian J f=u is in essence a Fourier transform matrix, the noise variance in u p can be obtained by a sum of the variances of its k-space components. Symbolically, F I G U R E 2 Graphical representation of the absolute values of the elements of the Jacobian matrix in a Cartesian case (where k x is not relevant, so k reduces to k y ). (A) J Γ=u , (B) The matrix J Γ=f . Part (C) shows the same content as (B), but with the indices rearranged so that the instance i and the property p are the inner indices, while k and l are the outer indices. Part (C) demonstrates that the 224 diagonal blocks, k ¼ l, are dominant. In the example, the blocks have a size of 5 Â 3. The black color corresponds to 0 with σ p the expected standard deviation to the map of property p (e.g., the T 1 -map).
From the values of σ p , we can calculate an overall figure of merit C as C ¼ 1 Optimal sequences can be obtained by minimizing C. This is elaborated in section 3.2.
Our approach also allows to map the performance of a sequence over a T 1 , T 2 ð Þ -landscape, thereby highlighting potentially suboptimal encoding for specific values of T 1 , T 2 ð Þ( Figure 3).

| Details of the applied MR-STAT sequence
In all our experiments, we applied Cartesian MR-STAT measurements. 15 An MR-STAT sequence, as shown in Figure 4C, consists of a repeated acquisition (T R ¼ 10 ms) of single lines through k-space; T R as well as the total gradient area between pulses is held constant (a.k.a. steady-state free precession), while the flip angle is varied from pulse to pulse. Experiments were performed on a Philips Achieva 3T system with a 13-element head coil (but precombined into a single virtual coil for our experiments, by taking the first principal component of the concatenated coil signals).
In all experiments, we measured a single slice with 4-mm thickness. The scan parameters were T E ¼ 5 ms with a voxel size of 1 mm Â 1 mm, 2 , and no parallel imaging. A sequence is always preceded by a single adiabatic inversion pulse.

| Sequence optimization
We used BLAKJac to calculate and minimize a figure of merit C given a time-varying sequence of 1120 flip angles. According to eqn. A2 in Sbrizzi As explained in the Theory section, σ p is calculated by inverting the Fisher Information Matrix. The inversion requires O N k N 3 p operations. In general, N k is of the same order as N locations , which in our example was 224 2 ; furthermore, we applied N p ¼ 3 (proton density, T 1 , and T 2 ). If we consider Cartesian sequences and we neglect the relaxation during a readout line, then all samples of that readout line share the same flip angle history, thus they share the same J Γ=f . In that situation, the inversion of the Fisher Information Matrix requires only O N phaseÀencodings N 3 p operations, which in our example was of the order of 224 Á 3 3 . Given the compact block-diagonal representation of BLAKJac, the processing time is dominated by the evaluation of the Jacobian for all 1120 timepoints (i.e., by evaluating the entries . This evaluation can be performed efficiently using a finite-difference method based on Extended Phase Graphs. 33 In practice, the evaluation and inversion of the Fisher Information Matrix required approximately 10 ms on a single-core Intel Xeon CPU E5-1620 v3 running at 3.50 GHz. (see the red crosses in Figure 3); these are spanning the range of relaxation values of frequently occurring human tissues. The combined criterion was taken as an average over the outcomes: In principle, S refers to all the relevant sequence settings, but in the sequel we will only consider the time-varying RF flip angles as being variable in time and consider all the other parameters as being constant. Consequently, approximately 70 ms was needed to evaluate C S ð Þ given a sequence S of 1120 time points.
Referring to Equation (3), and as motivated in the Discussion section, we slightly modify the estimation of σ p by weighing the 10 lowest spatial frequencies with a higher factor (3.0) than all other spatial frequencies.
Using the Nelder-Mead algorithm 34 from the Optim.jl package of the Julia programming language, 35 we minimized C S ð Þ over the RF flip angles of S. There were no explicit constraints on the flip angles, but to account for RF power constraints, we slightly modified the criterion by introducing a peak-power-correction factor, obtaining C overall S ð Þ¼C S ð Þ= . Frame colors correspond to colors used in subsequent graphs; "best" and "worst" refer to the performance in a 2D setting. Part (C) is, enlarged, the first sequence of (A), including a representation of the phase-encoding pattern. A sweep of 224 phase-encoding values is repeated five times, so the acquisition time of the whole sequence is 224 Á 5 ð ÞÁ10 ms ¼ 11:2 s Scope Crafts et al. 36 To reduce the sensitivity to local minima, we applied a multilevel approach: first optimizing on 10 equidistant points, using the result as a basis for 15 points and the result thereof for a 20-point optimization. In each of these three steps, the number of Nelder-Mead iteration steps was limited to 2000.
The optimization process was followed by a selection process because the function C overall S ð Þ has multiple local minima. By initiating the optimization with nine slightly different initial states, different local optima S opt,2D "2D" indicates that the optimization was performed in the context of phase encoding (see Figure 4C). The nine resulting shapes in Figure 4A seem mutually very different, indicating that not all of them can be the true optimum. Because we are dealing with a nonconvex minimization problem, the true optimum would require an exhaustive search, which would incur a prohibitive computational burden. But because the nine outcomes are mutually very close, this increases the confidence that all of them are close to the true optimum.
In a similar fashion, we also optimized on S, assuming an absence of any encoding, resulting in sequences S opt,0D i . Note that this 0D optimization approach is the most frequently adopted in the existing literature. [16][17][18][19][20][21]26,29 To simplify the experimental validation, the best and the worst of each set was chosen based on the C overall S ð Þ criterion. The flip angles of the resulting four sequences are called "2D-best", "2D-worst", "0D-best", and "0D-worst", where "best" and "worst" always refer to the performance in a 2D setting (i.e., including phase encoding) ( Figure 4A,B).

| MR measurements
Our measurements were intended to validate and apply the BLAKJac framework. In particular, we investigated whether: 1. Experimentally obtained noise levels and spatial-frequency noise spectra correspond to those predicted by BLAKJac.
3. Analytical tools (noise spectrum diagrams) derived from the BLAKJac framework successfully predict the spatial quantitative and qualitative behavior of a measured sequence. We performed MR-STAT measurements with the four sequences derived in subsection 3. , which is approximately 4%.
• A 2D axial slice of a phantom setup housing 16 vials of the Eurospin phantom.
• A 2D axial slice of the brain of two healthy volunteers (with approved consent according to the guidelines of the ethics committee).
The quantitative T 1 and T 2 maps for each experimental setup were obtained applying the MR-STAT reconstruction. 14,38,39 4 | RESULTS

| Optimized sequences
The result of the sequence optimization is shown in Figure 4. As explained in section 3.2, two sets of RF sequences were obtained, called S opt,2D i and S opt,0D i , with i ¼ 1,…,9. We observe that these must be local optima because the evaluated performance within a set is not identical. Within the 2D-set S opt,2D i , the ratio between the highest and lowest calculated noise level is 10%, that is,

| Single-voxel ("0-D") measurements
The measured standard deviations for the 0-D phantom reconstructions are reported in Table 1.
Within the expected precision, the measured ratios between the sequence performances were in accordance with those predicted by BLAKJac: for a 0D-measurement (i.e., a measurement without phase encoding), the 0D-optimized results were marginally better than the 2Doptimized results. Table 1 also shows that the difference between 0D-best and 0D-worst is marginal. This is not surprising because these sequences were all optimized in a 0D-setting.

| Phantom 2D measurements
The 2D reconstruction results of the Eurospin phantom are shown in Figure 5. A high level of artefacts in the 0D-worst sequence is clearly visible.
The standard deviations are displayed in Figure 6. The 0D-worst sequence resulted in noise levels that are approximately twice as high as those of the other sequences. This shows that there is a high performance spread among the 0D-optimized sequences, although they all performed almost equally well in a 0D setting. That spread is negligible among the 2D-optimized sequences.
From Figure 7, it is apparent that the 0D-optimized sequences, 0D-worst in particular, show substantially higher quantitative errors than the 2D-optimized sequences. These quantitative errors are attributed to poor rendition at the low spatial frequencies, as will be further elaborated in the Discussion section. Note: the measured standard deviation is compared with the BLAKJac estimate (the scaling factor was chosen to maximally match the measured T1 results). The obtained mean (T1, T2) values are reported in the last column. The main message is that the 0D-optimized sequences perform slightly better when applied on 0D than 2D-optimized sequences do when applied on 0D. This can be seen by the deviations in the last two rows being slightly lower compared with the "2D" rows. Abbreviation: BLAKJac, BLock Analysis of a K-space-domain Jacobian.

T A B L E 1 Single-voxel ("0D") measurement (Measured) results
F I G U R E 5 Reconstructed maps (left T1, right T2) of the phantom setup with the Eurospin vials. It can be seen that the 2D-optimized results have consistent image quality, while the 0D-optimized results do notalthough the 0D-optimized sequences score mutually very similarly in a 0D-setting

| In vivo 2D measurements
In vivo results for the four types of scan are shown in Figure 8. It is clear that the 0D-worst results are inferior to the others in several aspects. This is also reflected in the measured standard deviations (Figure 9). In addition, the 0D-sequences clearly create biased maps (with substantial differences between "best" and "worst"). This is compatible with Figure 8 (phantom bias) for T 1 ≈ 1000 and T 2 ≈ 70. Section 5 (Discussion) elaborates on a likely reason for the 0D-sequences to show more bias.
Each scan was repeated 10 times, providing a pixel-wise estimate of the standard deviation with a relative precision of 1= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 Â N repetitions À 2 p ≈ 0:24. To further increase the accuracy of the estimate, we averaged the standard deviation (a) over an ROI of 20 x 20 pixels of presumed pure white matter, and (b) over the whole image.

| Noise spectra analysis
As explained in section 2.3.1, BLAKJac can predict the spatial spectra of the resulting noise. The predicted spectra for all nine sequences are plotted in Figure 10. . The cross-hatched bars represent the predictions in the "2D" approach, that is, analyzed in the context of phase encoding, while the bulleted bars represent predictions according to the single-voxel ("0D") optimization approach. It can be seen that the 0D-approach (bulleted bars) predicts very similar noise performance for all sequences (even indicating that "0D worst" would be slightly better than the others); to the contrary, a 2D-analysis in BLAKJac (cross-hatched) is able to indicate that some sequences perform significantly worse in the actual scenario (i.e., with spatial encoding). These predictions match substantially better with the measured noise levels (solid bars). Consequently, the 2D-optimized sequences display consistent noise level, while the 0D-optimized sequences do not. BLAKJac, BLock Analysis of a K-space-domain Jacobian F I G U R E 7 Correspondence between reconstructed mean values and gold standard for the Eurospin vials (bullets). Each bullet represents a gel-filled vial. The dashed lines are the linear regression lines through the bullets. The full blue diagonal is the identity line. The 0D-worst optimization is clearly inferior, particularly on T2, where the relation to the gold standard value is lost or even slightly inverted. Also here, the message is that the 2D-optimized results are mutually consistent, while the 0D-optimized results are not depend on the selected T ref . The resulting figure shows a much higher consistency among the nine spectra of the 2D-optimized sequences than among the 0D-optimized. Most notably, the 0D-optimized sequences tend to have peaks on unexpected spectral locations; some of them, particularly the one denoted as "worst" (in gray), happen to have their peak at the low spatial frequencies (jk y j ≈ 0).  Figure 11 shows the BLAKJac-generated spectra (dashed lines) alongside the measured spectra for the phantom. Most spectral features from the BLAKJac prediction are clearly recognizable in the measured noise spectra.

| DISCUSSION
We present a methodology, called BLAKJac, which can evaluate the precision performance of the combined time-varying RF train and gradientencoding scheme in only a few milliseconds, allowing for flip-angle optimizations in the context of spatial encoding.
To show the added value of the BLAKJac analysis, we compared RF-flip-angle sequences that, on one hand, have been single-voxeloptimized (i.e., ignoring the spatial encoding, or "0D"), and, on the other hand, sequences that have been optimized in the context of phase encoding ("2D"). We have shown that a set of 0D-optimized sequences, which all lead to mutually very similar noise levels when analyzed in a 0D-setting, do lead to very unpredictable performance in the relevant scenario, that is, when spatial encoding is active (2D). On the other hand, the set of 2D-optimized sequences based on the BLAKJac criterion are mutually of very similar performance and, in general, outperform the 0Doptimized sequences. This is shown in phantoms as well as in vivo measurements (Figures 6 and 9, respectively); the actually achieved noise level  Figure 4. The dispersion of the curves in (C) and (D) shows that a single-voxel-optimized sequence has an unpredictable outcome when applied in a setup with phase encoding, while the outcome is much more predictable when applying 2D-optimization (A and B). Evaluation was performed assuming a seven-tissues mix (see section 3.2). BLAKJac, BLock Analysis of a K-space-domain Jacobian of the 0D-optimized sequences can be substantially higher than the noise level in 2D-optimized sequences. These findings are a natural consequence of the BLAKJac analysis: if the time-varying flip angle is evaluated in conjunction with the underlying gradient-encoding scheme (2D-optimization), the spatial noise spectra exhibit a much higher level of control ( Figure 10A,B). The performance of 0D-optimized sequences strongly depends on the initialization condition of the optimization process: sequences that are almost equally optimal in the single-voxel sense (e.g., the 0D-best and 0D-worst sequence) lead to widely varying outcomes in actual phase-encoded measurements.
Another feature of BLAKJac is its ability to calculate expected noise levels and spatial noise spectra in tissue-parameter maps. We show that spatial noise spectra reconstructed from measurements match the BLAKJac-predicted spectra ( Figure 11). The noise spectra analysis is helpful to gain insights into the 0D optimization process: the 0D-optimized sequences display a characteristically peaked noise spectrum ( Figure 10C,D).
This can be particularly unfavorable if a peak emerges in the low spatial frequencies (around jk y j ¼ 0, or around jkj ¼ 0): these spatial frequencies are not only affected by noise, but model errors-which will be discussed extensively below-also predominantly affect the low spatial frequencies. As a consequence, model errors are multiplied by the BLAKJac-calculated factor J H Γ=f J Γ=f À1 at low spatial frequencies; this factor happens to be particularly high for the 0D-worst sequence (as in the gray curves of Figure 10C,D). This, in turn, causes a significant error in reconstructed T 1 and T 2 values in 0D-optimized sequences. In contrast to 0D-optimized sequences, the spectra of 2D-optimized sequences do not display peaks.
While Figures 6 and 9 show that the noise standard deviation is highest with the 0D-worst sequence, Figures 5, 7, and 8 also show substantially increased quantitative errors for, in particular, the 0D-worst sequence. This is most clearly seen in the dashed gray line in the right half of Figure 7, indicating that the reconstructed T 2 of that sequence has a negative correlation with the true T 2 .
In essence, our approach focuses on precision. As in all inversion problems, accuracy will be affected by differences (i.e., "model errors") between the true underlying physics and the applied model. These model errors can have various causes, including erroneous knowledge of the B þ 1 field, erroneous slice profile, erroneously assuming a single-species model in each voxel, 40 neglecting magnetization transfer, 41 and neglecting diffusion weighting induced by the sequence, 42 etc. We do have some indications that the accuracy, for example, the extent to which diffusion influences the T 2 estimate, does depend on the sequence. This might explain the difference in bias in T 2 between the 2D-worst, 2D-best, and 0D-best sequences, which, according to BLAKJac, result in very similar performance (i.e., very similar performance in terms of precision).
Although it is beyond the scope of this paper, it is certainly indicated to study the relation between sequence and accuracy (e.g., T 2 -dependence on B þ 1 or diffusion), and to subsequently minimize these dependences. These dependences are a factor in the inaccuracy (or bias). In addition to these dependences, there is another factor influencing the bias that can be directly tackled by BLAKJac, namely, the value of for low spatial frequencies (jkj $ 0Þ. As can be seen from the gray curves of Figure 10C,D, this value is particularly high for the 0D-worst sequence. Empirically, for the 2D-optimization, we found it beneficial to penalize the matrices with a weighting factor of 3.0 for k ¼ 0, gradually decaying to 1.0 at jkj ¼ 20. This operation is irrelevant for the 0D-optimization: in the absence of encoding, each sample is considered to be at k ¼ 0. Particularly when optimizing a multiplicity of parameters-in our case, 20-optimization might lead to false minima. The nine resulting shapes in Figure 4A seem mutually very different, indicating that not all of them can be the true optimum. Because we are dealing with a nonconvex minimization problem, the true optimum would require an exhaustive search, which would incur a prohibitive computational burden. But because the nine outcomes are mutually very close, this increases the confidence that all of them are close to the true optimum. Likewise, the optimization outcomes of the 0D-sequences in Figure 4B are also mutually very similar-if evaluated in a 0D-context-but they differ substantially when evaluated in a 2D-context. That variation is not caused by finding false optima, but is the result of optimizing in a 0D context and applying them to the relevant scenario (i.e., 2D).
In principle, the BLAKJac methodology applies equally to non-Cartesian sequences, like those frequently used in MRF, which typically adopt spiral or radial encoding strategies. Yet, BLAKJac does not necessarily apply to any type of MRF reconstruction. The basic MRF reconstruction approach via a dictionary match of strongly undersampled images 13 is characterized by an additional perturbation term, which consists of aliasing artefacts. These kind of perturbations are not captured by BLAKJac, which studies the CRLB, as derived from purely stochastic (thermal) noise.
We hypothesize that BLAKJac should be applicable to MRF sequences reconstructed using advanced inverse problem techniques, 43,44 which better approach the ideal result expressed by the CRLB.
Although the BLAKJac methodology is applicable to a broad category of encoding patterns, in this work we focused on a Cartesian approach using MR-STAT. We chose to not take coil sensitivity into account, and thus we do not elaborate on parallel imaging. Intuitively, if the RF-flipangle patterns are smooth with respect to the phase-encoding index (as in our MR-STAT examples), one can regard undersampling-pattern optimizations 27 as being independent of the optimizations of the flip-angle pattern and therefore they can be treated separately. Parallel imaging with generalized encoding schemes (e.g., spiral and golden angle radial) can be included in principle, although the resulting formalism would substantially increase in complexity, and this would go beyond the scope of this paper.
Further potential applications of BLAKJac include sequences of RF flip angles with optimized phase variation, optimization of 3D sequences, and potential extension to other parameters such as B1+. Constraints such as maximum specific absorption rate can be easily incorporated into BLAKJac. Also, BLAKJac may be applied to optimize the encoding pattern given a flip angle sequence, or to optimize encoding and flip angle jointly. Finally, the clinical application of BLAKJac-optimized sequences will be the subject of future research.

| CONCLUSION
The proposed BLAKJac analysis is an efficient methodology for evaluating the precision performance of transient-state multiparametric MRI sequences (like MR-STAT) in the context of applied gradient encoding. The noise levels and spectra predicted by BLAKJac correspond to those measured in vivo and in phantoms. Because the calculations take only few milliseconds, BLAKJac can optimize sequences taking into account the phase encoding. This is shown to have clear benefits in MR-STAT over single-voxel optimization. BLAKJac, which uses the 2D-encoding information, is much better at predicting the actual noise level in images compared with a single-voxel approach.

ACKNOWLEDGMENTS
This work has been financed by Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) grant number 17986, which has partly been funded by the company Philips.

Miha Fuderer
https://orcid.org/0000-0002-5673-915X which can be written as A.1.1 | Factorization of the Jacobian Using the definition of the auxiliary variables defined in Equation (2), we can write From Equation (A6), it is clear that this tells that, if a system is linear, a sample with spatial encoding k will not be sensitive to a different spatial frequency l of the maps u 0 ⨀ up We point out that there is an exception to the above conclusion if the properties u p,r can be considered to be real-valued, as it is the case for the relaxation parameters T 1 and T 2 : in this case, dependence emerges between f p,k and f p,Àk , leading to a matrix that has both a diagonal and an antidiagonal component. This is further elaborated in Appendix C.

APPENDIX B B.1 | HIGHER-ORDER TERMS
In the hypothetical case where the signal model would be sufficiently described by a first-order Taylor expansion of the function G around u ref , it has been shown that the off-diagonal elements of J Γ=f are zero. This Appendix elaborates on these off-diagonal elements in the situation that G also has substantial second-order Taylor components.
For the sake of brevity, we simplify the situation, assuming that we have just one unknown property u 1,r next to the proton density u 0,r . Any additional second-order Taylor terms would follow an equivalent analysis, as well as Taylor terms associated with cross-products, for example, The off-diagonal elements of J Γ=f have the form ∂Γ k,i ∂f p,l with k ≠ l, where k represents the encoding vector in k-space and l represents the spatial frequency of the property map, which is also a vector in the same k-space.
Because f 0 and f 1 are heavily coupled, it only makes sense to jointly calculate the partial derivative to these (i.e., to calculate Because this derivative is formed via all of the values u 0,r and u 1,r , we have Using the relation Γ k,i ¼ P r u 0,r G k, i ð Þ u 1,r , …, u NpÀ1  , evaluated at k-space position k À l. Given the knowledge that such a spatial pattern exhibits strong correlation among neighboring voxels (like in, e.g., T1-weighted MR images), the influence of that higher-order k-space term can be expected to quickly decrease when considering off-diagonal components in the J Γ=f matrix (i.e., with increasing values of jk À lj).
APPENDIX C.

C.1 | THE CASE OF REAL-VALUED PROPERTIES AND CORRESPONDING SYMMETRY IN K-SPACE
In MRI, the data d k,i (modeled by Γ k,i ) are assumed to be complex. The unknowns u are obtained by a large-scale nonlinear inversion of the set of equations Γ k,i ¼ X r u 0,r G k, i ð Þ u 1,r , …, u NpÀ1 ð Þ ,r À Á e i2πkÁr : An unconstrained solution of this inversion problem would lead to complex values of the unknowns u p,r . Yet, many properties (e.g., the relaxation time T 1 ) are known to be real quantities.
Without loss of generality, we here consider all unknowns to be real. This has an impact on the structure of the Jacobian J Γ=f ¼ ∂Γ k,i ∂f p,l h i , because we now have a constraint on f. Given the relation between u and f, as specified in Equation (2), the assumption of a real-valued u implies that f p,k ¼ f Ã p,Àk . This reduces the number of (complex) unknowns (e.g., by considering only positive values of k for the unknowns). Consequently, the number of columns in the Jacobian reduces to roughly N p N k =2 (see sequel). The data, however, may still be available for all positive and negative values of k.
This symmetry implies a slight modification of the forward model of Equation (A6), which is reformulated as Γ k,i ≈ P Np À1 p¼0 w p,k,i Áf p,k ky ≥ 0 P Np À1 p¼0 w p,k,i Áf Ã p,Àk ky < 0 8 > < > : : (Here, the symmetry plane has been taken perpendicular to k y ; in principle, any other plane through k ¼ 0 could be taken). We make use of an auxiliary variable Γ 0 and we rearrange the vector Γ 0 such that positive and negative k y -indices are adjacent.
With this reformulation, and assuming a symmetric range of sampled encodings k, the matrix J Γ 0 =f consists of n ¼ ceil N k þ 1 ð Þ =2 ð Þnonzero blocks of size 2N instances Â N p blocks having content w p,k,i w Ã p,Àk,i ! (except for the outmost k-values and for k ¼ 0).
The Fisher Information Matrix thereby has n nonzero blocks of size N p Â N p .
This modification leads to better estimates of the noise (results not shown) compared with a block-diagonal Jacobian J Γ=f where the k-space symmetry is ignored.