A General Algorithm to Solve Linear and Nonlinear Inverse Problems

Um algoritmo geral para resolver problemas inversos lineares e não lineares, baseado em redes neurais recursivas, é discutido neste trabalho. O procedimento será aplicado a problemas físico-químicos modelados por equações integrais, diferenciais e de autovalor. As aplicações são discutidas em espectroscopia de aniquilação de pósitrons, cinética química e espectroscopia vibracional. O método é robusto com relação a erros nas condições iniciais ou em dados experimentais. A presente abordagem é simples, numericamente estável e tem uma grande aplicabilidade.


Introduction
Extracting chemical and physical information from experimental data is an important problem in science.This procedure, which characterizes an inverse problem, 1 is usually ill-posed in the sense that one of three conditions: (i) existence, (ii) uniqueness and (iii) continuity is not satisfied.When dealing with real data one has to handle ill-posed problems and this can be better clarified if a particular case is considered.For example, in an integral formulation of a physical problem, to be discussed along the paper, the solution of Kf=g has to be found.In such cases the matrix K is ill-conditioned, and consequently, its inverse will have very large values.Experimental errors will be amplified by the inversion procedure and condition (iii) is not fulfilled.Multiple solutions and existence conditions are also understood by analyzing the properties of the matrix K, as discussed in reference. 1For every inverse problem there is a direct problem that is considerably easier to solve.For example, in a chemical kinetics process, if initial conditions, rate constants and kinetic law are given, the concentration of the species can be calculated along the time.For a given set of concentration, measured at different observation times, one may ask what are the rate constants. 2Obtaining force fields from experimental data is another important inverse problem. 3 Calculation of vibrational frequencies, from a given set of force constants, will define the direct problem.On the other hand, force constants retrieved from the experimental frequencies will characterize the inverse problem.
Numerical methods to handle linear inverse problems are well established.For example, techniques such as Tikhonov regularization, 1 decomposition into subspaces 4 and dynamical neural network approach 5,6 have been used.A comparison between Tikhonov regularization and decomposition into subspaces has been discussed, 7 but the neural network framework has been shown to be more stable, simple and fast when applied to inverse problems. 8he dynamical neural network methodology was first suggested in reference 5 and applied to linear inversion problems by Vemuri and Jang. 9Several linear ill-posed problems were handled by this recurrent neural network technique, such as in: linear thermodynamics problem, 10 positron lifetime spectroscopy, 11 magnetic resonance multiple sclerosis diagnostic. 12All these applications have in common the integral form of the operator representing the direct problem.The original algorithm 9 has been extended to solve nonlinear integral problems 8 and to handle eigenvalue force field inverse problems. 13n the present work a general methodology to deal with inverse problems is presented.This algorithm will be able to solve linear integral inverse problem, nonlinear integral inverse problem, inverse problems in dynamic systems and identification problems.Some models will be used to exemplify the present approach.As a first example, inversion of positron annihilation lifetime spectrum represented by a Fredholm integral equation, will be taken.In a second example, chemical kinetics will be treated to illustrate inverse problems in differential form.The rate constants will be obtained from experimental concentration data.The final example will treat the eigenvalue inverse problem in which force constants are obtained from experimental frequencies.

Theoretical method
A Hopfield neural network is constructed by a recurrent layer, consisting of neurons fully connected.For example, in a Hopfield network formed by ten neurons, each neuron will have ten inputs, one for its own previous value and nine for the remaining ones.In this way, the output of a neuron u is a function of the input information, which is converted to another state by a activation function, f(u).The activation function propagates the required information and this occurs by following an energy-descent pathway.Therefore, the Hopfield neural network can be applied to solve optimization problems if an energy function is defined, which for practical applications, can be the error function for the physical problem.The learning property of the neural network is satisfied if f(u) is an increasing function of u, a condition to be fulfilled by functions of the kind tanh(u). 5 Both of these activation functions are used in the examples studied here.It is more appropriate to use the activation function results assume only positive values, as in the positron annihilation lifetime spectrum.This restriction will prevent the inverted function to have negative values.
For the same reason, this activation function was also used to invert data in inverse kinetics problem.In the force field inverse problem, the function f(u) = tanh(u) was used, for the inverted results represent force constants which can be negative.A more clarifying discussion about the recurrent neural network dynamics can be made by considering the error between measured and calculated properties.If the experimental and calculated properties are denoted, respectively by P EXP , P CAL and with o j = P j CAL -P j EXP the total error can be defined as, for m experimental data.Derivation with respect to the learning time, τ, gives, in which n is the quantity of neurons involved in the process.This is also equal to the number of variables to be inverted from experimental data.The total error will always decrease with respect to the learning time if two steps are taken.First, the condition Secondly, the relation is applied, which will imply in, Therefore, the present algorithm will have the property of always decreasing the error with respect to the learning time.Due to the nature of the ill-posed problems, multiply solutions will always be presented.One has to choose, based on physical and chemical grounds, among the several possibilities.Nevertheless, the present approach has the property of always decreasing the error between the inverted and experimental data.Therefore, among the infinite solutions, the one obtained in the present work is the solution which best reproduces the experimental property.The properties P EXP and P CAL are not necessarily the system variables, but can be a function of them.For example, recovering potential energy function from second virial coefficient will require solving an integral equation over the coordinates. 8In chemical kinetics, the measured property can be represented by a combination of the concentration, as in the absorbance measurement. 2For the force field inverse problem the property is the eigenvalue of a matrix.
The differential equations (3) were integrated by using a fourth order Runge-Kutta method with variable stepsize, 14 until the learning process, which defines the stopping condition, has finished.For the three examples to be discussed in this work, this stopping condition occurs if the error function reaches the minimum value, i.e., if there is no further decreasing in the network energy.Retrieval of the required quantities, that is, solution of the inverse problem, is obtained according to the following algorithm: (i) an initial guess to the required quantities is made and used to propagate the differential equations; (ii) since error will decrease, the time derivative of the neuron state will reach a constant value with zero derivative.This will be represent by also the final learning time; (iii) the converged values of the neuron states are activated to obtain the desired result.
Several inversion problems can be solved by these steps.Artificial errors were introduced into the simulated and experimental data to test the robustness of the algorithm.These errors were randomly generated with negative and positive values having the same weight.Activating the states u will have an effect on the quality of the inverted results.For example, if the function to be calculated happens to be only positively, this can be easily incorporated into the present methodology by imposing

Positron lifetime spectroscopy-inversion in integral equations
The inverse problem to determine s(y) in the Fredholm integral equation of first kind, 1 k x, y s y dy g x , from a given k(x,y) and g(x), will be tested for the present algorithm.As a representative example, the inversion of positron annihilation lifetime spectrum 6 will be discussed.This problem was previously discussed by using a singular value decomposition method, in a previous work. 11The density function of the positron annihilation rate spectrum, was taken as experimental result for the present analysis.
In equation ( 6), a 1 =0.42250, α 1 =11.6951,σ 1 =0.57710, a 2 =1.57270, α 2 =12.0688 and σ 2 =0.93530, parameters for lysozyme in water system. 10The theoretical data for this problem, calculated from max 0 ( ) ( ) , were generated by using a rectangular basis, n=128.Therefore, the error function in this example, will be Inverted results, with initial condition equal to zero and random errors of ± 0.3% added to the experimental data (the experimental error in this kind of measurement is about ± 0.3%) are compared with exact results in Figure 1.These results show the stability of the algorithm for this problem.Moreover, the algorithm discussed here overcomes the previous framework applied to linear inverse problems, 6 in which it was necessary to construct two matrices, derived from the available data and the linear transformation between input and output.For the algorithm proposed here, these matrices are unnecessary.This is an improvement in the inverse problem methodology since linear and nonlinear ill-posed problems can be treated on equal grounds.

Chemical kinetics-inversion in differential equations
The analysis of the fission and formation of C-C bonds in chemical reactions are often used as model systems in molecular kinetics.Experimental data consistent with the recombination-dissociation rate constants of several reactions are available in the literature.For example, trifluoromethyl radicals recombination, CF 3 +CF 3 → C 2 F 6 , has been studied experimentally 15 and theoretically to test other inverse algorithms. 16The ruby laser radiation has been used to produce photolysis of CF 3 NO. 15In this technique, which is effective to produce CF 3 radicals, the CF 3 and NO species are obtained, due to the photolyzing pulse of the laser.The CF 3 NO concentration was monitored by absorption in the visible spectral region and the reaction cell filled with a mixture of CF 3 NO, NO and a buffer gas.Therefore, concentration measurements of the CF 3 NO component were used to determine the rate constant k 2 for this recombination process.With experimental errors in the range 2% to 7% this rate constant was established to be (3.9 ± 1.3) × 10 -12 cm 3 s -1 .Using the notation x 1 =CF 3 , , the mechanism for the recombination process is described by the set of six differential equations, Thus, the inversion problem to be solved here will be: for given values of the rate constants k 1 , k 3 , k 4 and values of x 5 EXP along the time, what will be the rate constant k 2 ?The rate constants 16 will be taken as: k 1 = 1.3×10 -11 cm 3 s -1 , k 3 = 1.4 × 10 -13 cm 3 s -1 and k 4 = 2.0×10 -11 cm 3 s -1 .The CF 3 NO concentration was simulated in the interval 0.2 to1.0s, with step of 0.2s with 3.9× 10 -11 cm 3 s -2 for the k 2 rate constant.In this case, integration of equation (3) will provide the temporal evolution of the neurons which represent the inverted rate constant k 2 .With this rate constant it is possible to calculate the concentration x 5 CAL in an iterative process to minimize the error function 1 2 , according to equation ( 5).The experimental error was taken to be, at maximum, of 7%. 15hese results were taken as experimental data from which inversion is to be performed.Integration of the differential equations ( 7) was performed with x 1 (0) = x 2 (0) = 1×10 11 molecules cm -3 and x 3 (0) = x 4 (0) = x 5 (0) = x 6 (0) = 0 as initial conditions.Retrieved results are presented in Table 1.Even with the largest experimental error tested in the concentrations, the inverted rate constants are in agreement with the experimental value.The less favorable inverted rate constant, 4.40×10 -12 cm 3 s -1 occurs at an experimental random error of 7% in the concentrations.Also in this situation, the inverted rate constant is within of the experimental value, (3.9±1.3)×10 -1 cm 3 s -1 .The algorithm was very stable in relation to different initial conditions.
Even with an initial guess two times the correct value, rate constants shown in Table 1 are obtained.

Vibrational spectroscopy-inversion in eigenvalue problems
The eigenvalue inverse problem is based on the general equation AU = UΛ, in which A is a matrix related to the desired information, U are the eigenvectors and Λ is the diagonal matrix of eigenvalues.Although the present algorithm has been applied to this kind of problem, 13 it will also be applied here to another system.This additional example will exemplify the general applicability of the method for this problem.There are three main equations to formulate the problem: (i) the eigenvalue equation proposed by Wilson et.al., 17 GFL = LΛ in which G is a matrix related to mass and geometry of the molecule, F the matrix representing the forces constants, Λ a diagonal matrix with vibrational frequencies and L the normalized vibrational mode amplitudes matrix; (ii) equation ( 1) with the error function and (iii) equation (5), that guarantees a minimum error function, providing the calculation of reliable force constants.Therefore, the G and Λ matrices are known in formulating the problem.On the other hand, the F matrix is to be established.The error function in this problem is are the eigenvalues of the matrix GF calculated with the inverted force constants and Λ EXP are the experimental frequencies.The fluoroform molecule in its equilibrium configuration belongs to symmetry point group C 3v and this is considered to obtain the matrix G.The force constant regularized matrix, 3 with deviations up to 30%, was used as priori physical information, defining the initial condition for the network.The optimized results are presented in Table 2 and are in agreement with the regularized force field found in reference. 3The errors in the calculated frequencies, with respect to the experimental data, can not be attributed to the method itself, but to the harmonic force field model 17 considered.This specific model was used to carry on a comparison with previous results in the literature. 3oreover, for the experimental data used, the errors in the calculated vibrational frequencies are in agreement with the experimental error, as presented in Table 3. Table 4 shows the inverted results are stable and convergent for large experimental errors.Experimental deviations in the fundamental frequencies 18 which are about ±0.2cm -1 , corresponding to 0.03%, do not perturb the inverted results.The sensitivity of the inverted results with respect to different initial guesses is presented in Table 5. Calculation of retrieved force field presents low average error, up to 20% of deviation in relation to the exact initial guess.

Conclusions
A general technique to handle inverse problems is presented in this work.The algorithm can be applied to a problem in which there is not an explicit function to represent the data.For example, in the inverse chemical kinetics problem an analytical expression for concentrations was not available, but the inverted problem was still possible to be carried out.Experimental errors were also considered in the present approach.For the three examples studied, the corresponding experimental errors were taking into consideration.Errors in the initial conditions were used to test the robustness of the methodology.In both cases, either with experimental errors or initial conditions deviations, the present approach was stable, giving accurate physical results.Also, the method is simple and represents an extension of the previous used algorithm, 8 incorporating linear and non linear cases.The same approach can equally be applied to an integral, differential or eigenvalue inverse problem, for a general case of the measured property.
uniformly distributed were used to simulate the experimental error.

Table 1 .
Recovered rate constant with k 2 = 0 as initial guess

Table 5 .
Fluoroform force constants.Sensitivity to errors in initial guess

Table 2 .
Fluoroform force constants for the A 1 symmetry

Table 3 .
Fluoroform vibrational frequencies for the A 1 symmetry

Table 4 .
Range of errors in inverted fluoroform force constants