The SparSpec algorithm and the application to the detection of spatial periodicities in tokamaks: error weighting the penalization criterion to improve the performance of the algorithm^*

The problem of detecting frequencies from time series is common in a large number of fields in science, among which astrophysical observations are an important example. On most occasions, the analysis of such time series is the only feasible approach to the study of the stellar environment, such as in the context of astero-seismology, where the oscillation frequencies of stars are used to deduce their interior dynamics and verify proposed models [1, 2]. These time series are usually irregularly sampled—as a result of such constraints as the rotational period of Earth and/or practical limitations (meteorological conditions)—despite being often very large and long-term (time spans of months if not years). As a consequence, classic decomposition tools as the Fourier transform are ineffective and give wrong results. This has led to a remarkable effort in devising specific methods of frequency estimations such as phase-binning techniques [3–5], Lomb-Scargle periodogram [6–8], data-compensated discrete Fourier transform [9], and iterative methods based on pre-whitening steps [10, 11], just to mention a few.

Many similarities can be found between the astrophysical detection of temporal frequencies and the detection of the spatial periodicities of magnetic fluctuations (=modes) in a tokamak plasma and this has suggested that a similar approach to the data analysis could be successfully employed in the latter case. Nevertheless, a few important specificities of the plasma problem are to be stressed which complicate this task, namely:

• the number of detectors in the tokamak is limited and of the order of 10 → 50;
• the modes are both in the temporal and spatial domains, and their periodicities (=frequencies in the following) have rational values and are continuous in the time domain but have only discrete and integer values in the spatial domain;
• the input signals can be real-valued (a time series) or complex-valued (a frequency spectrum);
• there are practical necessities such as requiring that the computation be as fast as possible so as to conduct a real-time (RT) analysis of the modes developing during a run in the tokamak (contrary to the astrophysical application, where the data analysis can and usually does span a long period);
• finally, there is a different meaning, hence treatment, for the mean signal: in the astrophysical context this is usually discarded being normally the background noise or associated to the Earth's orbital movement, whereas in plasmas this could be associated with the spatial (n, m) =0 mode, thus retaining a physical importance, n and m being the toroidal and poloidal mode numbers, respectively.

Given the above properties pertaining to the detection of plasma modes, a fruitful approach can be found in the theory of the Sparse Representation of Signals (SRS), on which the SparSpec code [12, 13] (see: http://userpages.irap.omp.eu/~hcarfantan/SparSpec1.4/SparSpec_html.html) is based. This is its baseline and simplest version, and its application to the analysis of spatial periodicities in tokamaks has already been presented multiple times [1, 12–18]. Therefore only a rather short description is outlined in the next section to simplify the reading of this contribution.

This work aims at presenting a further development of the baseline SparSpec algorithm, dubbed SS-H2, as applied to the detection of spatial frequencies in tokamaks, namely using the measurement uncertainties to weight the data, the spectral window and thus the penalization criterion. This is part of a major update of the SparSpec algorithm, specifically intended to accelerate the RT processing, of which the second part, namely the introduction of a memory with relaxation (MwR) scheme to provide better input constraints for the analysis and thus hopefully speed up the RT calculations, is presented in a separate contribution for clarity. A third development, namely the application to a 2D spatial decomposition problem, is ongoing and will be reported separately at a later date.

Regarding the use of the measurement uncertainties to better constrain a penalization problem, other approaches obviously exist for these analysis, such as those deploying Bayesian inference techniques [19]. These may also be, perhaps even more naturally, suited for the quantification and propagation of the measurement uncertainties, but their discussion is outside the scope of this work.

Regarding the use of the measurement uncertainties within SparSpec, it has been often found in the analysis of JET and TCV magnetic data that not all input signals could be used with SS-H2, as they were characterised by rather different intrinsic errors and thus could not be assembled effectively as a unique input dataset. This partially reduced the accuracy of the algorithm, as we indeed found that sometimes different solutions could be obtained depending on which subset of input signals was used. A clear step forward is to introduce the measurement uncertainties in SparSpec, which has led to the development of the SparSpec-V5 algorithm (which also includes the MwR scheme), of which the SS-V5ν0 ⁴ version (without memory) is considered here.

This paper is organised as follows. Starting from existing literature, in section 2 we briefly review the mathematical foundations of the SRS and of the SparSpec code ⁵ , and it will become apparent that is quite straightforward to weight the input data with respect to the overall measurement uncertainties, and thus use them as well in the spectral window and penalization criterion. Section 3 then compares the results obtained with SS-H2 and SS-V5ν0 for an extensive number of simulated test cases, this being the most appropriate test as we know what the result should be. The clear improvements brought about by using SS-V5ν0 versus SS-H2 using simulated data have been verified on one test case using actual data from the JET tokamak, as shown in section 4. These performance improvement will have to be more systematically addressed in as many experiments (JET and TCV, and other devices as well) as possible, when we do not know what the results should be; however, this activity is left for future work. Section 5 will then present our summary and conclusions.

The classic way to deal with signal detection is the well-known spectral decomposition through Fourier analysis, whose aim is to find the series of frequencies f_l and the corresponding (possibly complex) coefficients x_l such that our periodic signal y(t_p) at the time point t_p can be written as:

$$\begin{eqnarray}&&y\left({t}_{p}\right)=\displaystyle \sum _{l=-{L}_{MAX}}^{{L}_{MAX}}\left[{x}_{l}\exp \left(i2\pi {f}_{l}{t}_{p}\right)+{\varepsilon }_{l,p}\left({f}_{l},{t}_{p}\right)\right].\end{eqnarray} \tag{ 1 }$$

In equation (1), ε is a possibly time-varying scalar quantity which includes the instrumental and random errors that could affect the measurement; ε generally depends on which sensors are used and on the frequency f_l: a more explicit expression for the error will be provided later; L_MAX is the total number of spectral components. In equation (1) we have used for simplicity the general conjugated variables time (=t) and frequency (=f) in the time domain, which could also represent conjugated variables in the spatial domain, i.e. toroidal (poloidal) angles and toroidal (poloidal) mode numbers in tokamaks.

Though immediate from a mathematical point of view, this method presents some rather serious drawbacks in its numerical implementation:

• the total number of frequencies 2L_MAX + 1 is not known a priori;
• the algorithm is non-linear with respect to the frequencies f_l.

A common alternative considers the theory of the Sparse Representation, in which a dictionary of frequencies, that in principle could be very large and redundant, is defined by the user according to some physical considerations. The problem, then, is to determine the coefficients x_l such that the solution vector [x] ⁶ is the sparsest possible, i.e. with the fewest non-zero components. In matrix notation, equation (1) then becomes:

$$\begin{eqnarray}&&{\bf{y}}=W{\bf{x}}+\varepsilon =\displaystyle \displaystyle \sum _{l\in {I}_{L}}{{\bf{W}}}_{l}{{\bf{x}}}_{l}+{\boldsymbol{\varepsilon }},\end{eqnarray} \tag{ 2 }$$

where [[W]] is the spectral window, defined as W_l,p = exp(if_lt_p), with [W_l] being its l-th column. Here [ ε ] is linked to the similarly-labelled, but scalar, quantity ε_l,p in equation (1) by [ ε _p] = Σ_l ε _l,p. For an evenly-distributed signal the matrix [[W]] satisfies the orthogonal condition [[W]^H[[W]] ∝ [[1]] (i.e. not necessarily orthonormal), while in the general case [[W]] is a Toeplitz matrix ⁷ .

The Sparse Representation approach, despite having a large number of unknowns (the coefficients {x}) offers the great advantage of being completely linear with respect to them, which makes numerical implementation much more efficient. The process of finding the sparsest solution of the under-determined linear system that best fits the data is based on a minimization procedure of a penalization criterion J(x) which, in its most general formulation, is written as follows:

$$\begin{eqnarray}&&J({\bf{x}})=\displaystyle \frac{1}{2}{\parallel {\bf{y}}-W{\bf{x}}\parallel }_{2}^{2}+\lambda R\left({\bf{x}}\right),\end{eqnarray} \tag{ 3a }$$

and [R(x)] is called the penalization function. Many different [R(x)] have been employed to address the minimization procedure. In principle, the one suitable to our problem is the L0 pseudo-norm ∣∣x∣∣₀ = Σ_i∈I0∣x_i∣₀ = Card{I₀} where I₀ = {i∣x_i ≠ 0}, which essentially corresponds to the number of non-zero components in the solution vector. However, since this norm lacks differentiability and convexity, it is computationally very costly to implement. Instead, it has been shown [12, 13] that under some assumptions the L0 norm is equivalent to the L1 norm, which has the advantage of being convex, so allowing only global minimizers. Equation (3a ) then reads:

$$\begin{eqnarray}&&J({\bf{x}})=\displaystyle \frac{1}{2}{\parallel {\bf{y}}-W{\bf{x}}\parallel }_{2}^{2}+\lambda {\parallel {\bf{x}}\parallel }_{1}.\end{eqnarray} \tag{ 3b }$$

The relaxing parameter λ plays a major role in our discussion since it is related to the freedom that we allow in the criterion to choose how many components to pick: the higher the value of λ, the fewer components we are allowing. Moreover, two important results can be proven:

• 1.  
  if the parameter λ is larger than a critical value λ_MAX, then the only solution is the zero vector;
• 2.  
  if the solution vector [x] has less than (2P + 1)/2 non-zero components, where P is the dimension of [y], then [x] is unique.

The value λ_MAX can be interpreted as the maximum amplitude of the residuals since, for any given value of it, the minimizer of the criterion given in equation (3b ) satisfies max(∣[[W]]^H([y] – [[W]]x_MIN)∣) ≤ λ. From the latter point, it immediately follows that we can uniquely specify λ_MAX = max([[W]]^H[y]), such that, therefore, for λ ≥ λ_MAX no solution can be obtained.

2.1. Definition of the input signal

To carry out our analysis of the SparSpec algorithm, an adequate signal vector S_IN must be built. Its adequacy lies in the fact that the signal should resemble, as accurately as possible, a real outcome of the data acquisition systems for magnetics fluctuations at a tokamak facility. Thus the vector of the measurement is a two-dimensional object (a K₁ × K₂ matrix) in which each entry represents the complex-valued signal ⁸ measured by the probe at the coordinate {ϕ_k1, θ_k2}, where {ϕ, θ} are the toroidal and poloidal angular coordinates, respectively. The index k₁ runs over k₁ = {1;2;...;K₁} and k₂ = {1;2;...;K₂}, with {K₁, K₂} both being of the order of 10 → 50. For simplicity, we will study the simplified case where the probes are placed along one spatial coordinate only (namely the toroidal coordinate ϕ), and thus the signal is a one-dimensional (1D) vector. The extension of the SparSpec algorithm to the joint analysis of a 2D {toroidal, poloidal} system is currently under development and will be presented later in a separate contribution.

If the signal was acquired without any error, it would be simply described by a sum over its Fourier components weighted by the Fourier coefficients, yielding a matrix [y] = [[W]]*[x] similar to that appearing in equation (2). Nevertheless, three sources of errors can be identified and have to be taken into account when building the signal:

• 1.  
  the instrumental error σ_k produced by the k-th sensor in making the measurement, which is determined by its precision, and the accuracy with which its position is known; considering also the electronic random noise, σ_k can be up to 15% of the total signal; note our convention whereby all σ's appearing in this (and its companion paper [20]) represent a relative error, and not a standard deviation as more commonly defined in statistical analyses;
• 2.  
  the delay δ_k (i.e. an equivalent phase shift) that the transmission line up to the data storage adds to the measurement made locally by each sensor;
• 3.  
  the error σ_l in the determination of the Fourier coefficients: considering the typical measurement accuracy, σ_l can also be up to 15% of the total signal.

As a result, a final expression for the signal vector at the time point t_p can be proposed as follows:

$$\begin{eqnarray}&&{S}_{0}\left({\phi }_{k}\right)=\displaystyle \sum _{l=-{L}_{MAX}}^{{L}_{MAX}}{W}_{l,k}\,\ast \,{x}_{l}.\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}&&{S}_{IN}\left({\phi }_{k}\right)=\displaystyle \sum _{l=-{L}_{MAX}}^{{L}_{MAX}}{W}_{l,k}\,\ast \,\left[{x}_{l}\exp \left(i\pi {\delta }_{k}\right)\left(1+{\sigma }_{l}\right)\right]\left(1+{\sigma }_{k}\right)\end{eqnarray} \tag{ 4b }$$

In equation (4) S₀(ϕ_k) is the nominal input signal, without measurement uncertainties and phase shifts, while S_IN(ϕ_k) is the signal as actually measured, which includes all sources of errors. We can now use this analytical form of S_IN(ϕ_k) to produce numerically a set of realistic input signals with which we can run our simulations.

It is now straightforward to see how SS-H2 can be modified to include the measurement errors: S_IN(ϕ_k) → S_IN(ϕ_k)/σ_TOT(ϕ_k), thus W_l,k → W_l,k/σ_TOT(ϕ_k), where σ_TOT(ϕ_k) = ∣S_IN(ϕ_k)-S₀(ϕ_k)∣/χ is the overall estimated measurement uncertainties from the simulations, i.e. the difference between the nominal and actual signal scaled by a normalization factor χ. This factor χ is needed in our case of simulated data as there is a certain degree of arbitrariness on how we estimate the measurement uncertainties. Setting χ = 1 correspond to using σ_TOT(ϕ_k)  =  ∣S_IN(ϕ_k)-S₀(ϕ_k)∣, namely the measurement error is exactly the difference between the nominal and the actually measured signal. Setting a higher value χ > 1 corresponds to reducing this error, while setting a lower value χ < 1 would correspond to increasing this error, which is unphysical ⁹ , thus this test case is not considered. We keep the same value of χ for all sensors ¹⁰ , and scan it as χ∈[1 → 10], as in actual (RT) experiments only good measurements are used, i.e. those for which the σ_k are the smallest and have similar values. We can then account for this σ_k in SparSpec by minimizing, instead of the cost function given in equation (3a ), the new cost function J_LS(x) = ([y] − [[W]]*[x])^T[[Γ]]⁻¹([y]-[[W]]*[x]) + λR([x]), where [[Γ]] is the diagonal matrix containing the σ_k ² on its diagonal ¹¹ .

2.2. The main features of the baseline SparSpec algorithm SS-H2

Because of the non-differentiability of the L1 norm at zero, a minimization procedure through the evaluation of gradients cannot be employed. Consequently, SparSpec uses an iterative Block Coordinate Descent (BCD) method that is based on subsequent scalar minimization steps for the complex-valued k-th component of the vector [x]. One of the features of this procedure is the low computational cost of each iterative step, which makes it suitable for applications where fast computations/results are needed. The BCD algorithm can be summarised as follows:

$$\begin{eqnarray*}&&{\rm{Iteration}}\,{\rm{t}}:{\rm{for}}\,{\rm{p}}=1:(2\times {{\rm{N}}}_{{\rm{MAX}}}+1):\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\bf{Update}}:{e}^{t}\left(:,p\right)=y-\displaystyle {\sum }_{l\lt p}W\left(:,l\right){z}_{l}^{t}-\displaystyle {\sum }_{l\gt p}W\left(:,l\right){z}_{l}^{t-1};\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\bf{Update}}:{z}_{p}^{t}={{\rm{\Psi }}}_{\lambda }^{s}\left[\overline{W}{\left(:,p\right)}^{\dagger }{e}_{p}^{t}\right].\end{eqnarray*}$$

It is possible to prove that the set of [x] minimizing the criterion J is equivalent to the condition:

$$\begin{eqnarray}&&\begin{array}{l}{x}_{k}^{\min }=\text{arg}\,{\min }_{{x}_{k}}J\left(\overline{x}\right)\Rightarrow {x}_{k}^{\min }={{\rm{\Psi }}}_{\lambda }^{s}\left({{w}_{k}}^{\dagger }{e}_{k}\right),\\ {{\rm{\Psi }}}_{\lambda }^{s}\left(z\right)=\left\{\begin{array}{c}\left(\rho -\lambda \right)\exp \left(i\varsigma \right)\,{\rm{if}}\,\rho \gt \lambda \\ 0\,{\rm{otherwise}}\end{array}\right..\end{array}\end{eqnarray} \tag{ 5 }$$

Here we have introduced the complex shrinkage function Ψ. In the above definitions, z = ρexp(iζ) is a complex number, [w_k ] is the k-th column of the matrix [[W]], and [e_k ] = y-Σ_l≠k[w_l ]x_l.

In words, the behaviour of the BCD algorithm can be explained as follows: at iteration t, it first computes the difference between the signal y and the Fourier Transform of the spectral components excluding x_k - particularly, x_l≠k computed at iteration t for l < k, and x_l≠k at iteration t-1 for l > k. If the absolute value of the inverse transform of the residual is lower than λ_REL × λ_MAX (up to a numerical tolerance), where λ_REL∈[0 → 1] is an user defined scaling factor for λ_MAX determining the penalized sensitivity of the solution to an increased number of output components (an higher λ_REL corresponds to less output components being sought), then x_k is set to zero; otherwise, it is shrunk by Ψ and used in the following iteration.

The only changes to be made in the algorithm to modify this new criterion J_LS(x) is to substitute [y] with [y/ σ _k ] (implying y_k/σ_k), and [[W]] with [[W/ σ _k ]] (i.e. W_k,l = W_k,l/σ_k). Then SS-V5 uses the weighting factor ss2 = Σ_k1/σ_TOT(ϕ_k)²∝χ² to obtain the {λ_MAX, ss2} penalized convergence in the successive iterations of the algorithm, see the shrink function described below. In SS-V5 λ_MAX,H2 becomes λ_MAX,V5 = max([[W/ σ _TOT ]]^H[y/ σ _TOT ]), hence λ_MAX,V5 = O(λ_MAX,H2/σ_TOT ²)»λ_MAX,H2.

An example of the implementation of the BCD algorithm in MATLAB is given in the following code, the shrinkH2 function (in bold in the following lines of code) being where the penalized convergence is tested.

t = 0;

lambda  =  max(abs(Wnorm' * y)) * lambdaREL;

z0 = zeros([2*Nmax + 1, 1]); z1 = zeros([2*Nmax + 1, 1]); e0 = zeros([K, 2*Nmax + 1]);

for p = 1:(2*Nmax + 1), e0(:,p) = y; end

e_t = zeros([K, 2*Nmax + 1]);

while t < 20

for p = 1:(2*Nmax + 1),

e_t(:,p) = y - lowSum(Wnorm, z1, p, K) - highSum(Wnorm, z0, p, K, Nmax);

z1(p) = shrinkH2(dot(Wnorm(:,p)', e_t(:,p)),lambda);

end

z0 = z1; t = t + 1;

end

function [phi] = shrinkH2(a,lambda);

if abs(a)<lambda, phi = 0; else, phi = (abs(a)-lambda) * exp(i * angle(a)); end

return;

function [a] = lowSum(W, v, k, K);

Slow = zeros([K,1]);

l = 1;

while l < k, Slow = Slow + W(:,l) * v(l); l = l + 1; end

a = Slow;

return;

function [a] = highSum(W, v, k, K, Nmax);

Shigh = zeros([K,1]);

for l = 1:(2 *Nmax+1-k), Shigh = Shigh + W(:,l + k) * v(l + k); end

a = Shigh;

return;

For SparSpec-V5, and in addition to using the error-normalized input quantities S_IN(ϕ_k)/σ_k, W_l,k/σ_k and λ_MAX,V5, the shrinkH2 function is simply modified as shrinkV5 and reads:

function [phi] = shrinkV5(a,lambda,ss2);

if abs(a) < lambda, phi = 0; else, phi = (abs(a)-lambda)/ss2 * exp(i * angle(a)); end

return;

Actually, a modified version of the BCD algorithm is used in the SparSpec code, which combines an Iterative Reweighted Least-Squares (IRLS) method: this ensures that J(x) does not grow over the different iterations, and if a component of the spectrum is zero at time t, it will also remain zero at time t + 1. Nevertheless, the BCD algorithm by itself as described above can be already applied quite successfully to spectral analysis. The BCD algorithm is the core of the SparSpec algorithm.

2.3. A few selected test cases for the SparSpec algorithms SS-H2 versus SS-V5ν0

For clarity of presentation and ease of reading, in this section we present a very limited number of selected test cases in order to illustrate the main features and the success of the SS-H2 and SS-V5ν0 algorithms. A more extensive presentation of the results obtained with SS-H2, both off-line and RT, and in diverse fields of science, can be found in the available literature, particularly [1, 12–18].

It is clear that multiple parameters need to be set to determine the accuracy and speed of SS-H2 versus SS-V5ν0. First, the input signal S_IN(ϕ_k) must be generated, and this is determined by:

• The number = N_SIG (labelled K₁ above) and the spatial distribution = {random, equi-spaced} of the sensors; as we mostly use SparSpec for the analysis of magnetic fluctuations in JET, TCV and ITER, and the number of magnetic sensors available for any such (toroidal, poloidal) analyses in these devices range from ~10 to ~50, which typically are not solely equi-spaced, for simplicity we set a fixed N_SIG = 35 and only use randomly distributed sensors.
• The number of non-zero input signal components =N_COMP and the value of L_MAX; magnetic fluctuations can be spatially analysed starting from a time series, hence multiple components are present, or at individual mode frequencies, hence typically only one component is present; thus we scan N_COMP = 1 → 5; in existing devices coherent fluctuations, for which discrete spatial periodicities (i.e. individual mode numbers) must be determined, usually have max(L_MAX) = O(5) for the toroidal and correspondingly max(L_MAX) = O(15) for the poloidal analyses, while for ITER we expect max(L_MAX) = O(15) and max(L_MAX) = O(35) for the toroidal and poloidal analyses, respectively; thus we set L_MAX = 35 to also cover the generic (toroidal, poloidal) ITER case; this choice of N_SIG = 35 with L_MAX = 35 is not at all facilitating the analysis, and is there to test at full the functionalities of SparSpec, as in a Nyquist sense one would want to have N_SIG≥2 × L_MAX.
• The randomly chosen values for the frequencies (=the spatial mode numbers) and amplitudes of these non-zero input components, thus determining the input set {x}.
• The phase shifts δ_k∈[0 → 1] × π and the errors {σ_k, σ_l}, all these being randomly distributed; for practical purposes, we set both {σ_k, σ_l}∈[0 → 1]×0.15 in relative units, i.e. with respect to the corresponding values of S_IN(ϕ_k) and x_l.
• The value of χ∈[1 → 10] to scan the sensitivity of SS-H2 (intrinsic, since ss2 does not appear in SS-H2) and SS-V5ν0 versus the overall measurement uncertainties.

Second, we need to define the analysis parameters required by the SparSpec algorithm, in addition to λ_REL∈[0 → 1]:

• The numerical set of atoms used to find the solutions: this has to be sparse, and it is practically limited to ≤∣FmaxS∣ only by RAM and CPU resources; previous analyses have indicated that FmaxS = [3θ → 5] × L_MAX would be typically sufficient for N_COMP = 5 and N_SIG = 35, , but since we do not have a problem with computing resources we set FmaxS = 10 × L_MAX.
• Since the solution set is sparse, it is possible that spurious, very-low amplitude elements are found due to the presence of noise in the input signal, and this can be eliminated using a cut-off value μ below which the amplitude of these spurious solutions is set to =0; this cut-off value is linked to the measurement uncertainties as μ = ss2/ξ, where ξ is user defined and can also be varied, but for simplicity ξ = 10 is always used.
• In principle n_IT = the total number of iterations allowed for the convergence procedure in the BCD algorithm and TR_IT = the relative (i.e. with respect to max(∣S_IN(ϕ_k)∣) tolerance threshold for checking the convergence conditions are user defined and can also be varied, but for simplicity n_IT = 1000 and TR_IT = 1 × 10^–5 are always used.
• Finally, as we randomize at each individual simulation the position of the sensors and the input signal, including the noise, for each value of N_COMP we need to run the analysis a very large number of times =N_ITER = O(10⁴), and calculate the RMS and STD value of our output results; with such a numerically intensive scheme, we believe that we should be avoiding all possible bias in the analysis of the performance of the SS-H2 versus SS-V5ν0 algorithms.

Our MATLAB routine randomly chooses K locations for the probes, the different spatial frequencies (modes) and their amplitude to avoid any bias in the analysis of the algorithm. The data that will be compared to assess and compare the accuracy and the speed of these two algorithms are the average elapsed time needed to converge, the relative error on the input versus output energy conservation, and the number of errors made in the result (hereafter, we will call them faults ¹² ). The rapidity of the computations has been measured through the MATLAB function timeit, which returns the clock time taken by the computer for the calculations. This quantity is then a relative, system-dependent, measure of how long it takes to perform one routine with respect to the other ¹³ .

Figures 1(a), (b) shows, as an illustrative example, two single test runs, always using N_COMP = 2 and χ = 1, but with λ_REL = 0.15 in figure 1(a) and λ_REL = 0.65 in figure 1(b). Note that the input signal (real and imaginary components) and the location of the sensors is different in the two runs presented in figures 1(a), (b), but obviously is the same for SS-H2 and SS-V5ν0 in the same run: the corresponding spectral window is only slightly different for SS-H2 and SS-V5, due to the weighting with respect to the σ_k.

Zoom In Zoom Out Reset image size

In general, the SS-H2 algorithm performs impressively well, and this sets quite a high target for the performance of the SS-V5ν0 algorithm. In both cases the two input components are detected with good accuracy by SS-H2 and SS-V5ν0. However, SS-H2 needs a ~6.15 × and a ~4.61 × longer run-time for the analysis, respectively. With λ_REL = 0.15 there are numerous spurious solutions with SS-V5ν0 and only one with SS-H2, while SS-V5ν0 has only one and SS-H2 none with λ_REL = 0.65. Apart from the very clear improvement in the required CPU run-time, it is difficult to choose which algorithm performs best from just these two test cases, as the differences between SS-H2 and SS-V5ν0 could simply be linked to the input signal and the location of the sensors being different. Therefore a statistical analysis is mandatory.

Finally, by inspecting the value of one of the outputs of the algorithms, i.e. the actual number of iterations (=N_IT) needed to obtain the solution, we find that in these two runs N_IT(SS-V5ν0) = [23 (figure 1(a)), 19 (figure 1(b))], while N_IT(SS-H2) = [130 (figure 1(a)), 82 (figure 1(b))]. This then correspondingly gives a ratio N_IT(SS-H2)/N_IT(SS-V5ν0) = [5.65 (figures 1(a)), 4.31 (figure 1(b))] which is close to the ratio in the required run-time. This result is similarly confirmed in all other analyses (for simulated and real data), and therefore we can state that the gain in computation time by SS-V5ν0 is essentially due to a much faster convergence, requiring significantly less iterations.

For illustrative purposes, figure 1(c) then shows a schematic layout of the distribution of the high-frequency (HF) magnetic probes in the JET tokamak, i.e. those later used for the analyses reported in this and its companion paper [20], together with a 3D folded view of the inner and the outer vessel structure. In the analysis reported here and in the companion paper [20], 11 magnetic sensors, not equi-spaced toroidally, were used with SS-V5 (all versions): these are the [T001, T002, H301(T003), H302, H303, H304, H305 (T004), T006, T007, T008, T009] sensors shown in this figure. Only 7 of these sensors were used in [21] with SS-H2: these were the [T001, H301(T003), H302, H303, H304, T007, T008] sensors. Not all input signals could be used with SS-H2, since they were characterised by somewhat different intrinsic errors and thus could not be assembled effectively as a unique input dataset.

When using simulated data, the main difficulty when testing the performance of SS-V5ν0 versus that of SS-H2 is setting the measurement uncertainties, through the value of the parameters {σ_k, σ_l, δ_k, χ}. While for actual experimental measurements these are (or should be ...) well known, in our case of simulated data there is a certain degree of arbitrariness since we could use different estimates of the errors {σ_k, σ_l} and phase shifts {δ_k} with which we construct S_IN(ϕ_k), and then correspondingly of σ_TOT = ∣S_IN(ϕ_k)-S₀(ϕ_k)∣/χ. This changes the penalization criterion for SS-V5ν0, first by the value of λ_MAX,V5 = O(λ_MAX,H2/σ_TOT ²)»λ_MAX,H2, and then by the complex shrinkage function Ψ, which is now scaled by ss2 = Σ_k1/σ_TOT(ϕ_k)²∝χ² in SS-V5ν0. The difficulty of this comparison is compounded by the fact the numerous additional quantities can (and should) be varied to generate the input signal, and then to perform the SparSpec analyses. It is however clear that if we fix {σ_k, σ_l}∈[0 → 1]×0.15 (in relative units), {δ_k}∈[0 → 1] × π, N_SIG = 35 (randomly spaced sensors), L_MAX = 35, N_COMP = 1 → 5, FmaxS = 10 × L_MAX, χ = 10, n_IT = 1000, TR_IT = 1 × 10^–5 as constant parameters that will be always used in our simulations, then we only need to scan χ∈[1 → 10] and λ_REL∈[0 → 1].

Therefore, in our view the best approach is to run a very large number of simulations N_ITER = 10'000 for each value of N_COMP = 1 → 5 and for each selected combination of {χ, λ_REL}, so as to remove possible bias due to the different sensors' geometry and input spectrum used in each individual run. For clarity of presentation, in the following we will only show first the results for a fixed χ = 1 scanning λ_REL = [0.15, 0.35, 0.65], thus progressively penalizing additional output components for the same overall measurement uncertainties, and then for a fixed λ_REL = 0.35 scanning χ = [1, 3, 10], thus evaluating the correctness of the penalization criterion as function of progressively smaller overall measurement uncertainties.

3.1. Comparison between SS-H2 and SS-V5ν0: λ_REL scan at fixed χ = 1

The first tuneable quantity for the comparison between SS-H2 and SS-V5ν0 is the value of the λ-parameter used in the penalization criterion, which appears in the shrinkH2 and shrinkV5 functions described above. The λ-parameter, controlled by λ_REL as λ = λ_REL×λ_MAX, keeping a fixed λ_MAX, directly appears in the L1-norm penalization criterion, see equation (3b ): an higher λ_REL corresponds to less output components being sought, as the λ-penalization for adding a new solution is higher.

The summary results of this λ_REL scan at fixed χ = 1 are shown in figures 2(a)–(c) for λ_REL = [0.15, 0.35, 0.65], respectively. For clarity of representation, the CPU run-time for SS-H2 has been rescaled by a factor /4 in these (and similarly, but by a different factor in some cases, in all the other) figures. The relative number of faults is similarly small for both SS-H2 and SS-V5ν0, always O(0.04) at most even for N_COMP = 5, with a similar scatter, and increase with N_COMP, as should be expected. The error on the input/output energy conservation, which is directly linked to the number of faults, is also similarly small for both SS-H2 and SS-V5ν0, O(0.10) at most even for N_COMP = 5.

Zoom In Zoom Out Reset image size

Note however that for λ_REL = [0.35, 0.65] the overall accuracy (number of faults and the error on the input/output energy conservation) is marginally better for SS-H2 than for SS-V5ν0. The most important result is the very significant reduction in the average CPU run-time, typically at least 4 × shorter for SS-V5ν0 than for SS-H2 for λ_REL = 0.15 and 2 × shorter for λ_REL = [0.35, 0.65]. For N_COMP = 1 → 5, SS-V5ν0 typically needs ~1.6 msec for N_COMP = 1 → 5 for λ_REL = [0.15, 0.35, 0.65]. Furthermore, while the average CPU run-time remains practically constant and has a small scatter for SS-V5ν0 for an increasing value of N_COMP, for SS-H2 it increases with N_COMP and has a progressively larger scatter. It is therefore clear from this λ_REL-scan at fixed χ = 1 that SS-V5ν0 out-performs SS-H2 in terms of speed for practically the same overall accuracy.

3.2. Comparison between SS-H2 and SS-V5ν0: χ scan at fixed λ_REL = 0.35

We now consider an actual test case from the JET tokamak, using data acquired with the Alfvén Eigenmodes Active Diagnostic (AEAD) system [22]. This test case is #79215, as extracted from [21], where the results were obtained using the SS-H2 version of the SparSpec algorithm. Looking at figure 3 in this previous work, we notice that in a number of instances toroidal mode numbers ∣n∣≥10 appear to be detected for stable, antenna-driven AEs, as resolved using SS-H2. Since this, and all similar discharges used for that work, are purely ohmic plasmas without sources of fast ions capable of driving (or supporting antenna-driven) such high-∣n∣ modes, and the antenna drive is also very small for these modes, it is rather difficult to believe that these are correct solutions. One good test is therefore to check whether these likely spurious solutions are still present if SS-V5ν0 is applied to the same set of data. However, as previously and repeatedly stated, we want to stress that we do NOT believe this to be the most relevant test for the validity of SS-V5ν0, simply because for actual data we do not know a-priori what the results should be, and our intuition might be erroneous. Only with simulated data, when we know what the results of our calculations must be, we can really be confident of the improvements brought about by any piece of software, in our case by SS-V5ν0 versus SS-H2.

Zoom In Zoom Out Reset image size

Figure 4 shows the results of this comparison, where we look at the occurrence of mode numbers associated to each individual frequency resonance for which damping rate measurements can be extracted, comparing the results of figure 3 in [21] obtained with SS-H2, and the results now obtained using SS-V5ν0. Details on how the analysis is performed can be found in [1, 21]. For the errors on the synchronously detected magnetic ∣δB_MEAS(ω,t)∣ and antenna current ∣I_ANT(ω,t)∣ measurements, producing the antenna-driven field ∣δB_ANT(ω,t)∣, used in SS-V5ν0, we consider those described in [23]. The relative uncertainty on the total normalized signal amplitude ∣δB_MEAS(ω,t)/δB_ANT(ω,t)∣ is due to propagation of the errors on the raw input magnetic signals and the antenna current, and their respective frequency calibration and precision of the synchronous detection, and is estimated to be around ~35%. In the analysis reported here 11 magnetic sensors (refer to figure 1(c): these are the [T001, T002, H301(T003), H302, H303, H304, H305 (T004), T006, T007, T008, T009] HF sensors), not equi-spaced toroidally, were used, but only 7 were used in [21] (again refer to figure 1(c): these were the [T001, H301(T003), H302, H303, H304, T007, T008] HF sensors): not all input signals could be used with SS-H2, since they were characterised by somewhat different intrinsic errors and thus could not be assembled effectively as a unique input dataset.

Zoom In Zoom Out Reset image size

The typical total antenna current in this experiment was ~15A. Therefore, in principle antenna-driven toroidal mode numbers up to ∣n∣≤30 can be detected (as illustrated by the labels in figure 3 in [21]), i.e. they have a nominal signal-to-noise (S/N) ratio S/N > 5 when considering only the random and bit-noise errors (in the absence of any antenna drive) on the synchronously detected signals. When considering the actual measurement errors, as SS-V5ν0 does, we expect that only modes up to ∣n∣ = 10 with ∣δB_MEAS(n)/I_ANT∣>10^–4mG/A can be correctly detected. All other higher-∣n∣ modes and lower mode amplitudes are estimated to be below the S/N > 3 detection threshold when including correctly the measurement uncertainties, as SS-V5ν0 does.

Looking now at the results shown in figure 4, practically all mode numbers ∣n∣>12 and/or for which the resulting amplitude ∣δB_MEAS(n)∣ < 10^–3mG are removed if the analysis is performed using SS-V5ν0, thus essentially confirming our initial intuition on spurious solutions. These spurious solutions are most evidently found during the large frequency sweep of the antenna (but not solely nor always either, see for instance between t  =  10.5 s and t = 11.2 s), i.e. when the AEAD system has not found a well-defined resonance to be tracked in RT (see [1] for more details on this). When tracking occurs, i.e. when the antenna performs small RT frequency sweep around a well-defined resonance, for instance between t = 2.2 s and t = 3.8 s, and then again between t = 6.2 s and t = 10.0 s, most solutions are found with both SS-H2 and SS-V5ν0. The typical CPU run-time needed to perform the toroidal decomposition for each detected frequency resonance is ~350 μs for SS-V5ν0, and it is around 4 × longer for SS-H2. Although just indicative, since not as powerful in our opinion as the extensive set of tests that can be performed using simulated data, this example on actual data from the JET tokamak, and the additional one reported in [20], confirms the significant improvements in terms of accuracy and speed brought about by SS-V5ν0 versus SS-H2.

This work starts from the analysis of the spatial periodicities of magnetic fluctuations in JET and TCV, where it was found that not all input signals could be used with the baseline version of the SparSpec algorithm, SS-H2, as they were characterised by somewhat different intrinsic errors and thus could not be assembled effectively as a unique input dataset. This then partially reduced the accuracy of the SS-H2 algorithm, as we indeed found that sometimes different results were obtained depending on which subset of input signals was used. In this work we have therefore compared, using mostly simulated data because we know what the answer must be, the performance of the SS-H2 versus the SS-V5ν0 algorithms. One test case is also run using actual magnetic data from the JET tokamak, and the results are in line with those obtained using simulated data. The development of SS-V5ν0 represents the first step towards the currently most advanced version of the SparSpec algorithm, SS-V5, which not only accounts for the measurement uncertainties but also uses a memory with relaxation scheme to further speed up the calculations. For clarity of presentation, the full version of the SS-V5 algorithm will be the subject of a separate, companion contribution [20].

When using an extensive set of simulated data, the main result from our work is that there is a significant improvement in the performance from SS-H2 to SS-V5ν0, both in terms of speed and overall accuracy. These improvements mainly come from correctly accounting for the measurement uncertainties. Additional sensors, with significantly different measurement uncertainties, can then be used by SS-V5ν0, and this also intrinsically improves the accuracy of the SparSpec algorithm. On simulated data, typically a factor at least 4× is gained in the average CPU run-time by SS-V5ν0 versus SS-H2 for an even slightly improved overall accuracy. While this is not an un-common result, it really shows the need of correctly introducing error tracking in the SparSpec code when implementing a memory with relaxation scheme. The results obtained on the test case run using actual data for the JET tokamak are in line with the simulations in terms of computational time, and confirms our initial intuition about the presence of spurious solutions in the dataset obtained with SS-H2.

Moreover, in our extensive tests using simulated data, we find that while for SS-V5ν0 the average CPU run-time remains practically constant, in the range from ~[370 ± 10]μs to ~[400 ± 10]μs and has a small scatter for an increasing value of input signal components, for SS-H2 it increases from ~[1.5 ± 0.1]msec for N_COMP = 1 to ~[2.0 ± 0.2]msec for N_COMP = 5, with a progressively larger scatter. The relative number of faults is slightly larger for an increasing error weighting in the penalization criterion, for both SS-V5ν0 and SS-H2. This occurs even if the related ss2 parameter does not appear in any form in SS-H2: SS-V5ν0 correctly weights the input signal, the spectral window and the penalization criterion with respect to the (simulated, in this case) measurement uncertainties, and this proves to be much faster, not only more accurate, than using an unweighted data set and thus a correspondingly unweighted penalization criterion. Finally, we have confirmed that the gain in speed between SS-V5ν0 and SS-H2 is indeed essentially due to a much faster convergence.

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 and 2019–2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. This work was also partly supported by the Swiss National Science Foundation. Finally, we would like to thank the Reviewer for his/her useful comments and suggestions to the first draft of this article.

The data generated and/or analysed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.