Iterative reconstruction methods and the resolution principle for fast-ion loss detector measurements

Fast-ion loss detectors (FILDs) are crucial for analyzing fast-ion dynamics in magnetically confined fusion plasmas. A core challenge is to derive an accurate ion velocity distribution, requiring treatment of thousands of remapped camera frames for a full discharge. The ill-posed nature of this task necessitates regularization with a well-chosen regularization parameter and computationally efficient methods. In this work, we introduce the ‘resolution principle,’ a novel criterion for selecting the optimal regularization parameter, providing a distinction between genuine features and artefacts smaller than the diagnostic resolution in the reconstruction, thereby preventing misinterpretations. This principle, coupled with three iterative reconstruction techniques—Kaczmarz’s method, coordinate descent, and Cimmino’s method—demonstrates enhanced reconstruction capabilities compared to conventional methods like Tikhonov regularization. Utilizing these techniques allows rapid processing of measurements from full discharges, removing the computational bottleneck and facilitating between-discharge reconstructions. By reconstructing 6000 camera frames from an ELMy H-mode discharge at ASDEX Upgrade, we capture the temporal evolution of gyroradii and pitch angles, unveiling a direct correlation between pitch-angle behavior and changes in the toroidal magnetic field for a specific subset of lost ions accelerated by edge-localized modes (ELMs) to energies approximately twice that of the injection energy.


Introduction
Fast ions are integral to the heating and stability of magnetically-confined fusion plasmas, transferring energy and momentum to the bulk plasma [1,2].These ions, produced primarily through neutral beam injection (NBI), electromagnetic wave heating in the ion cyclotron range of frequencies (ICRF), and from fusion reactions, have energies spanning from tens of keV to several MeV.Notably, alpha particles from D-T reactions, with birth energies of 3.5 MeV, are expected to be the dominant heating agent in future self-sustaining fusion reactors.Precise diagnostics of these fast ions, especially given their broad energy spectrum, are critical for optimizing plasma performance and advancing toward viable fusion energy generation.
Optimizing fusion device performance requires a comprehensive understanding of fast-ion loss mechanisms.Fast-ion loss detectors (FILDs) measure the velocity distribution of lost fast ions, providing data for interpreting the confinement and transport dynamics within the plasma [19,20].To identify the spatial origins of these losses, reverse orbit-following techniques are used, tracing the trajectories of ions measured by FILDs back into the plasma [21][22][23][24].Interpreting the velocity distributions obtained from FILDs involves solving an ill-posed inverse problem, typically approached through regularization techniques.
This paper is organized as follows.Section 2 establishes the relation between FILD measurements and the fast-ion velocity distribution.Section 3 describes Tikhonov regularization and outlines three iterative reconstruction methods used in this work.Section 4 evaluates the performance of these methods on synthetic FILD data.Section 5 assesses the iterative reconstruction methods using established stopping criteria, then introduces the resolution principle for selecting the optimal iteration number.Section 6 applies the iterative reconstruction methods to experimental data from a full discharge at ASDEX Upgrade, employing the resolution principle.Section 7 concludes the paper.

Velocity distributions from fast-ion loss detector measurements
In a scintillator FILD setup, fast ions enter through a pinhole in a probe located in the far scrape-off layer (approximately 3 cm from the plasma boundary) and strike a scintillator plate, triggering photon emission.These photons are then detected by a camera system.The captured image is the raw FILD measurement, denoted M raw .Our analysis utilizes measurements from FILD1 at ASDEX Upgrade; figure 1 illustrates its position in the toroidal and poloidal cross-sections of the device.
Each pixel in M raw corresponds to a distinct spatial location on the scintillator plate, making the raw FILD measurement a function of spatial coordinates (y, z).At ASDEX Upgrade, FILD measurements have traditionally been presented in terms of gyroradius, ρ, and pitch angle, α [5,6,15,20,25,26].This representation, denoted M f , is obtained by mapping M raw onto velocity space and is termed the 'mapped FILD measurement' [27].Subsequent discussions pertain solely to M f , hereafter referred to as 'FILD measurements'.
The gyroradius and pitch angle are related to the parallel and perpendicular velocities v ∥ and v ⊥ through the equations where v = v 2 ∥ + v 2 ⊥ , I p is the plasma current, and B t is the on-axis toroidal magnetic field.Here, ρ is explicitly expressed as a function of the ion's total velocity v, effectively serving as a normalized momentum.
Thus, FILDs quantify the velocity distribution of lost fast ions.However, ions with the same gyroradius and pitch angle may strike different locations on the scintillator due to differences in spatial coordinates and gyrophases as they traverse the pinhole.This leads to a smeared representation of the actual velocity distribution, necessitating a reconstruction process that delineates the true ion distribution from observed measurement distortions.

Reconstructions of FILD measurements
The relation between a FILD measurement M f and the velocity distribution f is modelled as a Fredholm integral equation of the first kind [25,27]: Here, W(ρ ′ , α ′ |ρ, α) represents the instrument response function, also referred to as the 'weight function' [28].Unprimed coordinates refer to the FILD measurement, and primed coordinates to the lost fast-ion velocity space.FILDSIM [25,27,29] simulations provide W, and with M f measured by the FILD, solving for f constitutes the inverse problem.
We discretize velocity space on a uniform grid and use a midpoint quadrature rule to approximate the integral in (2) [30], recasting the integral equation into a discrete matrixvector form: Here, y represents the noisy FILD measurement, ϵ the noise vector, and y = Ax the ideal noise-free signal such that y = y + ϵ.The system matrix A ∈ R m×n , composed of weight functions as rows, maps the discretized velocity distribution onto the measured signal.The dimensions m and n of A depend on the resolution of the velocity-space discretization.This model is analogous to those used for confined fast-ion diagnostics [17,28,[31][32][33][34][35][36][37][38][39][40][41].Thus, the system matrix A contains the forward model for each datum in y; specifically, the ith row of A represents how the fast-ion velocity distribution x generates the noise-free point y i .Hence, the task of reconstructing the velocity distribution translates to finding the optimal velocity distribution estimate x that yields the best agreement with the forward models across all measurements The noise in our model is represented as the sum of two independent normally-distributed noise components, ϵ i ∼ N 0, σ 2 i C ϵi , for i = 1, 2, and is assumed to be additive.Specifically, the measurement noise ϵ 1 and the background noise ϵ 2 are modelled as random vectors with covariance matrices C ϵ1 = diag y 2 , C ϵ2 = I, where σ 1 = 0.1 and σ 2 = 0.01 max y .The noisy data y is thus a combination of the noise-free signal and these noise components: (4)

Regularization methods
The inverse problem in (3) is ill-posed.As a consequence, the estimated solution x is highly sensitive to noise in the data.This sensitivity can be quantified by an upper bound on the relative error ∥x − x∥/∥x∥ [42], where x is the true velocity distribution and ∥ • ∥ denotes the Euclidean norm.In practice, this error bound is typically so large that a least-squares solution is divergent, necessitating the use of regularization techniques to stabilize the solution.

Tikhonov regularization
Tikhonov regularization is a well-established method for velocity-space tomography [34,[43][44][45][46][47][48][49] and for the reconstruction of velocity distributions from FILD measurements [25][26][27].The optimal estimate x is obtained by solving The matrix L imposes a penalty on deviations from certain solution characteristics, and the regularization parameter λ balances the trade-off between data fidelity ∥Ax − y∥ and solution regularity ∥Lx∥.The effectiveness of Tikhonov regularization in reconstructing fast-ion velocity distributions from FILD measurements has been established through both synthetic and experimental studies [25,26,50].

Iterative reconstruction methods
Iterative reconstruction methods can solve inverse problems as formulated in (3).These methods encompass various strategies, each distinguished by its computational focus.Row-action methods refine the solution by processing individual rows of A, column-action methods by processing individual columns, and simultaneous reconstruction methods update the solution using all rows of A. Each approach has unique advantages when solving (3).
To facilitate subsequent descriptions, let r i = (r ij ) and c j = (c ij ) denote the ith row and jth column of A. The matrix A can then be expressed as a composition of its row and column vectors as The velocity distributions obtained from FILD measurements are non-negative.Therefore, it is essential that the iterative reconstruction methods incorporate a non-negativity constraint.

A row-action method: the Kaczmarz algorithm
The Kaczmarz algorithm [51,52] is a row-action method that refines the estimate x(k) by iteratively satisfying individual equations in (3): the kth iterate satisfies the equation r i x(k) = y i , where k ≡ i (mod m).The update vector, ∆x (k) , is computed as Consistent with the non-negativity of FILD velocity distributions, the algorithm projects onto R n + after each iteration, making the algorithm non-linear.Consequently, the update equation is We implement the randomized Kaczmarz algorithm, selecting the ith row according to a uniform distribution, to leverage the exponential convergence of the reconstruction error [53].Cycling through m rows constitutes a single iteration.

A column-action method: coordinate descent
Coordinate descent iteratively updates x(k) using the columns of A in a cyclic order [54].For the jth canonical basis vector e j , the update is: where k ≡ j (mod n).To minimize the objective function J(x) = ∥Ax − y∥ 2 , the step length α k in the jth direction at the (k + 1)th iteration is calculated as To ensure physically viable solutions, updates are projected onto R n + , i.e. x(k+1 Cycling through n columns constitutes a single iteration.

A simultaneous reconstruction method: the Cimmino algorithm
The Cimmino algorithm [42], a special instance of simultaneous reconstruction methods, also minimizes J(x) = ∥Ax − y∥ 2 .Since the Hessian HJ x = 2A T A is positive semidefinite, the optimization problem is convex and amenable to gradient-descent solutions.The iteration step is given by Introducing a diagonal matrix D ∈ R n×n and a weighting matrix M ∈ R m×m , the iterative equation generalizes to For this study, we set D = I n×n and M = diag 1/m r i 2 , aligning with the traditional formulation of the Cimmino algorithm.Given the conditions n < m and the vector y lying outside the range of A, typical of our FILD data scenarios, the algorithm reliably converges to the weighted least-squares solution.At convergence of the Cimmino algorithm, the iterative process reaches a stable point, implying x(k+1) = x(k) = x.Thus, we have which can be rearranged as This represents the normal equations for a weighted leastsquares problem with solution This derivation outlines the convergence properties of the Cimmino algorithm, specifically for our FILD data scenarios.Ensuring physical validity by projecting onto R n + , the update equation is

Reconstructions of synthetic FILD measurements
In this section, Tikhonov regularization and the three iterative reconstruction methods, implemented using AIR Tools II [55], are evaluated on synthetic FILD measurements.The velocity distribution for these measurements, simulating three monoenergetic ion populations with ρ = 3.1, 4.3, and 5.4 cm, emulates prompt losses from deuterium NBI at 70 keV with a magnetic field strength of 1 T. The NBI power fractions are approximately 70%, 20%, and 10%, respectively [56].The pitch angle α is normally distributed with a mean of 55 • and standard deviation σ α = 7 • .The corresponding synthetic FILD measurement incorporates 10% normally-distributed noise and an additional background noise with an amplitude of 1% of the maximum ion count.Figures 2(a) and (b) depict the true fastion velocity distribution and the synthetic FILD measurement, respectively.For comparison with experimental FILD1 data, see figure 12.
In figures 2(c) and (d), we present the optimized reconstructions obtained using Tikhonov regularization and coordinate descent for λ = 1.1 × 10 −2 and k = 100.The determination of the optimal regularization parameter and iteration number involved an exhaustive search across a predefined range of λ and k values, selecting the values that yielded the best solution.The Tikhonov-regularized solution fails to adequately resolve distinct ion populations at ρ = 4.3 cm and 5.4 cm.Conversely, coordinate descent accurately reconstructs ρ and α, resulting in improved separation and precision.Reconstructions using the Kaczmarz and Cimmino algorithms are of comparable quality.It is a general finding of this study that iterative reconstruction methods outperform Tikhonov regularization in terms of precision and resolution for a broad range of FILD measurements, encompassing both synthetic and experimental data.Hence, further discussions in this paper will concentrate exclusively on the iterative reconstruction methods.For two examples showcasing reconstructions obtained using the iterative reconstruction techniques for k = 100 and k = 1000 iterations on synthetic FILD measurements, see figure 3.
The error histories ∥x (k) − x∥/∥x∥ for each method are illustrated in figure 4, indicating the relative convergence of the estimated solution x(k) to the true solution x.For k ⩽ 10 4 , the error decreases monotonically.The randomized Kaczmarz algorithm exhibits a lower error for k < 10 3 , after which coordinate descent yields a marginally superior error reduction.The Cimmino algorithm demonstrates a lower convergence rate than both the Kaczmarz algorithm and coordinate descent.
The computation time, t, for each algorithm scales linearly with the iteration number, as depicted in figure 5.The intercept β 0 of the linear regression t = β 0 + β 1 k represents the initialization time for the matrix A in computer memory and, for the Cimmino algorithm, validating the relaxation parameter value.The comparatively slower performance of the Kaczmarz algorithm, indicated by β 1 = 0.061 s/iteration, is a result of MATLAB's cache memory management inefficiencies.The coordinate descent and Cimmino algorithms, with β 1 = 0.016 s/iteration and 0.021 s/iteration respectively, are more computationally efficient.Given its computational efficiency and lower error for k > 10 3 , the coordinate descent algorithm is selected for subsequent reconstructions.
Employing iterative reconstruction techniques markedly reduces the time required for FILD data analysis.Tikhonov regularization typically requires several minutes to several hours per reconstruction of a single camera frame.A full discharge, comprising up to several thousand frames, can be processed in a matter of hours using the iterative reconstruction techniques.A subset of 100 frames can be computed in a few minutes.Thus, selecting a subset of FILD measurements, possibly at regular intervals from a discharge, enables between-discharge reconstructions that can be instrumental in fine-tuning experimental parameters.Consequently, the computational demand of inversion techniques no longer constitutes a bottleneck in experimental analysis.

Semi-convergence and the number of iterations
In the application of iterative reconstruction methods for solving linear systems of equations, a phenomenon known as semi-convergence is observed.Initial iterations tend to yield estimates x(k) that approach the true solution x, as they are less affected by the noise present in y.As k increases, x(k) approaches xnaive = A + y, leading to potential divergence due to the amplification of noise by the pseudoinverse A + .If the iteration count exceeds the optimal number, the algorithm risks entering a phase where it begins to fit the noise in the data, thereby deviating from the true solution and compromising the accuracy of the reconstruction.
To analytically examine the underpinnings of semiconvergence, consider the total error x(k) − x as the sum of two distinct components: where x (k) is the kth iterate corresponding to a synthetic noisefree measurement.The first term quantifies the noise-induced error, E n , defined as x(k) − x (k) , and the second term quantifies the iterative approximation error, E i , defined as The optimization problem is to minimize the total error by optimizing the trade-off between E n and E i to determine the optimal iteration k 0 at the point of semi-convergence.For k < k 0 , the iterations yield underfitted solutions characterized by oversmoothing; for k > k 0 , they result in overfitted solutions characterized by amplified noise.The semi-convergence behavior of the coordinate descent algorithm for synthetic FILD data is characterized by distinct phases, with the first two phases depicted in figure 6.For k ⩽ 10 3 , E i decreases at a faster rate than E n increases, indicative of improved solution fidelity.Beyond 10 3 iterations, the rate of reduction for both E i and E n diminishes, leading to a quasi-stable state which is a consequence of the coordinate descent algorithm's slow convergence on FILD data.However, the early emergence of E n imposes errors that may result in significantly flawed reconstructions.The determination of when to terminate iterations is therefore essential to mitigate the impact of the noise error.The following subsection systematically examines various stopping criteria to determine the optimal iteration k 0 for the analysis of FILD data.We note that while these criteria are fundamentally parameter choice techniques, their application extends to determining stopping criteria in our context, attributable to the direct relationship between the number of iterations and the regularization strength in our iterative reconstruction approach.

Stopping criteria for iterative reconstruction methods
The determination of the optimal iteration k 0 depends on factors such as the noise level, grid resolution, and the distribution of emissive regions in the FILD measurement.To address the resulting large variance in k 0 , various methods, primarily serving as parameter choice techniques, have been proposed and adapted as stopping criteria.Here, we examine the L-curve criterion [57], Morozov's discrepancy principle (DP) [58], and the normalized cumulative periodogram (NCP) [59], each facilitating the identification of k 0 by leveraging distinct aspects of the reconstruction problem.
The L-curve criterion is based on a log-log plot of the solution norm ∥x∥ versus the residual norm ∥y − Ax (k) ∥, which characteristically forms an 'L'-shaped curve.The 'corner' of the L is posited as the optimal iteration point.However, in the context of reconstructing velocity distributions, the criterion is prone to yielding overregularized solutions [43].
Morozov's DP determines the smallest iteration number such that the residual norm aligns with the noise level η = ∥y − y∥.Specifically, the optimal iteration k satisfies the inequality with τ typically set close to 1, acting as a multiplicative safety factor to account for uncertainty in the estimation of the noise level.The DP requires an accurate determination of η, which is its primary limitation.The NCP criterion analyzes the spectral content of the residual through its periodogram.The periodogram for the ith frequency component at the kth iteration is given by where x(k) i represents the ith Fourier coefficient of x(k) , and q is the largest integer satisfying q ⩽ ⌊m/2⌋.The jth component of the NCP vector, c x(k) , is then computed as The comparison is made against the standard reference c white , representing the idealized behavior of white noise, where c white,j = j/q, j = 1, 2, . . ., q.The optimal iteration k 0 is identified by which signifies that x(k0) has effectively minimized the deviation in the spectral content of the residuals from that of white noise.
Reconstructions of the synthetic velocity distribution in figure 2(a) determined using the L-curve, DP, and NCP criteria are displayed in figures 7(a)-(d), with the corresponding iteration numbers annotated on the L-curve in figure 8.The 'resolution principle' also appears in these figures, which will be discussed in detail in section 5.2.At a noise level of 10% and with τ = 1, the L-curve and DP criteria suggest k 0 = 15 and k 0 = 20, respectively.These small iteration numbers result in distributions deviating significantly from the true distribution, characterized by oversmoothing.Specifically, the two ion populations with ρ = 3.1 cm and 4.3 cm merge, and the ion population at 5.4 cm is too wide, suggesting that more iterations will yield more accurate reconstructions.
According to the NCP criterion, k = 1500 is the optimal iteration, consistent with the error histories in figure 3.However, at this iteration, pitch-angle artefacts within high fast-ion density regions becomes apparent, attributable to the measurement noise as the artefacts are absent in reconstructions from synthetic noise-free measurements.The artefacts are the cause of the increase in E n in figure 6 at low k values, suggesting a saturation beyond a certain iteration number since E n remains at approximately the same level.Such artefacts may lead to erroneous interpretations of the velocity distribution.Recognizing the need for a more reliable stopping criterion, we introduce the resolution principle.

The resolution principle: introduction
Consider two distinct lost fast ion populations, each characterized by their centers (ρ 1 , α 1 ) and (ρ 2 , α 2 ).The velocity distributions of these populations manifest as distinct peaks in velocity space, as illustrated in figures 9(a)-(d).In these examples, the challenge lies in resolving the two peaks from the corresponding FILD measurements.
The optimal reconstructions using coordinate descent are presented in figures 9(e)-(h).In subfigures (e) and (g), the two peaks are resolved, whereas in ( f ) and (h), they are not.This observation implies the existence of two resolution limits ∆ρ = |ρ 1 − ρ 2 | and ∆α = |α 1 − α 2 | in the ρ and α directions, below which two proximate peaks in the ground truth are indistinguishable in a reconstruction.Peaks in a reconstruction whose separation is below the resolution limits are considered erroneous, which we term artefacts.Identifying whether peak structures in a reconstruction correspond to actual ion populations or artefacts is essential for the integrity of our analysis.This distinction informs the degree of regularization necessary to produce a faithful reconstruction without overfitting.
For any point in velocity space, we seek to determine the resolution limits ∆ρ and ∆α.The resolution principle mandates that reconstructed peak structures be spaced apart by at least ∆ρ and ∆α in the ρ and α directions to ensure their physical plausibility.
The resolution principle provides a framework for selecting the regularization parameter.It advocates for the least amount of regularization while respecting the resolution limits across velocity space.This strategy calibrates the trade-off between adhering closely to the measured data and preventing the emergence of non-physical artefacts, thereby safeguarding the integrity of the reconstructed ion population peaks.

The resolution principle: definition
Let V represent velocity space.For each point x ∈ V and a given noise level η ∈ R + , we define the resolution map ϕ : where ∆ρ and ∆α denote the resolution limits in the ρ and α directions at x.These limits are functions of the noise level η.Define the set Λ = [a, b] ⊆ R + of regularization parameters λ, with a, b ∈ R + and a < b, such that for each λ ∈ Λ, a reconstruction obtained with this parameter ensures that all peaks are separated by at least ϕ(x, η) for every x ∈ V. Similarly, define the set K ⊆ N of iterations k where the reconstruction obtained after k iterations of an iterative reconstruction method satisfies the same separation criterion.The optimal regularization parameter λ 0 and iteration k 0 are determined as the minimal and maximal elements of these sets: Note that a smaller λ corresponds to less Tikhonov regularization, whereas a larger k captures more details of the data and is therefore equivalent to less regularization.

The resolution principle: algorithm
The resolution limit at each point in the discretized velocity space is determined by a post-reconstruction analysis of the spatial distribution of reconstructed intensities.This analysis entails the identification of local maxima, which correspond to the centroids of high-intensity regions, some of which are presumed to correspond to the original point sources.This identification is achieved by a multi-step process, which is schematically represented in figure 10.First, we discretize the relevant section of velocity space into a grid with uniform spacing.At a single grid point (i, j) in a velocity distribution x, we position a Dirac delta function.A second Dirac delta function is placed at a different grid point (i + k, j) or (i, j + k), spaced k grid units from the initial Dirac delta function.This velocity distribution x is mapped to the noise-free signal y via the system matrix A, after which noise is added according to (4).
The coordinate descent algorithm is then applied to the noisy data for 1000 iterations, resulting in a reconstructed solution x.This output is then reshaped into a 2D array, facilitating the transition from the algorithm's raw output to a format amenable to resolution limit analysis.The ensuing step verifies the accuracy of the reconstruction by identifying all local maxima in x, followed by isolating the two maxima with the highest intensity, which are expected to represent the positions of the original Dirac delta functions.Precisely locating these maxima is essential to determine the separation of the reconstruction peaks.
A point xij ∈ x is a local maximum if and only if its value exceeds the values of all adjacent points in its Moore neighborhood where the Chebyshev distance D Chebyshev between two points x and y with coordinates x i and y i is defined as equivalent to ∥x − y∥ ∞ .Thus, a local maximum satisfies To encompass scenarios where a local maximum spans multiple contiguous grid points, forming a plateau, we extend the definition of a local maximum to include cases where its value is equal to or greater than that of its neighbors.Therefore, we define a local maximum in x as satisfying The initial step of the algorithm involves transforming x into a binary image, where each local maximum is distinguished as '1', with all other elements marked as '0'.This transformation facilitates subsequent processing steps, which align with the methodology outlined in [60].Following this, run-length encoding [61], which compresses sequences of identical elements-termed 'runs'-into a concise format, is applied to the binary image.Each run is characterized by a value ('1' or '0') and its length, the count of consecutive occurrences.
After compressing the binary image using run-length encoding, each run indicating a local maximum (sequences of '1's) is assigned a preliminary label.Concurrently, an equivalence table is constructed to track and manage the relationships among these labels.When adjacency is detected, corresponding labels are marked as equivalent in the table, signifying their belonging to the same spatial cluster.
The next step revisits the encoded image, updating each run with its corresponding label.This results in the creation of a label matrix L. In this matrix, each element L ij is associated with a unique label.The set of all unique labels partitions x into a collection of non-overlapping regions C k , each characterized by a distinct label.
The final step involves computing the centroid and maximum intensity for each labelled region.The centroid of a labelled region C k is a vector κ k ∈ R 2 representing the weighted average position of the region, with the weights being the intensities of the points in x, so where χ C k (i, j) is the characteristic function of C k , defined as We further define M k , the maximum intensity within C k , as The regions corresponding to each label in L are ordered by descending maximum intensity to rank peaks most likely corresponding to true ion losses.
To establish the resolution limit, the two centroids with highest M k surpassing an intensity threshold of 0.15 are selected.The algorithm for the pitch-angle resolution-map calculation requires centroids to be at least one grid unit apart in pitch-angle and no more than three in gyroradius, to ensure that they correspond to distinct regions at the correct energies.For the gyroradius resolution-map calculation, the selected centroids must be separated by at least one grid unit in the gyroradius dimension and at most six units in the pitchangle dimension, to ensure that the centroids are representative of regions with approximately correct pitch angles.The distance between these two centroids in their respective dimensions is then calculated.If this distance is greater than or equal to the separation k of the original point sources, this distance is recorded as the gyroradius or pitch-angle resolution at that grid point.If the calculated distance does not meet or exceed the separation k, we increment k by one and reassess, iterating this process until the resolution criteria are satisfied.The procedure is repeated for all grid points in the discretized velocity space to compile the resolution maps.
Upon calculating the resolution maps under specific experimental conditions, the determination of k 0 for FILD measurements involves a post-processing evaluation of the reconstructions for a range of iteration numbers.Each iteration is scrutinized to ensure that the separations between discernible peaks in the velocity distribution comply with the resolution limits indicated by the resolution maps.Consequently, the optimal iteration k 0 is selected as the highest iteration number that meets the resolution criteria for all points in velocity space.

The resolution principle: example
This subsection discusses the resolution maps corresponding to the experimental conditions of discharge #35336 in ASDEX Upgrade.
Figure 11 illustrates the resolution limits for FILD1 during discharge #35336, calculated using k = 1000.The resolution limits vary across V, aligning with prior observations [25].Notably, the gyroradius resolution limit increases with increasing ρ, and the pitch-angle resolution limit increases with increasing α.Such trends align with prior expectations, given that the width of the 1D strike-point distributions in both ρ and α increases with increasing ρ and α [25,27], leading to higher resolution limits at larger values of the two variables.
The choice of k = 1000 stems from a detailed examination of the convergence behavior of the iterative reconstruction techniques, as illustrated in figures 3 and 6.Since the solutions exhibit slow convergence, they remain relatively robust against moderate changes in k.The selection of k = 1000 is based on achieving a balance between detail in the reconstruction of the two localized peaks and minimizing the amplification of noise-induced artefacts.It is pertinent to recognize that while this choice represents a midpoint, the exact value of k plays a role in the resolution limits deduced.Furthermore, we acknowledge the potential for the resolution map computation method itself to introduce artefacts, where noise-induced features may be falsely interpreted as part of the velocity distribution.
We observe localized spots in the resolution map for ∆α, particularly around pitch angles of 45 • and 60 • .These are indicative of the intricate interplay between the fidelity of the model simulations in the system matrix A and the particular characteristics of the reconstruction algorithms employed.The identification of the resolution limit at these peaks has been verified to be accurate within the constraints of our model and algorithmic approach.
Applying the resolution principle to the coordinate descent reconstruction of the synthetic signal in figure 2(b) indicates an optimal iteration number of k 0 = 100.The resulting reconstruction, illustrated in figures 2(d) and 7(c), demonstrate high accuracy in localizing the three distinct monoenergetic ion populations in ρ and α.Further, it captures the correct relative intensities and also significantly reduces the presence of artefacts.Among the evaluated stopping criteria, the reconstruction according to the resolution principle most closely approximates the true velocity distribution.
The implementation of the resolution principle serves to avoid potential misinterpretations from artefacts in the reconstructions.By establishing well-defined resolution limits for each point in the velocity space, the principle ensures that the reconstructed features represent physical phenomena according to the model rather than noise-induced artefacts.This is particularly valuable when discerning closely spaced peaks in the velocity distribution, cf section 6.As a result, the application of the resolution principle increases the reliability of the reconstruction.

Application to experimental data
We apply coordinate descent to reconstruct FILD1 measurements from ASDEX Upgrade experiments to evaluate its performance.Firstly, the convergence behavior of coordinate descent is examined using FILD measurements from discharge #35336, validating the resolution principle.Secondly, coordinate descent is employed on the entire FILD1 signal from discharge #34614, achieving reconstructions of all 6000 FILD measurements.

Validation of the resolution principle
Discharge #35336 at ASDEX Upgrade investigated TAE actuation via electron cyclotron current drive, using a combined ICRF and NBI power of 5 MW.The FILD1 measurement at t = 1733 ms, illustrated in figure 12, identified three distinct fast-ion populations.The population around ρ = 7 cm corresponds to ions accelerated by the ICRF, whereas the populations at ρ = 4 cm and ρ = 2 cm are associated with the full and half-energy NBI injected species, respectively, all with pitch angles around 65 • .Additionally, another population at ρ = 4 cm with a pitch angle around 50 • is observed.
Optimal iteration numbers, determined by the stopping criteria considered in the previous section, vary roughly over three orders of magnitude: k = O (10), O(100), and O(1000).Figure 13 shows reconstructions for these values of k for the FILD1 measurement.For k = 10, the reconstruction is overly smooth.For k = 100, the fast-ion population at ρ = 6.5 cm remains as a single resolved blob without evidence of pitchangle splitting.For k = 1000, however, splitting below the resolution limit becomes apparent within this population.At x = (6.5 cm, 60 • ), the resolution map is ϕ(x, 0.1) = (1.2cm, 7 • ), and the two peaks are separated by 2 • .Hence, this splitting is regarded as an artefact.According to the resolution principle, the optimal iteration number is identified as k = 100.
With increasing iteration number, structural and intensity variations in fast-ion populations for different k become evident, indicative of variation in the inferred ion phase-space densities within regions of high intensity as a function of k.These dynamics due to the regularization strength are captured in figure 14(a).For increasing iteration numbers, the relative amplitude at ρ = 4.5 cm diminishes.Conversely, at ρ = 3.5 cm, it increases, suggesting ions being redistributed by the algorithm from larger to smaller ρ as a function of k.Simultaneously, the amplitude associated with the high-energy population remains constant, while its distribution in ρ shows an increased degree of localization with more iterations.
Figure 14(b) shows the ion fraction at ρ = 4.5 cm as a function of the number of iterations, computed by integrating the relative intensity curve in figure 14(a) from ρ = 3.7-5.0cm, as indicated by the dashed lines in the inset.According to the resolution principle, k = 100 is the optimal iteration number, yielding an ion fraction of approximately 0.26 of the total measured lost ions.Thus, stopping the iterative reconstruction methods at the optimal iteration number, corresponding to choosing the optimal regularization strength, is crucial to mitigate artefacts induced by noise and to obtain the correct ion fractions.

Reconstructing a full discharge
We applied the iterative reconstruction techniques to the complete set of FILD camera frames for discharge #34614, characterized by a toroidal magnetic field ramp-down from 2.0 T to 1.7 T. These techniques have allowed for the reconstruction of the velocity distribution of lost fast ions for all FILD camera frames obtained throughout the discharge, revealing a direct correlation between the pitch angle of high-energy ion populations accelerated by ELMs and variations in the toroidal magnetic field.This detailed temporal analysis has uncovered behavior not observable with previous reconstruction techniques, marking a significant advancement in our diagnostic capabilities.
Figure 15 details key experimental parameters of discharge #34614.Neutral beam injectors Q7 and Q8, each delivering 93 keV, served as the primary fast-ion source with a combined power of 5 MW.For additional fast-ion diagnostics, NBI source Q3 was intermittently operated at   The discharge serves as a test case for evaluating the iterative reconstruction methods during a B t rampdown.A decrease in B t corresponds to a lower local magnetic field at the FILD probe, thereby increasing ρ of the lost fast ions.Given constant ion energy, this results in a spatial displacement of the strike point on the scintillator.The FILD captures both the full- and half-energy components of first-orbit losses from neutral beams Q7 and Q8, enabling an assessment of the reconstruction method's accuracy in representing these features for the changing ρ.As shown in figure 16, a decrease in the magnetic field corresponds to an increase in ρ.
The energy of the reconstructed primary NBI spot remains constant over time, consistent with the anticipated behavior of first-orbit losses.A high-energy feature, manifested as a small amount of fast-ion losses at higher ρ around 6-7 cm, presents a distinct behavior for ions from NBI sources Q7 and Q8.In this feature, ions from Q7 maintain constant energy, indicating their ELM-induced energy gain is not significantly affected by the decreasing magnetic field.Conversely, ions from Q8 demonstrate an increasing energy gain as the magnetic field is ramped down, increasing from 180 keV to over 200 keV.This finding corroborates findings in [26], but now presents the energy evolution in previously unseen detail.
Galdon-Quiroga et al [50] previously noted pitch-angle splitting in the high-energy features for Q7 ions during magnetic field ramp-down, with no corresponding feature for Q8 ions.Prior analyses were constrained to a select few FILD measurements; however, our application of iterative reconstruction across the entire discharge sequence has enabled a much finer analysis.Figure 17 demonstrates the temporal evolution of pitch angles, revealing the emergence of a dominant high-energy feature around 46 • for Q7 ions at approximately 2.2 seconds into the discharge.Furthermore, the ion population with a pitch angle near 50 • shows a clear dependency on the magnetic field strength, with an increase in pitch angle corresponding to a decrease in the magnetic field.This precise correlation, discerned through our new reconstruction methods, highlights a nuanced interaction between high-energy ion losses caused by ELM acceleration and the magnetic field-a phenomenon warranting further investigation.

Conclusion and outlook
work has used different iterative reconstruction techniques to shed light on the dynamics of lost high-energy ions during an ELMy H-mode discharge at ASDEX Upgrade, characterized by a toroidal magnetic field ramp-down.Our analysis has discerned distinct behaviors in the high-energy features from NBIs Q7 and Q8: while the high-energy feature from NBI Q8 exhibits acceleration to higher energies as the magnetic field decreases, the corresponding feature from NBI Q7 remains at a constant energy level.Moreover, a specific high-energy feature from Q7 demonstrates a notable increase in pitch angle in direct correlation with the decreasing magnetic field, suggesting a localized phase-space resonant interaction.
The iterative reconstruction methods employed-Kaczmarz's method, coordinate descent, and Cimmino's method-have yielded higher quality reconstructions than those achieved using Tikhonov regularization, a state-of-theart regularization method, and also significantly expedited the analysis process.The rapid computation time, approximately 20 ms per iteration, starkly contrasts with the extended durations required for Tikhonov-regularized reconstructions, which can exceed 10 min and take up to several hours.This efficiency enables the processing of measurements from full discharges, comprising more than 6000 FILD camera frames, and facilitates timely between-discharge reconstructions.Such advancements effectively remove computational speed as a bottleneck, thereby enhancing experimental planning.
Central to ensuring the accuracy of these reconstructions is the resolution principle, a novel criterion introduced for the determination of the optimal regularization parameter.By accounting for both the resolution limits of the instrumental system and the capabilities of the reconstruction method, the resolution principle aids in discerning genuine features from artefacts within the reconstructed data.This capability is crucial for avoiding misinterpretations and ensuring that the reconstructed distributions reflect the true plasma behavior.
Looking forward, the demonstrated efficacy of these iterative reconstruction techniques, combined with the resolution principle, holds promise for broader applications across various diagnostic systems in fusion research.

Figure 1 .
Figure 1.(a) Toroidal and (b) poloidal views of ASDEX Upgrade indicating the location of FILD1.The dashed lines in (b) represent the magnetic equilibrium from discharge #35336.

Figure 3 .
Figure 3. Kaczmarz's method, coordinate descent, and Cimmino's method applied to two synthetic FILD measurements with magnetic equilibrium from discharge #34614 at ASDEX Upgrade.The top row shows the ground truths, x, and the corresponding FILD measurements, y.The middle row shows the reconstructions, x, for the first example with k = 100 iterations.The bottom row shows the reconstructions for the second example with k = 1000 iterations.Observe the artefacts seen as localized high-intensity peaks occurring in the k = 1000 reconstructions.

Figure 4 .
Figure 4. Error histories for the Kaczmarz, coordinate descent, and Cimmino algorithms.All methods exhibit a monotonic reduction in error for k ⩽ 10 4 .The Kaczmarz algorithm maintains the lowest error for k < 10 3 , beyond which coordinate descent demonstrates a marginally lower error.The Cimmino algorithm consistently yields larger errors compared to the Kaczmarz and coordinate descent algorithms.

Figure 5 .
Figure 5. Computation times for the Kaczmarz, coordinate descent, and Cimmino algorithms as functions of the iteration k.The computation time increases linearly with k, modelled as t = β 0 + β 1 k, where β 1 has units of seconds per iteration.Occasional fluctuations in the time curves reflect minor variations in computer memory availability.

Figure 6 .
Figure 6.Norm of the total error, noise error, and iteration error as functions of the iteration number k for the coordinate descent algorithm applied to the FILD measurement in figure2(b).Here, x (k) denotes the kth iterate for a noise-free measurement.

Figure 8 .
Figure 8. L-curve for the coordinate descent method illustrating the norm of the solution ∥x (k) ∥ against the norm of the residual ∥y − Ax (k) ∥ for different stopping criteria.The L-curve criterion suggests an optimal iteration number of k = 15, the discrepancy principle (DP) k = 20, the resolution principle (RP) k = 100, and the normalized cumulative periodogram (NCP) k = 1500.

Figure 9 .
Figure 9. Top row: true velocity distributions consisting of two Dirac delta functions, each represented by a single nonzero pixel, separated by the indicated distances ∆ρ and ∆α.Bottom row: coordinate descent reconstructions for k = 100 iterations and 10% noise.

Figure 10 .
Figure10.Schematic representation of the algorithmic steps to generate resolution maps from a noisy signal y, originating from two Dirac delta functions.The top row illustrates a graph-theoretical interpretation of the algorithm; see appendix.The bottom row details the run-length encoding process and subsequent generation of equivalence labels leading to the label matrix L. The run-length encoding uses m to denote the run of zeros encompassing the discretization in ρ; n represents the count of leading zeros; p denotes the number of zeros separating two groups of ones, and q signifies the count of trailing zeros.

Figure 11 .
Figure 11.Resolution limits in (a) ρ (in cm) and (b) α (in degrees) for FILD1 at ASDEX Upgrade for coordinate descent with k = 1000 iterations and 10% noise.The color scales represent the resolution values.The depicted resolution maps are specific to the experimental conditions encountered during ASDEX Upgrade discharge #35336, and are constructed using simulations for those particular conditions.

Figure 14 .
Figure 14.(a) Normalized count for FILD1 as a function of ρ in the reconstructions of t = 1733 ms during discharge #35336.(b) Ion fraction at ρ = 3.7-5.0cm as a function of iterations k.An iteration number of k = 100 suggests an ion fraction of 0.26.

Figure 16 .
Figure 16.Discharge #34614: (a) magnetic field at the FILD probe, (b) ρ of the primary NBI spot, (c) energy of the primary NBI spot, (d) energy of the high-energy feature for NBI Q7 and (e) for NBI Q8.

Figure 17 .
Figure 17.Time evolution of the high-energy feature as a function of pitch angle.The data is obtained by integrating the reconstructions over ρ = 5 − 10 cm.The black and white squares indicate the Q7 and Q8 beam source distributions.