Pseudo-Inverses

In this figure the next definitions are used A - projection operator, B - pseudo-inverse operator for the image parameters a,( ), C - empirical posterior restoration of the FDD function w(a, ), E - optimal estimator. The projection operator A is non-observable due to the Kalman criteria [10] which is the main singularity for this problem. This leads to use the two step estimation procedure. First, the pseudo-inverse operator B has to be found among the regularization techniques in the class of linear filters. In the second step the optimal estimation d (n) for the pseudo-inverse image parameters d,(n) has to be done in the presence of transformed noise j(n). 

Because the pseudo-inverse filter is chosen from the class of additive filters, the regularization can be done without taking into account the noise, (n). At the end of this procedure the noise is transformed to the output of the pseudo-inverse filter (long dashed lines on Fig. 1). The regularization criteria F(a,a) has to fulfill the next conditions (i) leading to an additive filter algorithm, (ii) having the asymptotic property a, —> a, for K,M... 

The adaptive estimation of the pseudo-inverse parameters a n) consists of the blocks C and E (Fig. 1) if the transformed noise ( ) has unknown properties. Bloek C performes the restoration of the posterior PDD function w a,n) from the data a (n) + (n). It includes methods and algorithms for the PDD function restoration from empirical data [8] which are based on empirical averaging. Beeause the noise is assumed to be a stationary process with zero mean value and the image parameters are constant, the PDD function w(a,n) converges, at least, to the real distribution. The posterior PDD funetion is used to built a back loop to block B and as a direct input for the estimator E. For the given estimation criteria f(a,d) an optimal estimation a (n) can be found from the expression... 

C CT] is known as the pseudo inverse of C. Since the product of a matrix and its inverse is the identity matrix, [C CT][C CT] disappears from the right-hand side of equation [32] leaving... 

Next, we post-multiply both sides of equation [60] by [AproJ ATproJ], the pseudo-inverse of ATproJ. 

Furthermore, the implementation of the Gauss-Newton method also incorporated the use of the pseudo-inverse method to avoid instabilities caused by the ill-conditioning of matrix A as discussed in Chapter 8. In reservoir simulation this may occur for example when a parameter zone is outside the drainage radius of a well and is therefore not observable from the well data. Most importantly, in order to realize substantial savings in computation time, the sequential computation of the sensitivity coefficients discussed in detail in Section 10.3.1 was implemented. Finally, the numerical integration procedure that was used was a fully implicit one to ensure stability and convergence over a wide range of parameter estimates. 

In the 3rd run the porosity of the ten zones was estimated by using an initial guess of 0.1. Finally, in the 4,h run the porosity of all fifteen zones was estimated by using the same initial guess (0.1) as above. In this case, matrix A was found to be extremely ill-conditioned and the pseudo-inverse option had to be used. 

As expected, the estimated values were found to be closer to the correct ones compared with the estimated values when the water-oil ratios are only matched. In the 2nd run, the horizontal permeabilities of layers 5 to 10 (6 zones) were estimated using the value of 200 md as initial guess. It was found necessary to use the pseudo-inverse option in this case to ensure convergence of the computations. The initial and converged profiles generated by the model are compared to the observed data in Figures 18.25a and 18.25b. 

In the same way, the parameters kij and hij are joined to form a unique parameter matrix H. With these definitions a linear problem may be written like that of equation 5. The matrix H can then be estimated either by direct pseudo-inversion or by PLS. It is worth noting... 

As stated earlier, Matlab s philosophy is to read everything as a matrix. Consequently, the basic operators for multiplication, right division, left division, power (, /,, A) automatically perform corresponding matrix operations (A will be introduced shortly in the context of square matrices, / and will be discussed later, in the context of linear regression and the calculation of a pseudo inverse, see The Pseudo-Inverse, p.117). 

In general, non-linear problems cannot be resolved explicitly, i.e. there is no equation that allows the computation of the result in a direct way. Usually such systems can be resolved numerically in an iterative process. In most instances, this is done via a truncated Taylor series expansion. This downgrades the problem to a linear one that can be resolved with a stroke of the brush or the Matlab / and commands see The Pseudo-Inverse (p.ll 7). 

Matlab is, of course, aware of the fundamental importance of the pseudoinverse and created its own notation for it. In Matlab we could write a=inv (F F) F y but it is numerically much more efficient to use the appropriate Matlab back-slash command as in a=F y. It is to be read from the right to the left as y divided by F, implying, of course, the multiplication of the left pseudo-inverse of F with y as given in equation (4.30). 

Figure 4-13. Schematic representation of the matrix equations involving multiplication of y by the left pseudo-inverse F+=(FtF) 1Ft.

For the computation of the pseudo-inverse, it is crucial that the vectors f j are not parallel, or more correctly, that they are linearly independent. Otherwise, the matrix FlF is singular and cannot be inverted. Matlab issues a warning. We can gain a certain level of understanding by adapting Figure 4-10 ... 

In a strictly mathematical sense this matrix is not singular but numerically it is rank deficient and has effectively a rank of only 4. Calculation of its pseudo-inverse consequently is impossible, or at least numerically unsafe. What can we do about that ... 

It is important to realise that the fitting of a polynomial to a series of (x,y)-data pairs is really independent of the actual values in the x-vector. In other words, in the above example, it does not matter whether the x-values are between 140 and 150 or between 1 and 11 or between -5 and +5. What matters is the relationship between the x-values and their revalues. As the data are equidistant, any equidistant vector with the right number of values can be used to generate F. In the improved SavGol function we chose the values from -n to +n. The second important observation is that consequently we do not have to recalculate F each time and more importantly, we do not have to recalculate its pseudo-inverse F+, which is computed outside the loop for all points. The result is a routine which is much faster and numerically much sounder ... 

It is important to realise that each column a j can independently be calculated from the appropriate column y j, irrespective of all the other wavelengths, using equation (4.49). The pseudo-inverse C+ is the same for all. The equivalent of (4.49) for all wavelengths can be written as... 

Referring back to Matlab, it is very important to use the correct slash operator or / for the left and right pseudo inverse. Applying the wrong one will invariably result in an error message or worse, in a potentially undetected error. 

The pseudo-inverse for the calculation of the shift vector in equation (4.67) has been computed traditionally as J+= (J Jp1 J. Adding a certain number, the Marquardt parameter mp, to the diagonal elements of the square matrix J J prior to its inversion, has two consequences (a) it shortens the shift vector 8p and (b) it turns its direction towards steepest descent. The larger the Marquardt parameter, the larger is the effect. In matrix formulation, we can write ... 

J is the derivative of a matrix with respect to a vector. What is the structure of such an object and more disturbingly, what is its pseudo-inverse ... 

Note, that for the individual analyses, one of the two component spectra of species A or B needs to be set as colourless, i.e. non absorbing (see Known Spectra, Uncoloured Species, p.175), as the matrices Ci and C2, containing the concentration profiles, are rank deficient for this kinetic model and hence their pseudo-inverse C+ is not defined. 

For a linear fitting exercise, e.g. the calculation of the emission spectra A, we assume to know the lifetimes t and hence the matrix Csim, which we used for the generation of the measurement. The linear regression has to be performed individually at each wavelength. This is due to the fact that at each wavelength Xj the appropriate vector (Ty j is different and each weighted matrix Cw and its pseudo-inverse, needs to be computed independently. There is no equivalent of the elegant A=C Y notation. 

The spreadsheet in Figure 4-62 is heavily matrix based (see Chapter 2, for an introduction to basic matrix functions in Excel). It is the only way to keep the structure reasonably simple. The matrix C in cells A21 C31 is computed in the usual way, see equation (4.63) the parameters required to compute the concentration matrix are in cells Q4 S4, they include the initial concentration for species A and the two rate constants that are to be fitted. In cells E 16 018 the computation of the best absorptivity matrix A for any given concentration matrix C, is done as a matrix equation, as demonstrated in The Pseudo-Inverse in Excel (p.146). Similarly the matrix Ycaic in cells E21 031 is written as the matrix product CA. Even the calculation of the square sum of the residuals in cell R7 is written in a compact way, using the Excel function SUMXMY2, especially designed for this purpose. We refer to... 

The computation of the linear parameters b is easy, an orthonormal matrix is equal to its transposed as the pseudo-inverse of... 

The above equations are valid for orthonormal sets of basis vectors. They can be written in very similar ways for general non-orthogonal bases (e.g. F in Figure 4-12). The only difference is the computation of the pseudo-inverse, which can be numerically demanding, but is trivial for orthonormal bases. 

The advantage is that there is no pseudo-inverse to be calculated in this way. The computation of Tu, which comprises linear parameters is easier than usual as U is an orthonormal matrix, U+=Ul. As mentioned before, in equation (5.29), it is advantageous to compute the residuals as Ru = C - U(Ul C) it is considerably faster. 

---