ROAM: a Radial-basis-function Optimization Approximation Method for diagnosing the three-dimensional coronal magnetic field

The Coronal Multichannel Polarimeter (CoMP) routinely performs coronal polarimetric measurements using the Fe XIII 10747 $\AA$ and 10798 $\AA$ lines, which are sensitive to the coronal magnetic field. However, inverting such polarimetric measurements into magnetic field data is a difficult task because the corona is optically thin at these wavelengths and the observed signal is therefore the integrated emission of all the plasma along the line of sight. To overcome this difficulty, we take on a new approach that combines a parameterized 3D magnetic field model with forward modeling of the polarization signal. For that purpose, we develop a new, fast and efficient, optimization method for model-data fitting: the Radial-basis-functions Optimization Approximation Method (ROAM). Model-data fitting is achieved by optimizing a user-specified log-likelihood function that quantifies the differences between the observed polarization signal and its synthetic/predicted analogue. Speed and efficiency are obtained by combining sparse evaluation of the magnetic model with radial-basis-function (RBF) decomposition of the log-likelihood function. The RBF decomposition provides an analytical expression for the log-likelihood function that is used to inexpensively estimate the set of parameter values optimizing it. We test and validate ROAM on a synthetic test bed of a coronal magnetic flux rope and show that it performs well with a significantly sparse sample of the parameter space. We conclude that our optimization method is well-suited for fast and efficient model-data fitting and can be exploited for converting coronal polarimetric measurements, such as the ones provided by CoMP, into coronal magnetic field data.


INTRODUCTION
Modification to the polarization of light is one of the many signatures of a non-zero magnetic field in the solar corona, and more generally, in the solar atmosphere (e.g., Stenflo, 2015, and references therein). Several mechanisms producing or modifying the polarization of light have been observed and studied in the solar corona at different wavelengths including, but not limited to, the Zeeman and Hanle effects (see e.g., Hale, 1908;Hanle, 1924;Bird et al., 1985;White and Kundu, 1997;Casini and Judge, 1999;Lin et al., 2004;Gibson et al., 2016, and references therein). The former induces a frequency-modulated polarization while the latter induces a depolarization of scattered light (e.g., Sahal-Brechot et al., 1977;Bommier and Sahal-Brechot, 1982;Rachmeler et al., 2013;López Ariste, 2015). Both mechanisms allow us to probe the strength and direction of the coronal magnetic field. Coronal polarization associated with these two mechanisms is currently measured above the solar limb by the Coronal Multichannel Polarimeter from forbidden coronal lines such as the Fe XIII lines (10747 Å and 10798 Å; Tomczyk et al., 2008). For these two lines, the circular polarization signal is dominated by the Zeeman effect while the linear polarization signal is dominated by the Hanle effect (e.g., Judge et al., 2006).
Translating the polarization maps of CoMP into magnetic field maps is a challenging task. The main difficulty is that the solar corona is optically thin at these wavelengths (e.g., Rachmeler et al., 2012;Plowman, 2014). The observed signal is therefore the integrated emission of all the plasma along the line of sight (LOS). Hence, the polarization maps cannot, in general, be directly inverted into 2D maps of the plane-of-sky (POS) magnetic field. On the other hand, extracting individual magnetic information at specific positions along the LOS is extremely difficult without stereoscopic observations (e.g., Kramar et al., 2014). Another limitation is that the Hanle effect associated with the aforementioned forbidden infrared lines operates in the saturated regime (e.g., Casini and Judge, 1999;Tomczyk et al., 2008). Accordingly the linear polarization signal measured by CoMP is sensitive to the direction of the magnetic field but not its strength. Deriving the magnetic field associated with the polarization maps of CoMP therefore requires a different approach than the single point inversion that can be done with, e.g., photospheric polarimetric measurements.
The alternate approach we propose to follow is to combine a parameterized 3D magnetic field model with forward modeling of the polarization signal observed by CoMP. For that purpose, we take advantage of the Coronal Line Emission (CLE) polarimetry code developed by Casini and Judge (1999) and integrated into the FORWARD package. FORWARD 1 is a Solar Soft 2 IDL package designed to perform forward modeling of various observables including, e.g., visible/IR/UV polarimetry, EUV/X-ray/radio imaging, and white-light coronagraphic observations (Gibson et al., 2016). The goal is then to optimize a user-specified likelihood function comparing the polarization signal predicted by FORWARD to the real one and find the parameters of the magnetic field model such that the predicted signal fits the real data.
In the present paper, we develop and test a new method for performing fast and efficient optimization in a d-dimensional parameter space that may be used for converting the polarization observations of CoMP into magnetic field data. The optimization method, called ROAM (Radial-basis-functions Optimization Approximation Method) is designed to be general enough so that it can be applied independently of the dimension and size of the parameter space, the 3D magnetic field model, the type of observables (provided that one can forward model them), and the form of the likelihood function used for comparing the predicted signal to the real one. ROAM is introduced in Section 2. Section 3 describes the results of multiple applications of ROAM to a synthetic test bed as validation of the optimization method. Our conclusions are then summarized in Section 4.

METHOD
The goal of this paper is to propose a model-data fitting method to be used for near-real-time 3D reconstruction of the solar coronal magnetic field. This requires developing a fast and efficient method for searching for the set of values of the model parameters that optimize a pre-defined function quantifying the differences between the predicted (or forward-modeled) and real data. Although similar approaches are standard in engineering (e.g., Jones et al., 1998), we propose a simplified version and tailored to the context of solar physics. The proposed method, ROAM, combines the computation of a log-likelihood function on a sparse sample of the parameter space with function approximation and is based on the five following steps: Figure 1. Example of three designs generated via latin hypercube sampling (LHS) in 2D (red points). Note that each of these designs only possesses one point per column and per row, which is a special feature of LHS.
1.Sparse sampling of the parameter space is performed using Latin Hypercube Sampling (LHS; McKay et al., 1979;Iman et al., 1981). LHS is a statistical method for generating a random sample of the parameter values in a d-dimensional space. For a d-dimensional space of n d points (n is the number of points for each dimension), LHS creates a set, {x i }, of n independent points or d-vectors of the parameter space (an example is given Figure 1) that will be referred to as the design in the following. 2.The model is computed for each point, x i , of the design and used to generate the corresponding predicted observation, y (x i ), to be compared with the ground truth, y gt (which is either an actual observation or a synthetic one for test beds using analytical models or numerical simulations). 3.The set of predicted observations, {y(x i )}, is then compared to the ground truth by means of a userspecified log-likelihood function where L is the likelihood function, x i is a d-vector of the design and f is a general, user-specified, well-behaved, scalar function. Typically, the likelihood function simplifies to depend on the difference between the observations and the predicted values and function f reflects that. An explicit expression of f is given in the section of each test considered in this paper (see Section 3). 4.This log-likelihood function is then approximated using radial-basis-function (RBF) decomposition (see e.g., Powell, 1977;Broomhead and Lowe, 1988;Buhmann, 2003;Nychka et al., 2015) ( where ϕ j is the j-th RBF centered at point x j of the design, · is the usual Euclidean norm, m ∈ N is such that 2m − d > 1, and {ψ j } is a set of polynomials up to degree p in the dimension d of the problem with the constraint p ≤ m − 1. In the following, we always use p = m − 1. When periodic components of the d-space exist, the value of d must be modified for the RBF decomposition to take the periodicities into account (an example and further details on handling periodic components are provided in Appendix 2). Note that the particular choice of RBFs, ϕ j , in Equations (3) -(4) is called a Polyharmonic Spline (see e.g., Duchon, 1977;Madych and Nelson, 1990) and that the polynomial term in Equation (2) is not a regularization term but an additional term that directly comes from the definition of Polyharmonic Splines as minimizers of the energy functional V⊂R d |∇ m g| 2 dx (which is not modified by adding polynomials of order p ≤ m − 1 to g). Although required from the definition of Polyharmonic Splines, this polynomial term is particularly beneficial for improving the fitting accuracy and extrapolation away from the RBF centers x j , while also ensuring polynomial reproductibility. Note also that the a j and b j are coefficients determined from the set of n equations provided by the constraint (the detailed derivation of the coefficients is given in Appendix 1) 5.Finally, we compute the set of values of the model parameters optimizing the approximated loglikelihood function using the DFPMIN IDL routine and take it as the maximum-likelihood estimator of the set of values optimizing the exact log-likelihood function. To ensure the reliability of the maximum likelihood estimator (MLE; see Section 3.2) obtained with DFPMIN, we apply the latter from (i) the point of the design that possesses the largest likelihood function value prior to step (4), (ii) N d points spanning the entire parameter space and where N ( n) is a relatively low number of points (typically N 10), and (iii) the likelihood-weighted average position of these N d points (i.e., their center of mass). Starting from these N d +2 points ensures that at least one of them will lead DFPMIN to converge towards the global maximum when the approximated log-likelihood function contains multiple global and local maxima.
An RBF is a real-valued function that only depends on the Euclidean distance to a center whose location can be set arbitrarily. RBFs provide a class of functions that possess particularly interesting properties such as continuity, smoothness, and infinite differentiability. Their use is widely spread in various branches of applied mathematics and computer science including, e.g., function approximation (Powell, 1993;Buhmann, 2003), data mining and interpolation (Harder and Desmarais, 1972;Lam, 1983;Nychka et al., 2015), numerical analysis with meshfree methods for, e.g., solving partial differential equations in numerical simulations (Fasshauer, 2007;Flyer and Fornberg, 2011;Fornberg and Flyer, 2015;Flyer et al., 2016), computer graphics and machine learning (Broomhead and Lowe, 1988;Boser et al., 1992). Polyharmonic Splines (PHS) are a type of infinitely smooth RBFs that does not possess any free parameter requiring a manual tuning. PHS can therefore be easily implemented for automated calculations.
As previously stated, the goal behind combining sparse calculations of a log-likehood with an RBF decomposition is to limit the number of model evaluations / forward calculations (n) to reduce the computational cost while maintaining a good accuracy on retrieving the exact maximum likelihood. Through low number of model evaluations, we mean to keep n 100 − 300 regardless of the dimension of the parameter space, such that all model evaluations can easily be performed at once in parallel on a highperformance computing cluster. This provides us with a significant advantage as compared with more traditional sequential optimization methods since the effective computational time of our optimization method would only correspond to the computational time of one model evaluation (because steps 4 and 5 of the method only take up to 30 seconds as long as n 500). The optimization method we propose would, in general, also be more advantageous than a full grid search. Indeed, an accurate full grid search would typically require to sample each parameter of the d-space with about 50 − 100 points at the least. This rapidly leads to a number of model evaluations that is not practical even when using parallel computing. Finally, ROAM should be competitive with genetic algorithms. Genetic algorithms applied to small population samples, e.g., a few 100 points, typically require on the order of hundred generations to converge (e.g., Louis and Rawlins, 1992;Gibson and Charbonneau, 1998, and references therein), while faster convergence would require larger population sets. For ROAM, the equivalent of a population sample is a design of the parameter space and the equivalent of a generation would be an iteration of ROAM on a smaller parameter space region. For a population/design of n-points, ROAM should, in principle, be able to converge towards the solution without the need for iterations and, hence, we estimate would be at least 50-100 times faster than a genetic algorithm with the same population/design. In practice, preliminary tests of an iterative implementation of ROAM, which will be published in a subsequent paper, show robust and accurate convergence of ROAM within a few iterations (typically < 10). show the magnetic field lines of the twisted flux rope. The red lines correspond to the magnetic field lines of the embedding magnetic field. (b) Schematic of the 3 parameters considered for our first study. The black thin solid lines highlight the solar photosphere. (θ, φ) correspond to the angular coordinates of C, the photospheric center of the 3D box containing the magnetic field of our test bed, while Ω is the rotation angle of that 3D box around the solar radial direction passing by C.

RESULTS
In this section, we present a set of test cases performed on a synthetic test bed to validate ROAM (Section 2) prior to any observational application. The set of test cases aims at assessing the performance of our method in different circumstances and defining a framework of application that will make use of its strengths.

NUMERICAL SETUP FOR THE FORWARD CALCULATIONS
Our goal is to use the proposed optimization method for data-constrained modeling of the solar coronal magnetic field using, in particular, coronal polarimetric observations (i.e., the four Stokes parameters, (I, Q, U, V), where Stokes I is the total line intensity, Stokes V is the circular polarization, and Stokes Q and U are the two components of the linear polarization). All our test cases are therefore applied to a 3D model of magnetic fields chosen to represent scenarios typically observed in the solar atmosphere. The considered magnetic model is that of a 3D coronal magnetic flux rope generated from a 3D MHD numerical simulation of the emergence of twisted magnetic fields in the solar corona (panel (a) of Figure 2; Fan, 2012).
For the test cases, this magnetic field is assumed to depend on four parameters, i.e., height in the corona (h; monotonically depends on the time of the MHD simulation, though not linearly), co-latitude (θ), longitude (φ), and rotation angle 3 (Ω; panel (b) of Figure 2). A series of synthetic polarimetric data, referred to as the ground truth (GT) in the following, is generated for the flux rope associated with (h; θ; φ; Ω) = (0.16 R ; 45 • ; 90 • ; 30 • ) (see Figure 3; note that both Stokes Q and U are presented in a frame of reference relative to the local vertical, or radial coordinate). All synthetic data are computed using the We use 192 points along both directions for the POS and 80 points for each LOS, leading to spatial resolutions of 7.6 and 19.3 respectively. We limit the forward calculations of the polarization signals to a radial range of [1.03; 1.5] R , i.e., the FOV of CoMP. Although the spatial resolution of CoMP is 4.5 , we restrict ourselves to a spatial resolution of 7.6 to allow for relatively fast (about 4-5 minutes on a MacBook Pro with a 2.7 GHz Intel Core i7 processor) calculations of the polarization signals while maintaining a quasi-CoMP resolution. We impose this FOV and POS spatial resolution to show that CoMP data currently carry meaningful information that can be used to constrain 3D reconstructions of the solar coronal magnetic field.
Finally, it should be emphasized that the considered flux rope possesses a strong degree of symmetry, such that B(θ; φ = 90; Ω ± 180 • ) = −B(θ; φ = 90; Ω). We will exploit these symmetry properties to test ROAM when faced with a log-likelihood function containing multiple maxima.

LIKELIHOOD FUNCTION WITH A SINGLE MAXIMUM
We first apply ROAM in the context of a 3D likelihood function possessing a single maximum. The parameters considered for this study are the co-latitude, longitude, and rotation angle, i.e., (θ; φ; Ω). We then build a likelihood function that takes into account all four Stokes parameters, i.e., I, Q, U, and V. For the set {x i } of a design, we first define the log-likelihood function for a given Stokes parameter, S = {I, Q, U, V}, up to a constant, as Table 1. Characteristics of the test with a likelihood function possessing a single maximum. Note: n is the number of points per design. t elapsed is the elapsed time for forwarding the Stokes images associated with the n points of a design in series, while t full is the total elapsed time that would be required to compute Stokes images for the n 3 points of the 3D parameter space in series. Each test case contains 100 randomly-chosen different designs. The naming convention is such that "xD" indicates the dimension of the parameter space and "Nx" indicates the number of points per design (n). The polarimetric data for the ground-truth are associated with (h; θ; φ; where k is the k-th pixel of the Stokes, S , image. The final log-likelihood function is then constructed as where the weighting coefficients w S were chosen to ensure that I, Q, U, and V similarly contribute to the log-likelihood function, which behavior would otherwise be dominated by the quantity possessing the largest values (here, Stokes I). We use w I ; w Q ; w U ; w V = 1.3 × 10 −4 ; 1.9 × 10 −2 ; 9.2 × 10 −2 ; 1.2 × 10 4 .
With the log-likelihood function defined in Equation (7), we consider 3 test cases referred to as 3DN31, 3DN301, and 3DN31ZOOM (see Table 1). These 3 test cases each contain 100 different designs and differ by the number of points in the designs (31 or 301) as well as by the size of the parameter space to allow us to investigate their role on the performances of ROAM. These test cases are designed to allow us to determine the criteria required for the method to ensure robustness and reliability of the results, i.e., such that the method provides a maximum likelihood estimator (MLE) that gives a good approximation of the parameters of the maximum of the exact likelihood function independently of the design and number of points used.
For each test case, the parameters of the RBF decomposition are d = 3, m = 3 and p = m − 1 = 2. We choose the minimum m satisfying the condition 2m − d > 1 (see Section 2). Although θ, φ, and Ω all are periodic parameters, their corresponding range is smaller than half the associated period and, hence, no periodic effect is expected. As explained Appendix 2, disregarding the periodicity and curvature of the dspace should not significantly affect the results in such circumstances. We therefore ignore the periodicity of θ, φ, and Ω in all 3D cases considered in this section, but return to the issue of periodicity in Section 3.3. Figure 4 presents 2D dispersion plots of the MLEs obtained for each one of the 100 randomly-chosen designs of the 3DN31 (red), 3DN301 (blue), and 3DN31ZOOM (yellow) cases. For the 3DN31, the MLEs are fairly weakly dispersed for the θ parameter, spanning a range of roughly 10 • . As summarized in Table 2, the root mean square (hereafter, rms) of the MLEs, θ rms , is ≈ 47.1 • , which is only ≈ 2.1 • different from θ GT = 45 • . This suggests that θ MLEs is not overly sensitive to the design used for the RBF decomposition. These conclusions contrast with both the φ and Ω parameters. Although φ rms ≈ 92.8 • is very close to the ground-truth, φ GT = 90 • , the ensemble of solutions, φ MLEs , spans the entire φ-range considered for the 3DN31. Similarly poor results are obtained for the set of Ω MLEs , whose rms is ≈ 21 • off from the ground-truth, Ω GT = 30 • (see Table 2). Figure 4 further shows that there is a strong coupling between φ and Ω. In particular, we find that Ω MLEs provide a poor estimation of Ω GT whenever φ MLEs are themselves a poor estimation of φ GT (and vice-versa). Such results trace very poor performances of our optimization method for the chosen setup of the 3DN31 case. The MLE strongly depends on the design used to perform the RBF decomposition. Hence, the MLE obtained from applying our method to a single design is not reliable for the setup of the 3DN31 case.
When comparing 3DN31 to 3DN301 in Figure 4, we can see that increasing the number of points significantly improves the performances of ROAM (see dark blue crosses). For all three parameters, the  rms is only ≈ 0.2 • off the ground-truth for 3DN301. The range spanned by the ensemble of solutions is relatively smaller than for 3DN31, ≈ 5 times smaller for θ MLEs and ≈ 3 times smaller for both φ MLEs and Ω MLEs . Again, the strong coupling between φ and Ω is still present but its effect on their uncertainties is strongly reduced as compared with the case 3DN31. Such results are not yet perfect since the set of solutions for both φ and Ω is spread over 20 • , which is relatively significant considering the range of their values. However, they do show that increasing the number of points per design strongly helps in reducing (a) the dependence of the MLE on the design used for the RBF decomposition, and (b) the effect of (even strong) coupling of parameters on their uncertainty. Increasing the number of points per design therefore strongly helps in improving the reliability and robustness of the proposed optimization method.
Compared with 3DN31, the 3DN31ZOOM case is used to investigate the effect of focusing the parameter space around a region closer to the exact maximum while keeping the number of points per design constant. Figure 4 shows that reducing the size of the parameter space is also beneficial for reducing the MLEs dispersion (see yellow crosses). The rms values for all three parameters are as close to the ground-truth as for 3DN301 (see Table 2) and the solutions are spanning a range that is ≈ 4 times smaller than for 3DN301 and ≈ 12 times smaller than for 3DN31. The MLEs are now independent of the design used for the RBF decomposition for θ and very weakly dependent on that design for both φ and Ω. The effect of the strong φ − Ω coupling on their uncertainty is again strongly reduced and even smaller than for the 3DN301 case. These very good results prove that ROAM can perform very well and provide an accurate estimation of the ground-truth parameters when the setup is suitably defined. Figure 5 displays 2D cuts of the exact log-likelihood function and the approximated ones associated with the designs of 3DN31 and 3DN301 giving the best MLEs (referred to as best cases in the following), as well as the approximated log-likelihood function associated with the design of 3DN301 giving the worst MLEs (referred to as worst case). Here, the best (worst) MLE is defined as the MLE minimizing (maximizing) the distance to the ground truth in the parameter space. The best MLE from 3DN31 is (θ; φ; Ω) best−MLE = (44.6 • ; 92.9 • ; 30.1 • ) while the best MLE from 3DN301 is (θ; φ; Ω) best−MLE = (45.0 • ; 89.9 • ; 30.3 • ). The worst MLE from 3DN301 is (θ; φ; Ω) best−MLE = (46.5 • ; 103.7 • ; 43.0 • ). The figure shows that the approximated log-likelihood function of the best case from 3DN31 gives an overall rough approximation of the exact one both in terms of values and shape. Note, though, that the rms error on the log-likelihood is 0.24, which is rather small given that max(| (x)|) ≈ 3.5 for the considered parameter space. For the log-likelihood function of the best case from 3DN301, the results are very much better. The approximated log-likelihood function is able to accurately capture both the values and shape of the exact log-likelihood function; the rms error is 0.05, i.e., ≈ 5 times smaller than for the best case of 3DN31. For the worst case of 3DN301, the rms error on the log-likelihood is 0.25, which is very similar to that of the best case of 3DN31, and the MLE is far from the ground truth for both φ and Ω. However, we find that the worst case from 3DN301 provides a more accurate RBF decomposition of the exact log-likelihood function than the best case of 3DN31; the log-likelihood function surfaces display a similar pattern as for the best case of 3DN301 but shifted in the Ω direction. The difference with the 3DN31 lies in the density of points in the entire design, and in the vicinity of the exact maximum, with regard to the structuring, or gradients, of the exact log-likelihood function. This is because the goodness of the approximation is determined by that of the RBF decomposition, which depends on the number of constraints -and hence, points -brought by the design. In other words, the more structured the exact log-likelihood function, the stronger the effect of point density on the goodness of the RBF decomposition / log-likelihood function approximation, similar to what one would expect when discretizing a continuous functions that contains strong gradients. While not shown here, the combined effect of point density and log-likelihood function structuring on the quality of the RBF decomposition is further supported and illustrated by the best case of 3DN31ZOOM that provides the best approximation of the exact log-likelihood function in the vicinity of the exact maximum even though the corresponding design only includes 31 points.
The aforementioned results show that the RBF decomposition performed well from a sparse sampling of the parameter space, and hence ROAM, is able to capture both the values and variations of the exact log-likelihood function when suitable conditions are met, namely, when the design contains a high enough density of points in the surroundings of the exact maximum and in areas where the exact log-likelihood function is strongly structured. They further demonstrate that ROAM can perform well even with a very low number of points per design although not as robustly. The combined results from 3DN31 and 3DN31ZOOM indicate that an iterative application of ROAM with a smaller and smaller parameter space would be an interesting way to improve its robustness when used with a very sparse design. Such a robust approach has been successfully tested but is beyond the scope of this work and will be presented in a subsequent paper. In particular, the iterative implementation of ROAM strongly improves φ MLEs and Ω MLEs , leading to a better than 0.5 • accuracy on both of these parameters in typically 4-5 iterations with designs of 31 points. The solution is quasi-independent of the design used for the RBF decomposition and the strong φ − Ω coupling (previously mentioned and visible in panel c of Figure 4 and in the φ − Ω cut of the exact log-likelihood function shown in Figure 5) is comfortably reduced and overcome (note that such coupling could also be overcome by separately optimizing one of the coupled parameters, e.g., Ω, and apply ROAM to the 2D parameter space (θ; φ)).
Note that, with the goal of increasing speed, and considering our previous comments on the acceptable degree of roughness in the log-likelihood function approximation and an iterative implementation of ROAM, a less sophisticated approach might be conceived. For instance, one could first go through steps 1 to 3 of the method (see Section 2). Then, step 4 (i.e., the RBF decomposition) would be replaced by taking the point of the design associated with the highest log-likelihood function value as a temporary MLE and one would iterate the procedure by defining a smaller design centered around the temporary MLE until a convergence criterion is reached. There are several reasons for not making such a choice. The main reason is that such an initial guess can be far from the exact maximum likelihood, which would likely slow down the convergence by requiring unnecessary iterations and would make the final result more sensitive to local maxima. In addition, applying the RBF decomposition and the search for the maximum from the approximated log-likelihood function is computationally cheap when the number of RBFs is as small as for the cases considered in this study, i.e., typically takes less than 10 seconds for the designs of the 3DN301 case. The benefits of applying steps 4 and 5 of ROAM as proposed in Section 2 (i.e., the RBF decomposition and the search for the MLE from the RBFs approximated log-likelihood function) are illustrated in Figure 6. The figure displays Stokes images for the ground truth, the MLE obtained from fully applying our optimization method to one design of the 3DN301, and for the initial guess from that design. As one can see, the initial MLE guess from the design, (θ; φ; Ω) IG = (46.0 • ; 101.4 • ; 45.9 • ), has a φ IG and Ω IG that are far off both the ground truth ((θ; φ; Ω) GT = (45.0 • ; 90.0 • ; 30.0 • )) and the MLE obtained from the RBF decomposition (full application of ROAM; (θ; φ; Ω) MLEs = (44.6 • ; 90.0 • ; 28.0 • )). These strong differences in φ and Ω result in significantly different Stokes profiles. Iterations would then be needed for the results to be as close to the ground truth as the MLE from the full optimization, which (1) gives a very good estimation of the parameters of the exact maximum likelihood without any real need for iterations, and (2) only takes a few more seconds of calculations.
In practice, the current capabilities of the CoMP instrument and calibration software do not allow routine measurements of Stokes V since the signal-to-noise ratio is too small. We therefore perform an additional test to show that the current linear polarization signal from CoMP is sufficient to constrain the parameters of a magnetic model using ROAM. The log-likelihood function is defined as in Equation (7) keeping w I , w Q , w U as before, but now setting w V = 0. The results of that study are displayed in Figure 7 for 3DN301. The figure presents scatter plots of the MLEs of a design obtained when using all four Stokes versus obtained when using Stokes I, Q, and U only. In such plots, the points should form a line of equation y = x whenever the solutions obtained one way or the other remain the same. As one can see  from Figure 7, this is exactly the case for θ MLEs . Most of the points are also forming a straight line, y = x, for both φ MLEs and Ω MLEs , with only about 7-8 points (out of 100) being off the line. Such results indicate that a log-likelihood function built from Stokes I, Q, and U contains sufficient information to constrain the three spatial location and orientation parameters considered here. We therefore conclude that the current linear polarization measurements from CoMP contain sufficient observational information to constrain some of the parameters of a given magnetic model.

LIKELIHOOD FUNCTION WITH MULTIPLE MAXIMA
In this section, we test ROAM in the case of a log-likelihood function with multiple maxima having similar values. For that purpose, we only build the log-likelihood function with Stokes Q and U, setting the weight coefficients of Equation (7) to (w I ; w Q , w U , w V ) = (0.; 2.0 × 10 −2 ; 6.9 × 10 −2 ; 0.). Only the height of the flux rope in the corona, h, and the tilt angle, Ω, are considered for this test (see Table 3 for the range of values considered for each parameter).
Stokes Q and U signals are associated with the transverse magnetic field, i.e., the component of a magnetic field perpendicular to the LOS. For a single point in the solar corona, the transverse magnetic field diagnosed from either the Hanle or Zeeman effect is subject to a 180 • ambiguity (e.g., Casini and Judge, 1999;Judge, 2007). In terms of the parameters considered in our tests, it means that a single point magnetic field set with a rotation angle, Ω SP , will give the same Stokes Q and U signals as when set with Ω SP ± 180 • . Considering that φ GT = 90 • (that is, the flux rope is centered at the solar limb) and the strong symmetry of our flux rope (see Section 3.1), we expect the LOS integrated Stokes Q and U to be the same for Ω and Ω ± 180 • , resulting in a log-likelihood function with two maxima respectively located at Ω GT and Ω GT ± 180 • ; note that the symmetry of Stokes Q and U would be broken if the flux rope were not centered on the solar limb. This is indeed the case as shown in panel (b) of Figure 8   two additional maximum regions located at Ω = 150 and 330 • . These two solutions suggest a symmetry with regard to the plane Ω = 0 • that is not expected. We find that the corresponding Stokes Q and U images are, as expected, different from those of the ground truth. However, the differences are small as compared with other values of Ω, resulting in a local maximum in those two regions. Note, though, that these four maximum regions are only possible because φ = 90 • , whereas any other value of φ would break the symmetry of the Stokes Q and U images.
In the present test case, the periodic parameter Ω varies on a range of values larger than half its period. In such circumstances, we must consider its periodicity for the RBF decomposition (see Appendix 2). Accordingly, the parameters of the RBF decomposition are d = 3, m = 3 and p = m − 1 = 2. As in Section 3.2, we run our optimization method on 100 different designs whose properties are given in Table 3. The results are summarized in a 2D dispersion plot in panel (a) of Figure 8. The figure shows that the 100 MLEs are mainly, and almost equally, clustering around the two global Ω maximum regions, corresponding to the ground truth and its counterpart at 180 • . We further find that, out of these 100 solutions, only 4 are associated with one of the two local maximum regions, here Ω ≈ 150 • . As for the height of the MLEs, we find an average value of 1.6 × 10 −1 R with a 2σ dispersion level of 0.5 × 10 −1 R , meaning that the height is well constrained even from using Stokes Q and U only. The dispersion plot from panel (a) of Figure 8 therefore indicates that our optimization method is strongly sensitive to multiple global maxima and can be sensitive to local maxima. Note that the sensitivity to local maxima depends upon both the number of points used in the design and the value of these local maxima relatively to that of the global maxima. Panels (c) of Figure 8 displays a surface plot of the log-likelihood function from the best case of 2DN120. As one can see, the RBF decomposition is able to capture both the values and shapes of the exact log-likelihood function. We find an rms error of 0.04 on the log-likelihood. The RBF decomposition can therefore provide a good approximation of the exact log-likelihood function even with a periodic space and the presence of multiple maxima.
Finally, we show in Figure 9 that using the Stokes V signal to build the log-likelihood function removes the Ω ambiguities that were observed in the log-likelihood function constructed from Stokes Q and U only. When using Stokes V, the optimization leads to Ω rms ≈ 28.8 • . This means that some additional observables might be worth considering to remove ambiguities in parameters when they exist. Another alternative to remove ambiguities is to reduce the parameter space to regions having a single maximum. Then, one can either study each region separately or use prior constraints to eliminate regions that are very unlikely. For instance, one can use the photospheric magnetograms, or Hα observations, prior to or after the passage of the flux rope at a limb to estimate the rotation angle (i.e., Ω) and put strong constraints on the values of rotation angle to consider for the parameter space.

STABILITY WITH REGARD TO NOISE IN THE DATA
In practice, any real data is subject to measurement errors. Such errors may prevent the retrieval of any meaningful information about the polarization, and hence, the magnetic field in regions of weak signals and/or when the signal-to-noise ratio is weak. The results of ROAM might be sensitive to such noise and we therefore need to investigate that sensitivity. For that reason, we now test our method when the synthetic observations associated with the ground truth contain some noise. In this regard, we build the log-likelihood function with Stokes U images only where σ U is the root mean square of the noise in the synthetic Stokes U signal of the ground truth.
For a given value of photon noise, σ I , σ Q , σ U , σ V are all different. As a consequence, if one uses more than one Stokes component, then varying the noise further changes the relative contribution of each Stokes parameter to the log-likelihood function due to the weighting by 1/σ S . We need to be free of the variation of relative contribution of the different Stokes in order to isolate the sole effect of noise on the robustness of our optimization method, which then implies using only one Stokes parameter to define the log-likelihood function. Considering CoMP capabilities and the current magnetic model and ground truth, we performed several tests with different levels of noise (which can be added using FORWARD) and found that (1) Stokes V cannot be used for realistic exposure times because its values for our test bed are too weak and would require an unrealistic 4 days exposure time to reach a moderate level of noise for the particular choices of ground-truth parameters and pixel sizes (see corresponding values in Section 3.1), (2) Stokes I cannot be used because it is not sensitive enough to noise (even a 1 second exposure time leads to a very weak level of noise), and (3) Stokes Q and U are better suited for the noise test with exposure times of the order of 1 to 100 seconds. From this analysis, we chose Stokes U because it was slightly more sensitive to noise than Stokes Q for the setup considered in this paper (note that both Q and U are presented in a frame of reference relative to the local vertical, or radial coordinate).
FORWARD already implements a photon noise calculation for the infra-red lines under consideration (see e.g., Gibson et al., 2016). The noise is calculated according to the specifications of the instrument considered (telescope aperture, detector efficiency), the background photon level, and the exposure time  to obtain a forward calculation that includes the noise. For CoMP, the aperture is 20 cm, the efficiency is 0.05 throughput and the background is 5 parts per million of solar brightness. We perform three tests with different exposure time, t exp , hence noise level, i.e., t exp = (1; 10; 300) seconds that respectively correspond to strong, moderate, and weak noise cases for the considered setup. The synthetic Stokes U images of the ground truth for these noise levels are displayed in Figure 10. These synthetic ground truth are used with all designs of the 3DN31, 3DN301, and 3DN31ZOOM cases. Figure 11 presents scatter plots of the error on the MLEs obtained when noise is included in the ground truth Stokes U images as compared with the case when no noise is considered. The plots are only shown for the Ω parameter because all three parameters θ, φ, and Ω display very similar results. In Figure 11, one can see a nearly perfect horizontal line at y = 0 for Ω MLEs obtained with an exposure time of 300 seconds (yellow crosses) for all test cases (3DN31, 3DN301, and 3DN31ZOOM). This means that the t exp = 300 seconds case is equivalent to the no-noise case. For cases t exp = 10 and t exp = 1 second, the plots show some departure from the y = 0 line, which increases with the level of noise. The figure also shows that the noise effect on the robustness of the MLEs depends on the density of points in the designs, i.e., 3DN31 is the most affected by the noise while 3DN31ZOOM is the least affected. That being said, we find that only less than ≈ 10 − 15 points (out of 100) of 3DN31 exhibit a strong sensitivity to noise for the t exp = 1 second case, i.e., with an error larger than 5 • . This number drops to ≈ 5 when t exp = 10 seconds. For the t exp = 1 second, Figure 10 shows that the noise strongly masks the real Stokes U signal, although not entirely. We therefore deduce that our optimization method is very stable against the presence of noise in the data as long as the noise does not entirely mask the real signal. Considering that Stokes Q is similarly sensitive to noise as Stokes U and that Stokes I is much less sensitive to the noise, we conclude that our method can be robustly used with the Stokes I, Q, and U data provided by the CoMP.

CONCLUSIONS
In this paper, we introduced and validated a new optimization method for model-data fitting, ROAM (Radial-basis-functions Optimization Approximation Method). Our primary motivation for this work has been to develop a novel approach for diagnosing the solar coronal magnetic field by combining a parameterized 3D magnetic field model with forward modeling of coronal polarization. From various tests applied to the synthetic test bed of a coronal magnetic flux rope, we showed that ROAM allows for fast, efficient, and accurate model-data fitting in a d-dimensional parameter space. These test cases further enabled us to analyze and specify a framework for an optimal application of ROAM.
Applying our method with forward modeling of IR coronal polarimetry, we demonstrated that ROAM can be exploited for converting coronal polarimetric measurements into magnetic field data. The use of our model-data fitting method therefore opens new perspectives for the development and exploitation of coronal polarimetric measurements such as the ones routinely performed by CoMP (Tomczyk et al., 2008) and future telescopes such as the Daniel K. Inoue Solar Telescope 4 and the Coronal Solar Magnetism Observatory (Tomczyk et al., submitted), but also for a wider range of coronal observations including, e.g., UV (see e.g., Fineschi, 2001;Raouafi et al., 2009) and radio polarimetry (e.g., White and Kundu, 1997;Gelfreikh, 2004; see also Gibson et al., 2016, for discussion of multiwavelength magnetometry).
Beyond the analysis of coronal polarimetric measurements, ROAM offers interesting perspectives for magnetic field reconstruction models. Most of the current 3D diagnostics of the coronal magnetic field of solar active regions (ARs) are derived from the analysis of magnetic field reconstruction models including, e.g., force-free field extrapolations of the photospheric magnetic field (see e.g., Alissandrakis, 1981;Demoulin et al., 1989;Wheatland et al., 2000;Yan and Sakurai, 2000;Wiegelmann, 2004;Amari et al., 2006;Malanushenko et al., 2012, and references therein), and magneto-frictional methods (see e.g., van Ballegooijen, 2004;Valori et al., 2005Valori et al., , 2007Jiang et al., 2011;Inoue et al., 2012;Titov et al., 2014, and references therein). ROAM could, in principle, be used to perform model-data fitting with such reconstruction models that either already are (i.e., through the poloidal and axial flux for the magneto-frictional methods with flux rope insertion) or could be (e.g., through the photospheric forcefree parameter for both force-free field extrapolations and magneto-frictional methods without flux rope insertion) parameterized. The extensive work performed over the years in terms of forward modeling of various observables (see e.g., Gibson et al., 2016, and references therein) would then allow for using several types of different observations to constrain the parameters of the magnetic field reconstruction models. ROAM therefore opens new perspectives for including coronal polarimetric measurements into magnetic field reconstructions and, more generally, for data-optimized reconstruction of the solar coronal magnetic field. Such perspectives will be tackled in the framework of the Data Optimized Coronal Field Model 5 (DOCFM), a collaborative project that will make use of ROAM.
Finally, we wish to mention that ROAM is not limited to coronal magnetic field diagnostics and could be used for other optimization problems. The method will be of particular interest for model-data fitting for which a model evaluation (here, the evaluation of the model itself and/or the forward modeling of an observable if applicable) is computationally expensive.

DISCLOSURE/CONFLICT-OF-INTEREST STATEMENT
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

AUTHOR CONTRIBUTIONS
All the authors contributed to the building, writing, and editing of the content of the paper.

ACKNOWLEDGMENTS
We thank the two anonymous referees and Anna Malanushenko for a careful consideration of the manuscript and constructive comments. The National Center for Atmospheric Research is sponsored by the National Science Foundation.

APPENDIX: SOLVING THE COEFFICIENTS OF THE RBF DECOMPOSITION
The a j and b j coefficients of the RBF decomposition (Equation (2)) are determined from the constraint (x i ) = (x i ), which leads to the n following equations where N = p+d p ; note that N n. Equation (9) can then be re-written as where Equation (10) is the i-th row of the matrix equation where L and a are n-vectors, b is a N-vector, Φ is a n × n matrix, and Ψ is a n × N matrix.
Equation (11) is expressed as the sum of an affine, Ψb, and non-affine, Φa, term. The a and b coefficients can be solved by separating the affine and non-affine terms of Equation (11) through QR-factorization of Ψ. QR-factorization or QR-decomposition consists in the decomposition of a matrix M into a product of an orthogonal matrix Q and an upper triangular matrix R, such that M = QR. The QR-factorization of our rectangular matrix Ψ leads where R 1 is a N × N upper triangular matrix, Q 1 is n × N, and Q 2 is n × (n − N) (which implies n > N). In QR-factorization with n > N, the columns of both Q 1 and Q 2 are orthonormal, which implies Q T 1 Q 1 = I and Q T 2 Q 2 = I, where T denotes the transpose. Note that Q being an orthogonal matrix implies that each column of Q 1 is orthogonal to each column of Q 2 , which means that Q 2 is in the orthogonal complement of Q 1 and Q T 2 Q 1 = 0 and Q T 1 Q 2 = 0. With that in mind, it is then possible to separate the affine and non-affine terms in Equation (11) by taking the a coefficients in the space associated with Q 2 , such that a = Q 2 γ, where γ is a N-vector. Replacing this particular solution of a and multiplying Equation (11) by Q T 2 leads to which, after some algebra, results in The a and b coefficients are then given by In all the results presented in this paper, the QR-factorization and the a and b coefficients were computed using the DGEQRF, DORGQR, and DGETRF routines of the LAPACK fortran library.

APPENDIX: RBF DECOMPOSITION WITH PERIODIC COMPONENTS
For all applications considered for our optimization method, the parameter space, d-space, always defines a set of points in the Euclidean space, E d , where d ≥ d (although the d-space may itself be non-Euclidean). The d-space belongs to E d only when it does not possess any periodic parameter. Whenever they exist, periodic directions imply that the d-space is curved, which in turn means that the set of points of the d-space is actually defined in E d with d > d. For instance, this is the case of the points defining the surface of a cylinder and/or a sphere. Although they define a two dimensional d-space, these points actually form a set of R 3 . Hence, a transformation from R d to R d should be applied to the x d-vectors if one wants the periodicities of the d-space to be included for the RBF decomposition and one must then substitute d by d in Equations (2) -(4) (see Section 2). Let x T = x T n−per. , x T per. be a vector of the d-space such that x n−per. and x per. are the sets of non-periodic and periodic components of x. If all components of x per. are independent, the coordinate transformation from R d to R d is where the i-th component of u 1 and u 2 is where P i is the period of x i,per. . For a d-space with independent periodic components, then d = d n−per. + 2d per. (where d Y is the dimension of x Y ).
When some or all of the periodic components of the d-space are coupled, the coordinate transformation depends upon the relationship between the coupled periodic components. For instance in a 2D parameter space with x = (x 1 , x 2 ), if x 1 and x 2 are the longitude and co-latitude on a sphere (x 1 ∈ [0; 2π] and x 2 ∈ [0; π]), then the coordinate transformation from R d to R d is u 1 = cos (x 1 ) sin (x 2 ) , (22) u 2 = sin (x 1 ) sin (x 2 ) , (23) u 3 = cos (x 2 ) , leading to d = 3.
Considering the periodicities is not mandatory for applying the RBF decomposition. However, not considering it may significantly alter the accuracy of the RBF decomposition, leading to a bumpy approximation and a bad MLE. Such an example is displayed in Figure 12 for one of the 2DN120 case (see details in Section 3.3). The bumps in the case that disregards the periodicity are likely to occur whenever the range of any periodic component of the d-space becomes strictly larger than half the associated period. That hypothesis is suggested by the 3DN301 and 3DN31ZOOM cases of Section 3 that give good loglikelihood function approximations while disregarding the periodicities of θ, φ, and Ω whose ranges were smaller than half the periods. The aforementioned hypothesis can also be hinted at from the following. Let x A and x B be two points in R d , x A is held fixed while x B is being moved away from x A along the periodic component only. In R d where the periodicity is disregarded, the Euclidean distance between these two points always increases monotonically as x B is moved away. That is not the case in R d where periodicity is taken into account. Indeed, the Euclidean distance between these two points periodically goes through maxima and minima when the distance along the periodic component reaches half the period and the period, respectively (like what happens when one puts two points on a circle and moves one of the two points along the circle). As a consequence of its definition that only depends on Euclidean distances, an RBF in R d is expected to display a similar cyclic behavior along the periodic components. Therefore, Figure 12. Comparison between approximated log-likelihood functions when taking the periodicity of Ω into account (left panel) and when disregarding it (right panel). The approximated log-likelihood functions are associated with the design of 2DN120 giving the best MLE (shown panel (c) of Figure 8). The white "+" symbol indicates the position of the maximum log-likelihood. Notice the bumpy behavior of the approximated log-likelihood function when the periodicity of Ω is not accounted for and the MLE which is completely different from the MLE obtained when considering the periodicity. the behavior of an RBF will be consistent, i.e., monotonic, in both R d and R d as long as the distance to that RBF along any periodic component remains smaller than half the corresponding period, and hence, similar results are expected for the RBF decomposition whether it is performed in R d or R d . This is not true anymore when the range of one or more periodic components becomes larger than half the associated period.
Finally, let us emphasize that the intuitive approach to deal with periodic components and that consists in computing the contribution of an RBF, ϕ A , at x by taking the minimum distance, r min , between ||x − x A || and ||x − x A ± P|| (where P T = 0, 0, · · · , P ) is strongly ill-advised. Indeed, it would have the undesirable effect of creating an artificial discontinuity of the RBF at ||x A ± P|| caused by the fact that our RBFs are not compact. Taking the minimum distance, r min , would only work with compact RBFs, provided that their radius of influence (i.e., the radius below which they have a non-zero value) in the periodic direction be strictly smaller than the corresponding period.