Improving Levenberg-Marquardt Algorithm Inversion Result Using Singular Value Decomposition

Inversion is a process to determine model parameters from data. In geophysics this process is very important because subsurface image is obtained from this process. There are many inversion algorithms that have been introduced and applied in geophysics problems; one of them is Levenberg-Marquardt (LM) algorithm. In this paper we will present one of LM algorithm application in one-dimensional magnetotelluric (MT) case. The LM algorithm used in this study is improved version of LM algorithm using singular value decomposition (SVD). The result from this algorithm is then compared with the algorithm without SVD in order to understand how much it has been improved. To simplify the comparison, simple synthetic model is used in this study. From this study, the new algorithm can improve the result of the original LM algorithm. In addition, SVD is allowing more parameter analysis to be done in its process. The algorithm created from this study is then used in our modeling program, called MAT1DMT.


Inversion
MT forward modeling problems can generally be written with d = g(m) (1) With d is data, which can be apparent resistivity values ( ) or phase (Φ), and g (m) is the forward operator that relates model parameters (m) with its response.In 1D case, models can be defined the by the layer resistivity and thickness of the subsurface.Jupp and Vozoff (1975) suggested the value of a model in logarithmic to avoid a negative model value in the inversion result.Inversion process is essentially a matching process of the calculated model response with the observational data.The implementation of this matching process is performed automatically using mathematical and statistical algorithms.In general, the inversion solution which can be applied to non-linear cases, including MT, is the Gauss Newton solution.
= + ∆ (2) With m is the solution models, Δm is the model parameters update, is the initial model, and J is the Jacobian matrix (J ij ) or models sensitivity matrix.

Levenberg-Marquardt Algorithm
Levenberg-Marquardt Algorithm is an algorithm that applies the minimization of model perturbation to the Gauss Newton solution.It can be done by minimizing the objective function The solution now becomes With the parameter λ is showing of the effect of model perturbation.Small λ value will make the equation ( 6) becomes equal to Gauss -Newton solution in equation (3).The parameter λ is initialized with a large initial value.If the solution after iteration has larger misfit than the previous iteration λ value is increased, otherwise it is reduced.To ensure the calculation stability the equation ( 6) is modified to equation ( 7).Jupp and Vozoff (1975) are suggesting another solution to improve the stability of the results.This solution is obtained by modifying the algorithm using Singular Value Decomposition (Lanczos, 1958).Singular Value Decomposition (SVD) is a method to divide a M x N sized matrix into three arbitrary matrices that has the relation : With U is M x N matrix which is eigenvector from data, V is N x N matrix eigenvector from parameter space, and S is N x N diagonal matrix which value is the square root from eigenvalue of J.
In that study, Jupp and Vozoff (1975) defining the damping parameter in LM algorithm using a N x N diagonal matrix P that has component value defined in equation ( 9).
With k i is ratio between the diagonal values from S matrix (s i ) with its first value (s 1 ), H is an integer which is defined H = 2 in the algorithm and is error level which is defined by the user and is determined with a high value and then decreased as the error of the solution model is reduced in every iteration, the determination of this value will be further explained in the discussion section.Using this damping matrix, the modified solution of LM algorithm using SVD is then become equation ( 9) With S + is N x N diagonal matrix that has 1/s i (i = 1,2,3, …, N) as its component values.

Scaling Factor
In our inversion algorithm, we are adding one more extra parameter to the calculation which is called the scaling factor.The idea of the scaling factor parameter is to scale the solution model into a new model which has the same average calculated response as the observed data.The scale factor is a scalar value that is obtained by dividing the calculated apparent resistivity of the solution model in an iteration and then multiply the factor to the solution model which can be expressed mathematically as equation ( 11) and equation ( 12).

SF =
( 1 1 ) With m s is the scaled model and m is the solution model in an iteration.
In the case of 1D MT, the observed apparent resistivity is greatly reflecting the earth subsurface resistivity model of the measurement point.Therefore, it is only natural to choose the solution resistivity model which has an average value that match the average observed apparent resistivity value and scaling factor is used to fulfill this condition.However, in order to explore more solutions, the scaling factor multiplication is not applied to every solution model in the inversion process; but instead, an extra constraint is added in regard of this matter.The constraint included in the algorithm is the phase fitting condition, which means that only solution model, where the calculated value is in a certain range of phase misfit, is getting multiplied by the scaling factor.The phase fit criterion was chosen because phase data is reflecting the changes of resistivity in between layer boundaries of the earth layers.
In this study, the phase misfit criterion value is set at minimum to make sure that the scaling process is only executed when the solution model has shown the correct pattern of resistivity.However, this value is increased when a certain number of error level increments do not produce the desired result.The purpose of this treatment is to fix the solution model so the desired result can be achieved, but this treatment may not be always give a good result since increased tolerance of phase misfit means that more model is scaled regardless of its resistivity pattern.

Algorithm Test
To simplify the validation in comparing the algorithms, we will show inversion using synthetic data.First example is an inversion form synthetic data obtained from forward modeling calculation of homogeneous layer earth and synthetic data obtained from forward modeling calculation of simple three layered earth.We will also use varying starting models to understand the effect of starting model on inversion result.In addition to that, we will display the quality of the result using variety number of error level parameter and our proposed determination of error level in our algorithm.
The quality of the solution model is measured by the misfit value which is the root mean square (RMS) error value between the calculated and observed values of apparent resistivity and phase, that are normalized to the observed data values and then it will be presented in the percentage value.

Results and Discussion
In this study, we will present three examples that will show the quality of both inversion algorithms.The first example is an inversion from synthetic data produced from a homogenous layer with resistivity value of 200Ωm using both original and modified Levenberg Marquardt algorithm.The starting model used in this process is a homogenous layer with resistivity value of 100Ωm.The next example is an inversion case of 3 layered earth models with resistivities ( ) and thicknesses ( ) as follows:  As it can be seen on figure 2 that both algorithm inversions are producing a good result, where the both algorithms are providing a solution that has small error in response relatively to the data.In addition to that, the original LM algorithm is almost forming the exact synthetic model used in the forward modeling process, but the modified algorithm inversion model is not forming the exact synthetic model.However this may not be an issue, since the deviation from the true resistivity is not significant, where the layers may still be interpreted as the same lithology.From the result shown in figure 3, it can be seen that the inversion process is producing various results using different starting model.Unlike the results produced from inversion using 100 Ωm uniform layers, the result from original algorithm has a significant difference, where a part of the solution model is becoming a very resistive layer.This 'bad' result is caused by the instability in inversion process, where a small change in data leads to large change in the solutions (Jupp and Vozoff, 1975).However, this problem is solved in the inversion result using modified algorithm with SVD.It can be seen in figure 3 that the model result from inversion using the modified algorithm is showing a consistent pattern of resistivity compared to the result before, when 100 Ωm uniform layers is used as starting model.Therefore, it can be said that the modified LM algorithm using SVD can produce a more stable results than the original LM algorithm, regardless of the starting model used.
Even though the modified algorithm has been able to create the solution that has less dependency to its starting model, the stability problem still persists.One of the causes is the determination error level parameter in the inversion process.As it can be seen on the Figure 4, the misfit has a tendency to significantly increased as the error level value is increased, especially when the value is above 0.15 (15% error value) and the minimum value of misfit is located there.However, there is also a lower limit of error level value when the inversion is not giving the desired result and instead in the case with starting model of 10 Ωm (red dotted line in the figure), the misfit value is increased at rhe lowest value of error level.Therefore the suggested determination of error value by Jupp and Vozoff (1975), which is determined to be a high value and then lowered if the misfit were reduced or increased if the misfit were increased, has been the right decision because it will ensure the right error level value, which is not too high or too low.In addition to that, we suggest that the initial value of error level to be in the range of 0.10 -0.20 because it is ensuring that the error level value not become too low or too high.Another inversion scheme available is to use various fixed error level value for the entire inversion and take the result with smallest misfit value as the solution model like the result illustrated on the Figure 4, although this scheme may not be recommended because it will consume more computation time.

Conclusion
Based on all results discussed in the previous section, it can be concluded that the original Levenberg Marquardt inversion result has a high dependence on starting model used in inversion.In some cases, this algorithm can produces good results, but in other cases it can also produces bad results, which can lead to misinterpretation due to instability that may happen in the inversion process.However, this problem can be solved by the modified algorithm using SVD.This algorithm is providing a more consistent results than the original ones, which is can be said that it is more reliable to be used in the inversion process.In addition to improve the result further we suggest the determination of error level initial value in the inversion to be in the range of 0.10 to 0.20 and if the inversion still does not give the desired result, a change of starting model may be required as its response can be too far from the solution model and the inversion solution is stuck on the local optima of the misfit.Also the use of multiple error level values in one inversion program can become one of the inversion scheme alternatives to see which error level value is giving the best result.
Ω , = 100 Ω , = 1000 , = 2000 .Starting model used for all inversion is 5 layered earth model with same resistivity value of 100 Ω .These two cases of inversion can be called as a good case, where the model results only show a small ambiguity effect in the model.Another example, which is shown on figure 3, is an example of a bad case, where a solution model with significant difference from the synthetic model can produce the almost similar responses as the synthetic model responses.The starting model used in the next example is 1000 Ωm.These examples can be seen on figure 1 to figure 3.

Figure 1 .
Figure 1.Inversion results from the original and modified Levenberg Marquardt algorithm using 100 Ωm homogenous resistivity layer as starting model.Part (a) shows the apparent resistivity curve, (b) shows the phase curve, and (c) shows the 1-D resistivity models obtained.Synthetic model is indicated by blue line and its responses are indicated by green dots for apparent resistivity and yellow triangle for phase.Models and responses obtained from inversion using the original algorithm are indicated by green line and the result of the modified algorithm is indicated by red line The inversion model results from both inversion algorithms, which are displayed on figure 1, are forming the same model as the synthetic model.These results are indicating that the iterations involved in the algorithm calculations are leading the starting model to a model solution that has the same apparent resistivity and phase response as the data, showing that the algorithm used is giving the desired result.As it can be seen on figure2that both algorithm inversions are producing a good result, where the both algorithms are providing a solution that has small error in response relatively to the data.In addition to that, the original LM algorithm is almost forming the exact synthetic model used in the forward modeling process, but the modified algorithm inversion model is not forming the exact synthetic model.However this may not be an issue, since the deviation from the true resistivity is not significant, where the layers may still be interpreted as the same lithology.

Figure 2 .
Figure 2. Inversion results from the original and modified Levenberg Marquardt algorithm using 100 Ωm uniform resistivity layers as starting model.Part (a) shows the apparent resistivity curve, (b) shows the phase curve, and (c) shows the 1-D resistivity models obtained.Synthetic model is indicated by blue line and its responses are indicated by green dots for apparent resistivity and yellow triangle for phase.Models and

Figure 4
is illustrating the effect of various fixed error level values used in the inversion process of the synthetic data shown at Figure 3 with various starting model.

Figure 4 .
Figure 4.The misfit value of the inversion at the convergence condition.Blue, red, green, and magenta dotted line are indicating the results from 200 Ωm, 10 Ωm, 100 Ωm, and 1000 Ωm starting model.Each dot in the line is representing 0.01 interval of error level value