A SPIRED code for the reconstruction of spin distribution

In Nuclear Magnetic Resonance (NMR), it is of crucial importance to have an accurate knowledge of the sample probability distribution corresponding to inhomogeneities of the magnetic fields. An accurate identification of the sample distribution requires a set of experimental data that is sufficiently rich to extract all fundamental information. These data depend strongly on the control fields (and their number) used experimentally. In this work, we present and analyze a greedy reconstruction algorithm


Introduction
Quantum Control (QC) is nowadays a well-recognized area of research [1][2][3][4][5] with many applications ranging from magnetic resonance [6][7][8] and atomic and molecular physics [9][10][11][12] to quantum technologies [7,13,14].In this context, the control of spin dynamics by means of magnetic fields is today a well-established tool in Nuclear the control (see however [15] for a closed-loop example).A good agreement between theory and experiment is achieved if all the parameters of the model system are perfectly known within a given range of precision.In the context of NMR, the latter may be due to the spatial extension of the sample which means that in practice all its elements are not subjected to the same external electromagnetic field.The values of such parameters can be estimated experimentally but can also be actively found by using specifically adapted controls.To this aim, different approaches using quantum features have been developed recently with success [16][17][18].Among others, we can mention inversion techniques [19], selective controls [20][21][22][23], the maximization of quantum Fischer information [24][25][26][27][28] and the fingerprinting approach [29,30].Such methods allow one to estimate the value of the Hamiltonian parameter as well as its variation range.However, the latter is not the only interesting quantity and the probability distribution is also a key feature of the experimental sample.When controlling an ensemble of quantum systems, this distribution can be interpreted as the number of individual systems having a given value of the parameter.The distribution can have a simple form such as a Gaussian or a Lorentzian one.In this case, the identification is quite straightforward and can be done using standard techniques.However, the identification is much more difficult when the distribution has a complex structure with, e.g., several peaks.
In a previous work [31], we introduced a Greedy Reconstruction Algorithm (GRA) to identify in a systematic way the probability distribution of one given Hamiltonian parameter.This was based on the framework presented in [32,19].In particular, we focused on an ensemble of spin 1/2 particles in liquid state NMR subjected to an inhomogeneous radio-frequency magnetic field [33,6,34,8,35,36], where the algorithm was successfully applied to identify the distribution of the scaling factor corresponding to the magnetic field inhomogeneity.More precisely, we considered an inhomogeneous ensemble of uncoupled spin systems in liquid state whose dynamics are governed by the Bloch equation, but subjected to a different radio-frequency magnetic field.Applying a given number of constant control fields to the sample during a fixed time, the algorithm was able to determine the probability distribution of the inhomogeneous parameter from the measure of the final ensemble magnetization.We showed numerically in [31] that GRA is able to identify with a good accuracy nontrivial probability distributions with several peaks or with a step variation.An optimized version of this algorithm was also proposed to further improve the identification process and to limit the number of required controlled processes.Notice that a convergence analysis of the algorithm was only briefly sketched in [31], without rigorous proof.
The goal of the present paper is to extend the work [31] from three different points of view.Thus, the first novelty of this work is the extension of the GRA presented in [31], designed for the reconstruction a single parameter distribution, to the case of joint probability distributions of two distinct inhomogeneous Hamiltonian parameters.We consider numerically the case of both inhomogeneous static and radiofrequency magnetic fields applied to the sample.The second novelty is that we provide full MATLAB codes (SPIRED) implementing our GRA and its optimized version (denoted by OGRA) to find spin distributions.These codes can be directly used to solve the problems presented in [31] (in addition to the ones investigated in the present paper).The third novelty is that we prove theoretical convergence results for GRA.These cover also the ones only stated in [31].As a result, this paper is a selfcontained work that not only considers a more general problem than the one presented in [31], but also provides a full MATLAB code and a detailed and rigorous convergence analysis.
The paper is organized as follows.The identification problem of the spin distribution in NMR is presented in Sec. 2. The different variants of the greedy reconstruction algorithm are described in Sec. 3. Section 4 is dedicated to the description of the structure of the code SPIRED and its use.A convergence analysis of the algorithm is provided in Sec. 5. Numerical results are presented in Sec. 6.Conclusion and prospective views are drawn in Sec. 7. Additional results are presented in Appendix A.

Identification of spin distribution
The framework of our SPIRED code is illustrated in a standard control problem in liquid-state NMR, i.e., a spin ensemble subjected to both inhomogeneous static and radio-frequency magnetic fields [33,37,8].In a given rotating frame, each isochromat is characterized by a Bloch vector M = [  ,   ,   ] ⊤ , evolving in time according to the equations Notice that the components of M satisfy  2  +  2  +  2  =  2 0 , with  0 the equilibrium magnetization.Here,   and   are time-dependent controls corresponding to the components of the radio-frequency magnetic field along the and directions.The parameters  and  correspond to offset, i.e., to the inhomogeneity of the static magnetic field applied along the direction, and radio-frequency magnetic field inhomogeneities, respectively [6].In standard experiments, we have

2𝜋
are expressed in Hz.We consider a typical field amplitude  0 that can be fixed, for instance, to  0 = 2 × 100 Hz.We introduce normalized coordinates as In what follows, we omit the prime to simplify the notations.We deduce that the differential system can be expressed in normalized units as with  2 +  2 +  2 = 1.The initial state of the dynamics for each spin is the thermal equilibrium point, i.e., X 0 = [0, 0, 1] ⊤ .We consider a control time of the order of 100 ms, that corresponds to a normalized time  ′  of the order of 10.The range of variation of the parameter where Δ 0 is a frequency value that can be used to shift arbitrarily the interval.For the purpose of this paper, we assume that Δ 0 ≥ 0.4, meaning that Δ ≥ 0.
The goal of our SPIRED code is to estimate simultaneously the distributions for the parameters  and Δ by designing specific controls (  ,   ).We consider an ensemble of  spins whose dynamics are governed by Eq. ( 1).We assume that the control amplitudes (  ,   ) belong to the admissible set where   is the maximum amplitude of each component.A simple way to proceed can be described as follows.We consider that the system of  spins is divided into  Δ groups, and we associate with the -th subgroup a certain value Δ  and the corresponding probability The probability  Δ ⋆ () is unknown, which means that the number of elements  Δ, of each group is to be found.Similarly, for the parameter , we have   groups with the probabilities   ⋆ () = This problem can be viewed as a natural extension of the work [31] and leads to the identification of two independent discrete distributions.However, this approach has two main drawbacks.First, the two random variables Δ and  are assumed to be independent.This is a limitation when trying to reconstruct the two unknown distributions, since any possible correlation is a priori neglected.Second, the final identification problem is nonlinear, since the product of the two distributions would appear.This is in contrast with the case of the reconstruction of one single distribution, where the identification problem is quadratic [31].For these reasons, rather than considering two independent distributions, we work directly with the joint distribution, i.e., the system of  spins is divided into  groups and we associate to each subgroup a pair (, Δ)  and the corresponding joint probability  ⋆ () =    ,  = 1, … , , with ∑  =1  ⋆ () = 1.Now, the joint probability  ⋆ () is unknown, namely the number of elements   affected by the pairs (, Δ)  .This approach has the advantage of taking into account correlation effects and the final identification problem remains quadratic.It should be noted that these are acquired at the cost of an increase in the dimension of the unknown object(s), i.e., from two one-dimensional functions to a two-dimensional function.Finally, we point out that two independent distributions can also be treated as a specific case of joint distributions.
Since we are dealing with an inverse problem, we need to define what quantities can be observed in an experimental setting.In NMR, only the first two coordinates of the magnetization vector can be directly measured.We do not have accessed directly to the  component due to the strong constant magnetic field applied along this direction [6].We denote by Y u,(Δ,) () = [(), ()] ⊤ the projection of the Bloch vector onto the first two coordinates.Here, the dependence on The relation between the theoretical description of the dynamical system to the experimental outcome can be expressed as in which the two sides of the equation crucially depend on the control u.
In general, one control protocol is not sufficient to obtain an appropriate identification of the unknown  ⋆ , but a set of K control processes with K different control functions denoted u  ,  = 1, ⋯ , K, needs to be used.
where the vector  = (  )  =1 is taken in R , a subset of ℝ  , so that  = ∑      is a probability distribution.
We show in this study that GRA allows us to design a set of controls u  that makes (4) solvable and well conditioned.The algorithm is composed of two steps, namely an offline and an online steps.In the offline step, GRA computes the controls u  .In this step, only the theoretical model is needed without any experimental input.The derived controls are used in the online step in which the different magnetization vectors are measured and the minimization problem (3) is solved.Note that the controls are the same for any probability distribution to identify and only depend on the model system under study.Finally, we point out that in our algorithms the duration of each control pulse is considered as a variable to be optimized together with its amplitude.In particular, we assume that the controls are constant in time, i.e., u() ≡ u ∈ ℝ 2 , and that we can freely choose the control time up to a fixed maximum value   .Since the initial state is an equilibrium point, this is equivalent to turning on the control at a time  ≥ 0. We show in Sec. 5 that these hypotheses are sufficient for the different examples to identify the probability distributions.The generality of GRA allows one to tackle this situation in a straightforward manner.

Greedy reconstruction algorithms
We present in this section the GRA in its classical form and in an optimized extension called optimized GRA (OGRA).
GRA computes the controls u  and the corresponding control times   by solving a sequence of fitting-step and discriminatory-step problems.The goal of the fitting step is to identify a defect of the system, namely a nontrivial kernel of a certain matrix  introduced below, while the discriminatory step designs a new control which is aimed to correct this discrepancy and to eliminate the identified nontrivial kernel.The explicit formulation of GRA is presented in Algorithm 1 and is given in terms of the function h () defined as for any  in ℝ  .Notice that the fitting step minimizes over the full space ℝ  , meaning that ∑      does not have to be a probability dis-
1: Compute the control u 1 and the control time  1 by solving where    +1 is the ( + 1)-th canonical vector in ℝ  .

4:
Discriminatory step: Find u +1 and  +1 that solves 5: end for Algorithm 1: GRA consists of three steps.First, the initialization problem (6) computes the control u 1 and control time  1 .Afterwards, the algorithm iterates between the fitting and discriminatory step.The fitting step attempts to find two different probability distributions such that their final ensemble magnetization can not be distinguished by any previously computed control.Then, the discriminatory step finds a new control u +1 and control time  +1 that maximize the difference between these final ensemble magnetizations.A graphical representation of the workflow of GRA can be found in Fig. 1. tribution.However, this is a restrictive condition.On the contrary, it allows the algorithm to find and correct more nontrivial kernels than might be necessary.
One main characteristic of GRA is that the set Φ and its order have to be fixed a-priori.However, the choice and order of Φ can have a crucial impact on the outcome of the algorithm as shown in [32,Sec. 5.3].Hence, the idea of OGRA, which is stated in Algorithm 2, is to make the algorithm independent of the choice and order of the set Φ. Additionally, it aims at avoiding the computation of unnecessary control functions.This is achieved by two adaptations in GRA.The first one is that in each step one does not only consider the next element (the next canonical vector    +1 ) in the set, but all remaining elements (the canonical vectors    + for all 1 ≤  ≤  − ) in parallel.Hence, it is also possible to enlarge the set Φ (and thus enlarge ) by additional functions   which do not have to be linearly independent.In order to progressively remove linearly dependent functions in the set and to avoid scaling issues, all remaining basis elements are orthonormalized against the already selected ones after each iteration.The second adaptation is the introduction of two tolerances tol 1 , tol 2 > 0. The first tolerance tol 1 is used as a stopping criterion.The algorithm terminates if the function value in the initialization or any of the discriminatory steps (denoted by   in Algorithm 2) is too small, thus not adding new information.The second tolerance tol 2 is used to skip the computation of a new control field in the discriminatory step, if the minimum cost function value computed by the fitting step is not small enough.If this function value is large, then there already exists a control function in the current set of controls that discriminates between the two distributions  + and ∑  =1      .Notice that setting tol 2 to a very small value is reasonable if the final identification problem is quadratic.In this case, one can prove that a nonzero cost function value in the fitting step implies that one does not need to compute a new control for the corresponding set element (compare with [32]).However, if the final identification problem is not quadratic, then it can make sense to set tol 2 to a larger value.Orthonormalize all elements ( +1 , … ,   ) with respect to ( 1 , … ,   ), remove any that are linearly dependent and update  accordingly. 12: Set k = k + 1.

17: end while
Algorithm 2: OGRA consists of the same three main substeps as GRA, with the difference that OGRA considers all remaining basis elements simultaneously in each substep.Additionally, OGRA is able to skip discriminatory steps if no new control is required (line 8-9), or stop if no control that adds new information is found at all (line 3).The orthonormalization of the set Φ in line 4 removes redundant elements and prevents numerical instabilities.A graphical representation of the workflow of OGRA can be found in Fig. 1.Routine that orthonormalizes all remaining basis elements after each iteration of OGRA.

SVD_solver
Routine that solves the fitting step problem using the singular value decomposition (SVD).

Table 2
Routines related to the reconstruction of the probability distribution.

Reconstruction routines Description generate_data
Routine that generates the experimental realizations for all computed controls.reconstr Routine that either solves the final identification problem (4) using a second order interior-point algorithm, or the compact form (compare (12) in Section 5) using a solver based on the SVD.
In conclusion, the main adaptations of OGRA in lines 3 and 8-9 allow the algorithm to reduce the number of computed controls K, meaning that K < .On the other hand, as we mentioned before, GRA is de- signed to always compute exactly K =  controls.
Flowcharts describing GRA and OGRA can be found in Fig. 1.The numerical implementation of GRA and OGRA is presented and discussed in the following sections.

Structure of the code
In this section, we provide a full list of all MATLAB functions contained in the SPIRED code.Inside the SPIRED folder the user can find the main routine used to run the SPIRED code, as well as the routines that run GRA and OGRA, and that solve their sub-problems (see Table 1).Additionally, the SPIRED folder contains the routines that generate the synthetic experimental data for the true parameter probability distribution, and that solve the final identification problem (3) (see Table 2).Notice that the both the discriminatory_step and the reconstr routine use MATLAB's fmincon-solver, which requires MATLAB's Optimization Toolbox to be installed.
There are also three subfolders labeled "Test1", "Test2" and "Test3".These contain three test problems the user can choose from."Test1" corresponds to the problem discussed in this paper."Test2" is the same as "Test1", but only considers a control in the x direction (in other words u  = 0 for all control fields).Finally, "Test3" corresponds to the problem investigated in [31], where the resonance offset Δ is fixed and one attempts to reconstruct only the control inhomogeneity parameter .
Each of these "Test" folders contains routines to set the input variables, describe the cost function and its gradient for the discriminatory-step problem, and solve the corresponding dynamical system (see Table 3).They also each contain two routines used to plot the results and condition number of the reconstruction process (see Table 4).

Usage of the code
Here we illustrate the working procedure of the SPIRED code with an example.The user needs to define the test problem in the function main, which is used to initialize the procedure.f u n c t i o n [ c o n t r o l s , bases , model , r e s u l t s ] = main  OGRA, model , o p t i o n s ) ; % STEP 4 : Compute ( s y n t h e t i c ) experimental data ; Y_exp .GRA = g e n e r a t e _ d a t a ( c o n t r o l s .GRA, model ) ; Y_exp .OGRA = g e n e r a t e _ d a t a ( c o n t r o l s .OGRA, model ) ; % STEP 5: Solve the f i n a l i d e n t i f i c a t i o n problem ; ... In particular, at the "STEP 1" the user needs to define the path of the folder containing the test routines.Then, the user can define the input variables in the function starting, which is listed exemplary for the first test problem in the following.% o p t i m i z a t i o n method f o r the f i n a l i d e n t i f i c a t i o n problem s o l v e r = ' fmincon ' ; At "STEP 1" in this function, the user can define the input variables and the path to the .matfile containing the true probability distribution  ⋆ .
The input parameters related to the problem are  • um: bound   for the absolute value of the control amplitudes; • tf: maximum normalized control time   ; • Delta0: frequency shift Δ 0 for the normalized resonance offset interval; • Delta1: width Δ 1 of the normalized resonance offset interval; • Delta_interval: interval boundaries for the normalized resonance offset Δ; • alpha_interval: interval boundaries for the control field inhomogeneity parameter ; • nr_alphas: number of grid points in the direction of  for the joint discrete parameter probability distribution of  and Δ; • nr_Deltas: number of grid points in the direction of Δ for the joint discrete parameter probability distribution of  and Δ; • nr_spins: number of spins in the system; • iterations: (maximum) number of iterations performed by GRA and OGRA; for any full basis of the discrete parameter space, the obvious choice is the total number of discretization points, which is the product of nr_alphas and nr_Deltas; • Display_GRA Display option to print information about the current iteration of GRA and OGRA in the command window; can be set to 'off' to display no output, 'iter' to show the current substep of GRA and OGRA, or 'iter-detailed' to also show the current optimization problem during the substeps of OGRA; • flag_orth: flag variable that turns the orthonormalization of the remaining basis elements during OGRA on or off; • tol_OGRA_fit: tolerance tol 2 for OGRA; • tol_OGRA_discr: tolerance tol 1 for OGRA; • tol_svd: tolerance for the SVD solver, used in the fitting step and (optionally) for the final identification problem; • solver: optimization method used to solve the final identification problem (4); can be set to 'fmincon' to solve (4) using the second-order interior-point algorithm of MATLAB's fmincon-solver, or to 'svd' to solve a compact form of the problem (compare (22) in Section 6) using the SVD solver; The .mat file has to contain the variable P_star, which is the vectorized true joint probability distribution  ⋆ .If the user is considering a true experimental (laboratory) setup, meaning that they perform real experiments for the different controls to obtain the experimental data and that the true probability distribution is truly unknown, they should replace "STEP 4" in the "main.m"file with a load command to fetch the real experimental data.
Finally, to run the code the user has to write on the MATLAB prompt the following > > [ c o n t r o l s , bases , model , r e s u l t s ] = main After the computations, the routine saves the results in the MATLAB workspace (as documented in the code) and plots the reconstructed probability distributions and their difference to the true one, as well as the condition numbers for different mesh sizes.
In particular, the results obtained by running "Test1" are the true and reconstructed probability distributions for the control fields generated by GRA and OGRA, shown in Fig. 2. In Fig. 3 we show the difference between the true and reconstructed distributions for the two control field sets.Additionally, the routine provides a table containing the exact condition numbers corresponding to GRA and OGRA.Since the solver for the discriminatory step problem is initialized with a random vector, there may be small variations in some results, without changing the overall outcome.
Examples of all figures produced by the different test problems are also provided in the folder "Results" that can be found in the corresponding "Test"-folder.There one can also find the .matfiles containing the set of random controls for each test problem, loaded in the main.

Convergence analysis
In this section, we prove that the controls generated by GRA and OGRA make possible the identification of the unknown probability distributions of the parameters Δ and , i.e., they make problem (3) uniquely solvable.
Since the set of vectors  is a convex subset of ℝ  , we deduce that the problem is uniquely solvable if the matrix  is positive definite.In the case  has a non-trivial kernel, infinitely many solutions may exist which lead to wrong probability distributions different from the experimental one  ⋆ .We stress that the non-triviality of the kernel depends completely on the choice of the controls u  and the corresponding control times   .
Using the notation ( 13)-( 14), we can now also rewrite the subproblems of GRA in terms of the matrix  .The initialization problem (6) can be written as The fitting step problem ( 7) is equivalent to where   ∶= ∑  =1  (u  ,   ) and v  ∶= [ ⊤ , −1] ⊤ .Finally, the discriminatory step problem ( 8) can be written as A direct interpretation of these reformulated problems is that each control u  generated by GRA at iteration  ensures that ⟨    | |    ⟩ > 0.
Iteratively, this implies that ⟨w| |w⟩ > 0 for any w ∈ ℝ  which is equivalent to  being positive definite.
In more details, the first control u 1 and the control time  1 are chosen by the initialization (15) such that the first upper left entry of where we used that  (u, ) is positive semi-definite for any u and .
Proof.For brevity, we identify  with 1 +  for the remainder of this proof.We start by writing ) Y u,(Δ,)  ().
We conclude our analysis by the following theorem.Proof.Let   be the solution to the fitting step problem (16) for  = 1, … ,  − 1.By Theorem 2, the vector Thus, we obtain by Theorem 1 that the matrix  = ∑   (u  ,   ) is positive definite.Hence, problem (12) is uniquely solvable by  =  ⋆ .By equivalency of problems ( 12) and (4), we obtain the result.□ Notice that, in the notation above, OGRA simply reorders rows and columns of the matrix   while attempting to find and correct its kernel.In fact, the second improvement in lines 8-9 in OGRA skips the discriminatory step only if there exists a row and column of   with index  +1 such that, by swapping  +1 and   +1 , the matrix   [1∶+1,1∶+1] is positive definite.Thus, if tol 1 is sufficiently small, one can also prove convergence of OGRA analogously to GRA.

Numerical results
We test GRA and OGRA on the setting described in Sec. 2. We choose a maximum control time of 160 ms, which corresponds to a normalized time   = 16.The shift of the parameter Δ is set to Δ 0 = 4 and the width of its interval to 4Δ 1 , with Δ 1 = 0.2.We consider two different probability distributions  ⋆ , a simple Gaussian one (see panel on the left in Fig. 2) and a step distribution with three peaks (see panel on the left in Fig. 4).They are discretized by a uniform mesh of 100 points (10 points in each direction).Similarly, we discretize the set of linearly independent functions {  }  =1 by setting  = 100 and   =     ∈ ℝ 100 the -th canonical vector in ℝ 100 .Finally, we fix the tolerances for OGRA to be tol 1 = 10 −14 and tol 2 = 10 −4 .Now, let us briefly discuss how we solve the sub-steps of the algorithms numerically.The initialization and discriminatory step problems are solved by a second-order trust-region method.For the fitting step, we use the equivalent compact form (16).The corresponding first-order optimality system is given by Since the matrix [  ] [1∶,1∶] is symmetric and positive definite, any solution to Eq. ( 21) is a global solution to Eq. ( 16).Hence, we solve the fitting step problem by solving the linear system (21) using a solver based on the SVD.This solver first computes the SVD of [  ] [1∶,1∶] , i.e., two orthogonal matrices ,  ∈ ℝ × and a diagonal matrix Σ ∈ ℝ × such that  Σ ⊤ = [  ] [1∶,1∶] .To make the method more robust against numerical instabilities, it then removes all singular values that are smaller than a given tolerance, and the corresponding columns of  and  .Finally, it computes  by setting β =  ⊤  and solving After running the algorithms, we reconstruct  ⋆ by solving problem (4).Notice that, using the notation ( 13)-( 14), the gradient of the cost function in ( 4) is given by   − ∑  Γ(u  ,   ) ⊤ Y exp u  (  ), where the columns of Γ are given by the   (u  ,   ) defined in Eq. (14).We can also immediately see that the Hessian of the cost function in Eq. ( 4) is exactly  , which is guaranteed to be positive definite by our analysis in Sec. 5. Hence, the global solution to Eq. ( 4) is given by the (unique) solution to However, in order to ensure that the coefficients of the computed solution correspond to a probability distribution (i.e., belong to R ), we add the necessary constraints and solve Eq. ( 4) with the second-order interior point algorithm of MATLAB's fmincon-solver.Nonetheless, the code includes an option to solve directly Eq. ( 22) using a SVD solver (see Section 4).Now, we run both GRA and OGRA on the canonical set {  } 100 =1 of hat functions.In contrast to [31], we do not include any additional random vectors in the canonical set for OGRA and also do not remove any elements from the set during OGRA (but still reorder them).The reason for this is that we experienced for the problem of this paper that additional random elements do not improve the results and removing elements from the canonical set does not reduce the number of controls, but is more likely to make the final identification problem numerically unstable.While GRA computes 100 controls, OGRA only designs 51 by skipping 48 discriminatory steps.We then choose  ⋆ as the Gaussian distribution in Fig. 2 (left) and compute the corresponding experimental realizations {Y exp u  (  )} K =1 for the two resulting sets of control fields, with K = 100 for GRA and K = 51 for OGRA.Recon- structing  ⋆ as described above, we obtain the coefficient vectors   and thereby the distributions   = ∑ 100 =1  ,   corresponding to GRA and OGRA, shown in Fig. 2. Looking at the errors with respect to the true distribution  ⋆ shown in Fig. 3, we observe that GRA outperforms OGRA by one order of magnitude.However, the difference is so small that it is not visible in the reconstructed distributions.Similar results are obtained for a step distribution with three peaks in Fig. 4.
To investigate the dependence of the results on the choice of parameters, we repeat the experiment for different maximum control times, widths of the Δ-interval and  = 400 mesh points.First, we take a look at the number of control fields generated by OGRA in Table 5.We observe that the number of generated control fields is increasing with decreasing maximum control time and decreasing width of the Δinterval.We also observe that the ratio between the number of GRA controls (which is equal to the number of mesh points ) and the number of OGRA controls is decreasing with an increasing number of mesh points.To validate this point, we plot the number of controls for both algorithms, for different total numbers of mesh points in Fig. 5.
An explanation of this behavior is given by the condition number of the corresponding matrices  , defined in Eq. ( 13), representing the compact form (12) of the final identification problem.The condition numbers corresponding to GRA and OGRA for the settings in Table 5 are shown in Tables 6 and 7. Based on our theoretical results for GRA and OGRA proving that K =  controls are sufficient, we also add a set of fully random controls (randomized within the given bounds     for different total numbers of mesh points.To highlight the ratio between the amount of controls, we also plot half the amount of GRA controls (dotted squares).
and   ) that has the same number of controls as GRA (i.e., K = 100 and K = 400, respectively).We observe that the condition number shows the same correlation with respect to the maximum control time, width of the Δ-interval and number of mesh points, as the number of OGRA controls.In particular, the condition number of OGRA is below 115 for all settings where OGRA computed less than 60% of the number of GRA controls.
Regarding the condition numbers, GRA and random controls show the same behavior as OGRA.The reason can be found by taking a closer look at the entries of the matrix  .It can be shown that the difference between two adjacent rows or columns of  is bounded in norm by   ,   and the mesh size for the probability distribution, i.e.,  +1 −   and Δ +1 − Δ  .The interested reader can find more details about this result in Appendix A. We conclude that, if the control bound, the maximum control time, or the mesh size (or equivalently the width of the Δ-interval) is too small, the difference between two adjacent rows/columns of  can become numerically equal to zero, implying that  has a nontrivial kernel.
In order to investigate the impact of this numerical instability on the reconstructed results, we consider again the setting of the beginning of this section (i.e., Δ 1 = 0.2 and   = 16), but for  = 400 mesh points.
The results for a Gaussian and a step distribution with three peaks are plotted in Figs. 6 and 7, respectively.We observe that all three control field sets are able to fully reconstruct the step distribution and, at least partially, the Gaussian distribution.This is because the admissible set of solutions for the final identification problem ( 4) is restricted to R .Thus, a bad condition number does not necessarily imply that it is impossible to (at least partially) reconstruct the true probability distribution.However, a good condition number guarantees stability of the numerical solver and improves the accuracy of the results.In this context, notice that if we would sufficiently increase either the control bound   , or the maximum control time   , both GRA and OGRA would show better condition numbers and be able to perfectly reconstruct also the Gaussian distribution in Fig. 6.
We observe also that, if one knows the number of sufficient control functions K = , then even completely random control fields can be able to perform similarly to GRA and OGRA controls.However, while OGRA finds automatically K (reduces the number of control fields to a sufficient amount), there is no indicator for a sufficient amount of random controls in general.Additionally, the corresponding condition numbers are in many cases worse than for GRA and OGRA, as seen in Tables 6 and 7, meaning that they are more likely to show numerical instabilities.Thus, the recommended strategy is clearly OGRA, since it is able to reduce the number of control fields by up to 50%, while accurately reconstructing the probability distributions.Lastly, we remark that making the tolerance tol 2 smaller can generally lead to even fewer controls being computed by OGRA.However, this in turn can lead to less accurate results in the reconstructed solution, meaning the user has to decide for themselves if such a trade-off is desirable.

Conclusion
In conclusion, we introduce SPIRED, a Greedy reconstruction algorithm to estimate spin distribution in NMR.We show that this approach can be used to jointly find the distribution of two Hamiltonian parameters, namely the offset term and the magnetic field inhomogeneity.We discuss the accuracy and limitations of this method through experimentally relevant numerical simulations.We provide and describe the codes allowing to reproduce the results of this paper.A proof of the algorithm convergence is also given.
This paper opens the way to a series of interesting questions in quantum control.A first step is to apply this algorithm to other areas in which an ensemble of quantum systems is used.Among others, we mention Bose Einstein Condensates in an optical lattice [40,11] or molecular rotational dynamics in gas phase [10,41].The greedy reconstruction algorithm can in principle be applied to these examples, but specific constraints related to the corresponding experimental setups would be to take into account and would require adaptations of the SPIRED code.A final stage concerns the experimental implementation of this approach which seems realistic in the near future in view of the current state of the art.

Table 6
Condition number of  for different control sets, maximum control times   and widths 4Δ 1 of the Δinterval.The total number of mesh points is  = 100 and the bound on the control is   = 10.Bold numbers indicate that the condition number is smaller than 115.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

u
and (Δ, ) has been explicitly mentioned.The corresponding experimental realization of this controlled dynamic is obtained at  =   and leads to Y exp u (  ) = [ exp u (  ),  exp u (  )] ⊤ , where Y exp u (  ) can be viewed as the average at time   of the experimental measures of all the spins of the set subjected to the control u.The coordinates  exp u and  exp u are those of this measured magnetization vector.

Fig. 1 .
Fig. 1.Graphical representation of the workflow of GRA (left) and OGRA (right).The main difference between the two strategies are the while-condition that allows the algorithm to stop when it can not separate any more states by a new control, and the if-condition that allows it to skip the discriminatory step if one of the previously computed controls already separates the two states computed by one of the fitting step problems.
Optimized Greedy Reconstruction Algorithm (OGRA) Require: A set of  functions Φ = { 1 , … ,   } and two tolerances tol 1 , tol 2 > 0. 1: Compute u 1 ,  1 and the index  1 by solving the initialization problem max ∈{1,…,} max u∈ ∈[0,  ]‖h(1) (    , u, )‖ 2 , (9) 2: Swap  1 and   1 in Φ, and set  = 1, k = 1, and   = ‖h(1) (    , u 1 ,  1 )‖ 2 .3: while  ≤  − 1 and   > tol 1 do 4: f u n c t i o n [ model , bases , o p t i o n s ] = s t a r t i n g ( ) % STEP 1: Input v a r i a b l e s ; % c o n t r o l bounds and maximum c o n t r o l time ; um = 10; t f = 16; % v a r i a b l e s f o r the unknown parameters Delta0 = 4 * pi ; Delta1 = 0 . 2 ; D e l t a _ i n t e r v a l = Delta0 + 2 * pi .* [ − Delta1 , Delta1 ] ; a l p h a _ i n t e r v a l = [ −0.2 , 0 . 2 ] ; % number o f g r i d p o i n t s f o r the unknown parameters n r _ a l p h a s = 10; n r _ D e l t a s = 10; % open the f i l e to get the input p r o b a b i l i t y d i s t r i b u t i o n load ( ' Test1 / D i s t r i b u t i o n s / Gaussian .mat ' , ' P _ s t a r ' ) % number of spins ; n r _ s p i n s = 100000; % options for GRA and OGRA i t e r a t i o n s = n r _ a l p h a s * n r _ D e l t a s ; Display_GRA = ' o f f ' ; f l a g _ o r t h = 1; % numerical parameters f o r OGRA; t o l _ O G R A _ f i t = 1e −4; tol_OGRA_discr = 1e −14; % t o l e r a n c e f o r the SVD s o l v e r in the f i t t i n g s t e p ( and o p t i o n a l l y f o r the r e c o n s t r u c t i o n s o l v e r ) t o l _ s v d = 1e −10;

Fig. 2 .
Fig. 2. The plot on the left shows the true Gaussian probability distribution for  = 100 uniform mesh points.The plots in the middle and on the right contain the reconstructed distributions for the control sets generated by GRA (containing 100 control fields) and OGRA (containing 51 control fields).

Fig. 3 .
Fig. 3. Difference between the true probability distribution  ⋆ and the distributions   reconstructed using the control sets generated by GRA (left) and OGRA (right).In brackets are the number of control fields for each set.

Fig. 4 .
Fig. 4. Same as Fig. 2 but for a step distribution with three peaks.In brackets are the number of control fields for each set.

Fig. 5 .
Fig. 5. Number of controls for GRA (dashed circles) and OGRA (solid crosses)for different total numbers of mesh points.To highlight the ratio between the amount of controls, we also plot half the amount of GRA controls (dotted squares).

Fig. 6 .
Fig. 6.Same as Fig. 2 but for  = 400 and including the reconstructed distribution for 400 random control fields.In brackets are the number of control fields for each set.

Fig. 7 .
Fig. 7. Same as Fig. 6 but for a step distribution with three peaks.In brackets are the number of control fields for each set.
Hz and  ∈ [−0.2, 0.2].For the purpose of this paper, we assume that the probability densities of  and  are unknown.The

On the basis of the experimental outputs, a straightforward way to determine 𝑃 ⋆ is to solve the minimization problem min 𝑃 ∈ℙ 𝐾 ∑ 𝑘=1 ‖Y exp u 𝑘 (𝑡 𝑓 ) −
(

Table 1
Routines related to the greedy reconstruction algorithm.

Table 3
Routines related to the test problems.

Table 5
Number of controls computed by OGRA for a control bound   = 10, and different numbers of discretization points , maximum control times   and widths 4Δ 1 of the Δ-interval.Bold numbers indicate that the number of OGRA controls is less than 60% of the number of GRA controls.Notice that GRA always generates  controls.

Table 7
Same as Table6but for a total number of mesh points  = 400.