{'fraction_non_alphanumeric': 0.036612351639509155, 'fraction_numerical': 0.01460470730235365, 'mean_word_length': 4.404653185475103, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multi-Channel Three-Dimensional SOLA Inversion for Local Helioseismology

Inversions for local helioseismology are an important and necessary step for obtaining three-dimensional maps of various physical quantities in the solar interior. Frequently, the full inverse problems that one would like to solve prove intractable because of computational constraints. Due to the enormous seismic data sets that already exist and those forthcoming, this is a problem that needs to be addressed. To this end, we present a very efficient linear inversion algorithm for local helioseismology. It is based on a subtractive optimally localized averaging (SOLA) scheme in the Fourier domain, utilizing the horizontal-translation invariance of the sensitivity kernels. In Fourier space the problem decouples into many small problems, one for each horizontal wave vector. This multi-channel SOLA method is demonstrated for an example problem in time-distance helioseismology that is small enough to be solved both in real and Fourier space. We find that both approaches are successful in solving the inverse problem. However, the multi-channel SOLA algorithm is much faster and can easily be parallelized.

Essentially all inversions that have been employed in local helioseismology are linear. These inversions are based on the assumption of a linear relationship between perturbations to a reference model for the solar interior and the corresponding changes in the helioseismic measurements. The assumption of linearity is reasonable for inversions in the quiet Sun (e.g. Jackiewicz et al., 2007b;Couvidat and Birch, 2009). Within the context of linear inversions, the two main approaches are optimally localized averages (OLA: Backus and Gilbert, 1968) and regularized least squares (RLS: Paige and Saunders, 1982). In OLA methods, the goal is to produce spatially localized estimates of conditions in the solar interior while also controlling the associated random noise. In RLS methods, the goal is to produce a model of the solar interior that provides the best fit to the data under particular smoothness conditions.
The first three-dimensional (3D) inversions in local helioseismology were based on the RLS formalism (Kosovichev, 1996;Couvidat et al., 2005) and carried out using the LSQR algorithm (an iterative method, Paige and Saunders, 1982). RLS corresponds to Tikhonov regularization (Tikhonov, 1963) in the mathematical literature. This approach continues to be used extensively in time-distance helioseismology (e.g. Zhao, Kosovichev, and Duvall, 2001;Zhao and Kosovichev, 2004;Zhao et al., 2007). Jacobsen et al. (1999) introduced the Multi-Channel Deconvolution (MCD) approach to solving the RLS equations. In MCD, the (assumed) horizontal-translation invariance of the kernels is exploited to decouple the full RLS problem into a set of easily solvable problems, one for each horizontal wave vector. This method has been used by, e.g., Jensen, Jacobsen, and Christensen-Dalsgaard (1998); Jensen et al. (2001); Couvidat et al. (2004).
The 3D-OLA approach is computationally impractical for typical time-distance inversions due to the size of the matrices involved, as we will show in Section 4. An improved OLA variant, termed Subtractive OLA (SOLA) and introduced by Pijpers and Thompson (1992), allows one to perform fewer matrix inverse computations, yet does not reduce the sizes of the matrices. SOLA corresponds to what is known as the method of approximate inverse in the mathematical literature on regularization of inverse problems (Louis and Maass, 1990;Schuster, 2007). In problems where the kernel functions are separable as products of functions of horizontal position and functions of depth, it is possible to reduce the 3D-SOLA problem to a set of 2D-SOLA problems followed by 1D depth inversions (Jackiewicz et al., 2007a;Jackiewicz, Gizon, and Birch, 2008a). However, this separation is not always possible.
In this article we will show an efficient Fourier method for carrying out 3D-SOLA problems for local helioseismology. This method requires horizontaltranslation invariance of the reference model for the solar interior and is based on an MCD approach. We will show that, unlike direct solution of the SOLA equations, this method is computationally feasible. We will also demonstrate that the multichannel approach is many orders faster than the real-space method with example computations.

Setup of the Inverse Problem
The solar interior is filled with numerous scatterers such as flows, magnetic fields, hot and cold spots, density and pressure anomalies, etc. The scattering mechanism for each of these perturbations is physically different. Consider P such perturbations acting on the wave field. Not accounting for the magnetic field, the thermodynamic and flow perturbations would sum to P = 5 independent quantities. For example, the three components of flow velocity, sound speed, and first adiabatic exponent.
In time-distance helioseismology (Duvall et al., 1993), the measurements consist of travel times between points and concentric annuli or between points and quadrants. These travel-time measurements are performed for different choices of Fourier filters (e.g. Gizon and Birch, 2005;Jackiewicz, Gizon, and Birch, 2008b). Taken together, we assume that we have M such measurements, each denoted by the running index a, where 1 ≤ a ≤ M . In this article we are concerned with the local helioseismology of a small patch of the Sun near disk center, and we make the approximation that sphericity can be ignored. Thus we adopt a Cartesian coordinate system: where r is the horizontal position vector on the solar surface and z is height. The statement that a number of scatterers acts on the wave field to create small shifts in travel times [δτ a (r)] may be expressed as follows (e.g. Gizon and Birch, 2002): where δq β (r, z) represents the P perturbations in the various physical quantities that describe the solar interior, indexed by β. The sensitivity of a travel-time measurement [δτ a (r)] to a localised change [δq β (r , z)] is given by the traveltime sensitivity kernel [K a β (r − r, z)]. For point-to-annulus or point-to-quadrant measurements [δτ a (r)] the position vector r usually denotes the center of the annulus (although there is some freedom in this convention). Note that in Equation (2) we have explicitly assumed that the background solar model as well as the sensitivity kernels are invariant under horizontal translations. Sensitivity kernels result from a forward modeling under the single-scattering Born approximation (Gizon and Birch, 2002). The integral is taken over the volume of the Sun. The term n a (r) is the stochastic noise of the travel-time measurement [δτ a (r)] due to the forcing of waves by turbulent convection. The travel-time noise covariance matrix [Λ] has elements Λ ab (r i − r j ) = Cov n a (r i ), n b (r j ) . ( Details about the computation of Λ can be found in Gizon and Birch (2004). The general OLA inversion problem for time-distance helioseismology seeks to find an estimate of δq α at any chosen target position [r 0 ; z 0 ], given a set of δτ , K, and Λ. In other words, we are looking for a linear combination of the travel times such that is an estimate of δq α (r 0 ; z 0 ). The weights w α a (r i −r 0 ; z 0 ) are the unknowns of the problem. In the above equation, N = (2n + 1) 2 is the total number of horizontal position vectors [r i ]. Throughout the article we assume that the travel times are given on a uniform square grid with sampling h x in both horizontal directions and with n x = n y = 2n + 1 pixels on each side.

Subtractive Optimally Localised Averaging
Using the Equations (2) and (4) we have: We can rewrite Equation (5) as: where the functions K α β are averaging kernels given by where β is a running index between 1 and P that labels the physical quantities.
The term (6a) on the right-hand side of Equation (6) is what we are searching for, i.e. the quantity δq α convolved with the averaging kernel K α α . If the averaging kernel K α α (r, z; z 0 ) is well localized in the horizontal (around r = 0 and vertical (around z = z 0 ) directions, we will recover a smoothed estimate of δq α .
The term (6b) is the leakage from the other perturbations β = α into δq inv α . Ideally, one would like all averaging kernels K α β with β = α to be zero. The term (6c) represents the propagation of random noise from the travel times into the inverted δq inv α . The SOLA method consists of searching for the inversion weights w α a (r; z 0 ) so that the averaging kernels K α β resemble user-supplied target functions [T α β ], while keeping error magnification as small as desired. This can be achieved by minimising the cost function Equation (8) has two main components: The first term (8a), the "misfit", is a measure of how the averaging kernels [K α β ] match the target functions [T α β ]. To minimise the cross-talk, the target has only one non-vanishing component [T α α ], which we write as the product of a 2D Gaussian in the horizontal coordinates by a 1D function of the vertical coordinate. Thus we write where δ βα is the Kronecker δ-function. Typically, the function f peaks at a desired target depth z = z 0 . The constant C is taken so that the spatial integral of T α α is unity. The parameter s controls the width of the Gaussian. The second term (8b) of Equation (8) is proportional to the variance [σ 2 α ] of the random noise in q inv α due to the propagation travel-time noise: The trade-off parameter µ in Equation (8) is chosen to provide a satisfactory trade-off between the misfit and the noise; this choice is somewhat subjective.
The minimization is also subject to the constraints which ensure that the inverted quantity is normalised appropriately.

Linear System of Equations
The solution to the SOLA problem defined in Section 3 has been traditionally solved in real space for the one-and two-dimensional cases (e.g. Pijpers and Thompson, 1994). Below we write the problem for the 3D case. We recast the optimization problem of Equation (8) subject to the constraints (11) by minimising the function with respect to the vector of weights w α and a vector of Lagrange multipliers λ. There are M × N unknown weights w α a (r j ; z 0 ) and P unknown Lagrange multipliers λ β . The minimization consists of solving the following set of M × N + P equations: and Equations (13) and (14) imply where we define Equations (15) and (16) form a system of M × N + P linear equations. The weights w α a (r i ; z 0 ) can be obtained by matrix inversion, or some equivalent linear solver. Inspecting Equations (17), (18) and (19), one sees that the dimension of the matrix that is to be inverted is (M N + P ) × (M N + P ). An advantage of the SOLA method (see Section 3) compared to the OLA method is that the systems of equations for different target functions differ only in the right-hand sides and not in the matrix. Therefore the large matrix in Equation (17) has to be set up and inverted only once. For example, after inversion, we can infer the physical quantity at any depth using a new t α a (r i ; z 0 ). Let us now discuss the computational cost of the SOLA scheme for typical time-distance inversions. The computational costs depends very much on the size of the matrix A to be inverted. Typically, A may have over 10 14 elements (over a petabyte!), which do not fit in computer memory. In practice, we do not solve the full problem, but truncate the sums over j in Equations (15) and (16). We consider a restricted number of convolution "shifts" where h x is the horizontal sampling. The total number of shifts [N shifts = (2n shifts + 1) 2 ] must be much less than N so that the problem can now be solved. The minimum number of shifts that is acceptable depends on the size of the target function and the horizontal extent of the sensitivity kernels. In the example of Section 6, we have n shifts = 45. The smaller the number of shifts, the worse the approximation of the problem and the worse the localisation and the noise of the answer. If one takes the same number of parameters as in the (2+1)D flow inversion of Jackiewicz, Gizon, and Birch (2008a) then P = 3 (the three components of velocity), n shifts = 10, so that N shifts = 441, and M = 3 × 5 × 20 for three geometries (waves propagating in "outward minus inward", "West -East", and "North -South" directions), five ridges (f, p 1 , p 2 , p 3 , and p 4 modes), and twenty different radii. The problem would then require inverting a matrix of size ≈ 100 000 × 100 000 that occupies tens of gigabytes of memory. Furthermore, it is well-known that matrix inverse operations scale as O(N 3 ), where N = (N M + P ) is the length of the matrix on one side. Even without memory issues, this calculation becomes intractable very quickly.
Another issue that should be mentioned is the computation of the kerneloverlap integral in Equation (17). This computation is extremely expensive in real space. However it can be sped up very significantly by transforming to horizontal Fourier space, where it becomes a simple multiplication. Convolution operations are O(N 2 ) in the real space but reduce to O(N ln N ) when performed in Fourier space. In order to describe this step explicitly, let us define the discrete Fourier transform. Any function f and its Fourier transformf are related according tof where h k = 2π/(n x h x ) is the sampling in Fourier space. The horizontal wave vector k takes the discrete values k pq = ph kêx + qh kêy , where p and q are integers in the range [−n, n]. Using this definition of the Fourier transform, we obtain where we used the fact that the Fourier transform of K a β (−r) isK a * β (k) since K a β is real. In this form the kernel-overlap integral is computed much faster than in the Equation (17). The SOLA inversion examples presented later were computed using these equations. In order to avoid the edge effects resulting from the implicit periodicity assumed by the Fourier transform, we padded the sensitivity kernels and noise covariance matrices with zeros over a zone as wide as the size of the widest sensitivity kernel.

Solution in Fourier Space
In this section we fully exploit the horizontal-translation invariance of the sensitivity kernels and rewrite the entire problem in Fourier space. Using the definition of Equation (22), the Fourier transforms of Equations (15) and (16) are and where k takes the discrete values k pq = ph kêx + qh kêy , with p and q in the range [−n, n]. This set of equations can be written conveniently in matrix form where the vectort α (k) = t α 1 (k; z 0 )t α 2 (k; z 0 ) . . .t α M (k; z 0 ) T and the matrix A(k) = Ã ab (k) have the following elements where z bot and z top are the bottom and top heights of the computation box.
In Fourier space the problem decouples into many small problems, one for each horizontal wave vector k. These small problems are completely independent and therefore can be solved in a parallel fashion. The solutionw α a (k) is constructed for each k separately. By analogy to the RLS Multi-Channel Deconvolution (Jacobsen et al., 1999), we call the current approach multi-channel SOLA or MCD SOLA.
For each wave vector, the matrix to be inverted is much smaller than in the real-space case. For each wave vector k = 0, the matrix is of size M 2 . Taking the same parameters as in Section 4, the Fourier approach would only need 441 inversions of matrices of size 300 × 300. This would result in an increased speed by more that five orders of magnitude over the real-space method for this realistic example. Note that there is no need to truncate the problem anymore (n shifts = n).
We provide below expressions for the averaging kernel [Equation (7)] and the variance of the noise [Equation (10)] in terms of the Fourier transform of the weights: We emphasize that the averaging kernel is now computed as a matrix multiplication instead of a convolution.  The inferred solar property δq inv α at position (r, z 0 ) is where δτ a (k) is the Fourier transform of the travel-time maps.

Setup
We now show a rather simple example of a time-distance helioseismic inversion to demonstrate the Multichannel SOLA method and compare it to its realspace counterpart. For simplicity, we will only consider one mean (mn) point-topoint travel-time measurement with the distance between the two observation points fixed at ∆ = 15 Mm. For consistency with the measurement, we have computed a point-to-point Born-approximation sensitivity kernel according to Birch, Kosovichev, and Duvall (2004). This kernel gives the sensitivity of mean travel times to the sound speed perturbation [δc 2 /c 2 ]. No prior filtering has been done, i.e. the whole model power spectrum is used. We are not interested in specific types of kernels for this example problem, only a comparison between inversion methods and to prove that the MCD inversion works. This soundspeed kernel, denoted K mn c 2 according to the conventions in Section 2, is shown in Figure 1. This kernel has 91 × 91 elements in the horizontal direction and 80 elements in the vertical direction.
The second input quantity kept fixed for our example inversion is the target function, shown alongside the kernel in Figure 1. This 3D function has a Gaussian horizontal structure with a full width at half maximum of 20 Mm [see Equation (9)]. Since we are only working with one single sensitivity kernel (one ∆) in this example, it would be futile to attempt to obtain an averaging kernel peaked at some chosen z = z 0 , since the kernel itself possesses no such depth properties. Therefore, again for simplicity, we choose a depth profile of the target function by horizontally integrating the sensitivity kernel at each depth coordinate to obtain a one-dimensional curve according to f (z; z 0 ) = d 2 r K mn c 2 (r, z).
The 1D-curve f (z; z 0 ) is then combined with the horizontal Gaussian to construct the 3D target as in Equation (9). This choice for f (z; z 0 ) keeps the example as simple as possible. Note that the target in Figure 1 has a weak negative lobe beneath a depth of 10 Mm as a result of the depth profile of the sensitivity kernel.
The final input quantity to the inversion is the noise covariance matrix defined in Equation (3) and denoted in this case as Λ mn,mn . We compute the covariance from the model power spectrum according to Gizon and Birch (2004) and show the results in Figure 2. What this matrix tells us is how two mean point-topoint travel-time measurements are spatially correlated due to noise as the pairs of observation points are moved around with respect to each other. Note that for this case there is a significant correlation only when the measurements are made within about 10 Mm of each other.

Comparison of the Real-Space and Fourier-Space Solutions
We compute a set of real-space inversions and one Fourier inversion using the input quantities. The Fourier inversion is not computed in parallel for this example. For each inversion we generate a trade-off, or "L" curve (Hansen, 1998) by choosing ten values of the parameter µ [see Equation (28)]. The values of µ are chosen to span the space of misfit and noise. One trade-off curve is generated for the Fourier inversion, but several are generated for the real-space inversion, each corresponding to a different number of shifts employed. The possible n shifts range from 1 to 45, with 45 being the maximum due to the size of the kernel for this example (where n x = n y = 91). Since the MCD-SOLA inversion, in some sense, utilizes all possible shifts once, the real-space method with 45 shifts and the multi-channel method should agree.
In Figure 3 we show the results for these inversions. The top panel of Figure 3 shows the trade-off curves, with red lines indicating the real-space inversion for varying shifts indicated by the numbers at the bottom of the curves. The thick blue line is the MCD inversion trade-off curve. These curves are typically plotted as the square of the random noise level versus the misfit on a logarithmic scale. We see that for increasing n shifts the real-space inversion solution tends to the MCD solution. In fact, the L-curve for the 45-shift inversion falls on top of the one for the MCD inversion. Also shown in Figure 3 is a particular inversion weights w c 2 mn for each inversion, chosen from the first (topmost) point on each trade-off curve when µ = 0.01 Mm −3 . These are the points where the averaging kernel and target match best, i.e. smallest misfit. This is a reasonable choice since our main concern here is not the noise, which is still quite small anyway. The last two weights in Figure 3 demonstrate what the L-curves already suggest: the solutions of the two types of inversions are perfectly comparable when we take the maximum allowable number of shifts in the real-space method. The weights for inversions with a smaller number of shifts are quite ill-behaved due to edge effects, and in practice one actually never uses all possible shifts since it is computationally impractical to do so. This suggests that standard 3D-OLA inversions might have undesirable properties in the solution due to the necessary truncation of the problem. Cuts through all weights are shown in Figure 4, reinforcing this point when only a subset of shifts is used.
We recorded the computation time for each inversion. For the 45-shift case, the convolution matrix size is 8281 × 8281. The real-space inversion took two orders of magnitude longer to compute than the Fourier inversion (100 seconds compared to 1 second). This distinction only gets larger as the problem gets larger. Simply stated, the Fourier inversion takes a fraction of the time for small problems; for large problems, the real-space inversion is computationally intractable.  To show that the inversion does indeed work, in Figures 5 and 6 we provide comparisons between averaging kernel and target from the MCD solution. In the horizontal direction the agreement is quite acceptable, especially considering we have used only one input sensitivity kernel. Since the vertical profile of the input target function was constructed to match that of the sensitivity kernel, the good agreement there is expected.

Discussion and Conclusions
The example considered here is a very simple, toy inversion to demonstrate the usefulness of this new Fourier-based MCD-SOLA method. We have also experimented with a larger problem whereby we consider point-to-point measurements of various orientations of the observation points with respect to the x-axis. We input the same kernel as the one shown in this work, as well as horizontal rotations of it to match the measurements. Using the MCD method, we found that only five rotations, spaced evenly between 0 and 90 degrees and keeping ∆ fixed, are needed to find a very good averaging kernel that does not change with the addition of more rotations. The computation time for this inversion was about 1.5 seconds. Had we attempted to solve the same problem with the real-space OLA inversion, in addition to consuming 40 gigabytes of memory, it would have taken weeks to compute. We also solved our toy problem with several kernels of various distances [∆]. This allows one to obtain some resolving power in depth since the sensitivities differ. A target function was chosen as a 3D Gaussian peaked at a depth of 4 Mm beneath the surface. The MCD-SOLA inversion was able to find, as expected, an almost identical averaging kernel as the standard SOLA method. For a realistic application of the MCD-SOLA method with various target depths from 5 Mm to the surface, we refer the reader to the recent work ofŠvanda et al. (2011), who inverted for vector flows using synthetic travel-time observations as input.
In conclusion, in this article we have extended what was originally done for RLS inversions (Jacobsen et al., 1999) to a SOLA inversion. A toy example inversion problem was solved with this new approach to compare and contrast to the more standard real-space SOLA method. The example proved that the MCD-SOLA works completely satisfactorily while the real-space counterpart may be intractable for all but the smallest problems. In fact, we demonstrated that for a realistic helioseismic problem, the MCD-SOLA method can be orders of magnitude more computationally efficient than the corresponding real-space method. We focused here on applications to time-distance helioseismology, but this approach is completely generalizable to any local helioseismic method requiring inversions, such as ring diagram analysis and acoustic holography. With the vast amounts of seismic data from the Solar Dynamics Observatory (SDO), it is imperative to have efficient and consistent local helioseismic OLA inversion procedures for studying the solar interior.