DCIP3D manual: 
1.6 Compression of Sensitivity matrix

Both DC resistivity and IP inversions involve a large dense sensitivity matrix. Since we are faced with 3D problems, this matrix can become very costly to store and to apply to a vector during the inversion. We adopt our wavelet compression method developed in the 3D magnetic inversion to DCIP3D programs. In this method we form a sparse representation of the sensitivity matrix using a wavelet transform based on compactly supported, orthonormal wavelets. For more details, the users are referred to Li and Oldenburg (1997). In the following, we give a brief description of the method. Understanding of this material is necessary for anyone wanting to use the DCIP3D library.

Each row of the sensitivity matrix in a 3D inversion can be treated as a 3D image and a 3D wavelet transform can be applied to it. By the properties of the wavelet transform, most transform coefficients are nearly or identically zero. When the coefficientswith small magnitude are discarded (the process of thresholding), the remaining coefficientsstill contain much of the necessary information to reconstruct the sensitivity accurately. These retained coefficientsform a sparse representation of the sensitivity in the wavelet domain. The need to store only these large coefficientsmeans that the memory requirement is reduced. Further, the multiplication of the sensitivity with a vector can be carried out by a sparse multiplication in the wavelet domain. This greatly reduces the CPU time. Since the matrix-vector multiplication constitutes the core computation of the inversion, the CPU time for the inverse solution is reduced accordingly. The use of this approach increases the size of solvable problems by nearly two orders of magnitude.

Let be the sensitivity matrix, and be the symbolic matrix-representation of the 3D wavelet transform. Then applying the transform to each row of J and forming a new matrix consisting of rows of transformed sensitivity is equivalent to the following operation,

     (26)

where is called the transformed matrix. The thresholding is applied to individual rows of by the following rule to form the sparse representation ^s,

     (27)

where _i is the threshold level, and _ijand ^s_ijare the elements of and ^s, respectively.

The threshold levels _i are determined according to the allowable error of the reconstructed sensitivity, which is measured by r_i(_i), which is the ratio of the norm of the error in each row to the norm of that row. It can be evaluated directly in the wavelet domain by the following expression:

     (28)

Here the numerator is the norm of the discarded coefficients. For each row we could choose _i such that r_i(_i)=r*, where r* is the prescribed reconstruction accuracy. However, this is a costly process. Instead, we choose a row, i₀, corresponding to a datum nearest the mesh centre, to be a representative row, and calculate the threshold level _i0 using eq.(28) and a prescribed value of r*. This threshold is then used to define a relative threshold The absolute threshold level for each row is obtained by

     (29)

This compression procedure is implemented in three programs in the DCIP3D library. They are DCINV3D, DCAIM3D, and IPSEN3D. The user needs to specify the relative error r* and the program will determine the relative threshold level . Usually a value less than 0.05 is appropriate for r*. For experienced users, the program also allows the direct input of the relative threshold level . A conservative choice of would be 0.001, which will ensure a satisfactory reconstruction accuracy and still produce a reasonably large compression ratio.