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 *d*_{i} 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 Hz Hz.
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.

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 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 , with
elements *c*_{ij} 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 and imposing the constraint
, weighted by a parameter
, i.e. we minimize the
quadratic form *s* defined by

Therefore, the estimated solution is given by

The parameter 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 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 , for two values of :
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 ,
we need to define the goodness of the fit quantitatively.
To do this,
we compute the quadratic deviation as a function of the radial
point *r*_{i} and for a given :

where and *r*_{i}=(*r*_{1,i}+*r*_{2,i})/2.

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

**Figure 2:** as a function of spatial period
for *r*_{i}=0.3*R* (solid line) and for *r*_{i}=0.61*R*
(dashed lines)

We define the cut-off frequency as the value of
for which .
The spatial resolution at a point *r*_{i} 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.

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 *f*_{i} 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 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 at each point in the solution. We give now the general expression for the penalty function and then describe how our method compute .

In general, the minimization of relation (16 (click here))
can be expressed by a matrix , where *Q* is
defined by

in which case

Notice that since the are included in *Q*, and hence in *H*,
Eqs. (4 (click here)) and (5 (click here)) are replaced by

The parameters 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 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 *r*_{i} of the mesh has its own contribution to the eigenvalue
matrix . Let us introduce the matrixes defined by

where the elements *w*_{i}^{jk} of *W*_{i} are related to the elements *v*^{jk}
of *V* by the relation

It is straightforward to demonstrate that

The contribution of each point *r*_{i} is taken as the maximum
value of the corresponding diagonal of the matrix ,
and the weighting factors 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.

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 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 versus
the error propagation for different values of *c*. The curve has two
trends: for small values of *c*, 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:** 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 of computing time is dedicated to obtain the
spatial resolution, while 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*.

web@ed-phys.fr