next previous
Up: A new strategy

2. A new method: Optimal Mesh Distribution (OMD)

We have already pointed out the main problems when helioseismic data are inverted with standard techniques. The method we have developed, called Optimal Mesh Distribution (OMD), is a mixture of both the aforementioned techniques, although it is based mainly on RLS.

OMD tries to solve the main disadvantages of RLS by taking care of looking for the optimal mesh and defining the smoothing function by a deep analysis of the properties of the basis function of the problem. In particular, a different smoothness to the solution at each point of the mesh is used.

In order to illustrate the method developed here, we will take as di the rotational splittings parameterized as Clebsh-Gordan coefficients, so f(r) will be a measurement of the solar rotational rate (e.g. Ritzwoller & Lavely 1991). To test the method, a set of artificial data has been created. These consist in s=1 Clebsh-Gordan coefficients, obtained from the rotation curve shown in Fig. 6 (click here), for 1380 p-modes with l ranging from 1 to 150 and frequencies in the interval tex2html_wrap_inline1351Hz tex2html_wrap_inline1353Hz. Noise and error levels are taken from actual observations. We will call this set test 1.

In Sect. 3 we will apply the method to a more complex problem, where instead of one function f(r), we work with two, sound speed and density.

2.1. Determination of the mesh

The radial resolution that can be achieved in a given inverse problem depends on the mode data set used through the associated kernels. Since the kernels scan the solar interior in a significantly non-homogeneous way, the resolution strongly depends on r. In fact, one of the main problems when a solution is obtained, is that it can show oscillations due to an excessive number of points in the mesh. But then, it seems a good idea to decrease the number of points according to the spatial resolution; that is, we will have a non-equally spaced mesh point distribution.

The basic idea in obtaining the spatial resolution is as follows. Since the result is obtained by means of linear combinations of kernels, the solution at a given radial point can contain oscillatory patterns only with larger periods than those given for such a linear combination. Thus, in the procedure we will consider a given spatial oscillatory pattern (a sinusoidal function) and search for the combination of kernels that better fit this behaviour. If the fit is good at a given radial point, then the solution at this point is allowed for having such an oscillatory pattern. The smallest period for which a good fit is obtained at a given r gives the cut-off spatial frequency, and hence the spatial resolution. Finally, the mesh is built by taking only one radial point within the interval given by the spatial resolution.

To calculate the spatial resolution given by the kernels, a variant of the OLA method has been developed. As a measurement of the width of the averaging kernel, OLA uses a quadratic form. A variant of OLA, called SOLA (Pijpers & Thompson 1994), uses the quadratic deviation of the combination of kernels from a "target'' function, in the ideal case a Dirac delta function. In our case, the "target'' functions are sine waves of different frequencies; namely,
where tex2html_wrap_inline1361 is the spatial frequency. The amplitude a is obtained by assigning the area of the sine function as unity.

Let us call b the coefficients of the linear combination of kernels, the vector of dimension M we want to determine with this analysis. We want to minimize the quantity
where the vector P is defined by
The corresponding normal Gaussian equations are
where C is the kernel's covariance matrix of dimensions tex2html_wrap_inline1373, with elements cij given by
and y is the "target'' function vector,
Then the unknown vector b can be obtained by inverting relation (10 (click here)). Since C is a quasi-singular matrix, we obtain b as in the RLS technique, by keeping tex2html_wrap_inline1385 and imposing the constraint tex2html_wrap_inline1387, weighted by a parameter tex2html_wrap_inline1389, i.e. we minimize the quadratic form s defined by
Therefore, the estimated solution is given by
The parameter tex2html_wrap_inline1389 is introduced in order to avoid solutions with undesirable oscillatory components. Nevertheless, its value is very small, and there are a wide range of values of tex2html_wrap_inline1389 for which the fit to the target function is almost invariant and very good.

Let we consider test 1. As illustrative cases, Figs. 1 (click here)a and 1 (click here)b show the best fit of the kernel combinations tex2html_wrap_inline1397, for two values of tex2html_wrap_inline1399: 0.018 R and 0.11 R. The sine with spatial period of 0.11 R is well fitted by the combination of kernels at almost all the radial points. However for a period of 0.11 R, only points with r/R > 0.75 can be fitted properly.

Figure 1: a) Best fit to a sine with period 0.11 R. b) The same but for a period of 0.018 R. test 1 has been considered

Once b is obtained for a given spatial frequency tex2html_wrap_inline1361, we need to define the goodness of the fit quantitatively. To do this, we compute the quadratic deviation tex2html_wrap_inline1419 as a function of the radial point ri and for a given tex2html_wrap_inline1361:
where tex2html_wrap_inline1425 and ri=(r1,i+r2,i)/2.

In Fig. 2 (click here) is represented the value of tex2html_wrap_inline1429 as a function of the spatial period tex2html_wrap_inline1399, and for two points with ri=0.3R and ri=0.61R. The fit is good for values of tex2html_wrap_inline1361 for which tex2html_wrap_inline1439 is smaller than a given tolerance factor tex2html_wrap_inline1441. In particular, we have taken tex2html_wrap_inline1443, but as can be seen in Fig. 2 (click here), there is a sharp change in the behaviour of tex2html_wrap_inline1419, so there is a wide range of values for tex2html_wrap_inline1441 for which the result remains almost unchanged.

Figure 2: tex2html_wrap_inline1429 as a function of spatial period tex2html_wrap_inline1399 for ri=0.3R (solid line) and for ri=0.61R (dashed lines)

We define the cut-off frequency as the value of tex2html_wrap_inline1361 for which tex2html_wrap_inline1459. The spatial resolution at a point ri is defined as the inverse of the cut-off frequency. Figure 3 (click here) shows the spatial resolution obtained for test 1. Since only p modes are used, the resolution is worse at deeper points.

Figure 3: Spatial resolution of the solution for test 1

We can now build the inversion mesh by including only one radial point within the interval given by the spatial resolution. The resulting number of points is too low, in particular for test 1 we have N=89. Therefore the discretisation of relation (1 (click here)) could not be accurate enough. To avoid this problem, we have increased the actual spatial resolution by a factor of four, so that the number of points in the mesh grows in the same amount. With such a mesh, the discretization is very good, but some oscillations will appear in the solution.

2.2. Radial dependence of the smoothing function

As in the RLS method, to avoid oscillations in the solution, it is necessary to define and give a weight to the smoothing constraint. The other novelty of OMD is the way the smoothing function is defined directly from the spatial resolution analysis, and how it is weighted differently for each radial point.

Concerning the penalty function, we know that the mesh has four times the optimal number of points, so we will apply a fourth difference smoothing function H. To do this, we have to minimize the quadratic difference between the value of the solution at one point fi and the average value of the four surrounding points

Let we now consider the weighting of the penalty function. When the RLS technique was explained, a constant weighting parameter tex2html_wrap_inline1313 was defined, so the smoothing constraint is applied in the same way to all the solution. But then the weighting parameter is mainly determined by the result at the points that have a higher sensitivity to the kernels and, hence, the smoothing at other points can be dragged. To avoid this problem, we have applied a different weighting parameter tex2html_wrap_inline1483 at each point in the solution. We give now the general expression for the penalty function and then describe how our method compute tex2html_wrap_inline1483.

In general, the minimization of relation (16 (click here)) can be expressed by a matrix tex2html_wrap_inline1487, where Q is defined by
in which case
Notice that since the tex2html_wrap_inline1483 are included in Q, and hence in H, Eqs. (4 (click here)) and (5 (click here)) are replaced by


The parameters tex2html_wrap_inline1483 have been calculated following a method developed by Ruiz Cobo & del Toro-Iniesta (1994) of inverting the Stokes polarization line profiles. The matrix A, defined in Eq. (2 (click here)), can be written by using Spectral Value Decomposition (SVD) as
where tex2html_wrap_inline1501 is a diagonal matrix that contains the eigenvalues of A. U and V are orthonormal matrixes, with V containing the eigenvectors of A.

Each point ri of the mesh has its own contribution to the eigenvalue matrix tex2html_wrap_inline1501. Let us introduce the matrixes tex2html_wrap_inline1517 defined by
where the elements wijk of Wi are related to the elements vjk of V by the relation
It is straightforward to demonstrate that
The contribution of each point ri is taken as the maximum value of the corresponding diagonal of the matrix tex2html_wrap_inline1517, and the weighting factors tex2html_wrap_inline1483 are given by
where c is just a proportionality factor. In Fig. 4 (click here), the contributions of points for test 1 are presented. As expected, the contributions for outer points are higher than for inner points.

Figure 4: Point contributions for test 1

In summary, we have obtained a method, similar to RLS, that smoothes each point in the solution independently. The form of the smoothing function is given by the spatial resolution, which also gives the distribution of points in the inversion mesh. The only free variable of the method is the proportionality factor c, introduced in Eq. (23 (click here)), which will play the role of the trade-off parameter between error magnification and accuracy of the solution.

2.3. Solution for test 1

Here, we are going to see the performance of the method by inverting artificial data corresponding to test 1. To obtain the solution, we give different values to c, the only free variable of the method, and calculates tex2html_wrap_inline1439 and the error propagation for each solution. The error propagation of the solution is defined as a quadratic average over all the points of the mesh. Figure 5 (click here) shows tex2html_wrap_inline1543 versus the error propagation for different values of c. The curve has two trends: for small values of c, tex2html_wrap_inline1439 is almost invariant, while the error decreases very fast. For high c values the invariant behaviour is obtained for the error propagation (oversmooth solution). We have chosen as the optimal c the one where the two different trends concur (the "elbow''of the curve).

Figure 5: tex2html_wrap_inline1543 versus error propagation for different values of the trade-off parameter c

Figure 6 (click here) shows the solution obtained by our code for test 1. The result is quite good and fits the actual solution below the error level.

Figure 6: Result for test 1. The dashed line represent the actual solution and solid line shows the inverted solution with its error bars

OMD is a mixture of standard methods and is neither the fastest nor the slowest one. Over tex2html_wrap_inline1559 of computing time is dedicated to obtain the spatial resolution, while tex2html_wrap_inline1561 is used by the inversion analysis. The rest of the time is used to read files. The time ratio between RLS (the fastest), OMD and OLA is, approximately 1:12:40 for test 1.

next previous
Up: A new strategy

Copyright by the European Southern Observatory (ESO)