## 4 Numerical Procedures  
Based on the discussions in the previous sections, we illustrate the following numerical workflows for the inversion of fluid content from NMR measurements of whole cores.  
Finding the optimal solution of the \(L_{p}\) norm (\(p<2\))-based optimization (9) can be achieved by applying the weight in Equation (16), followed by finding the solution to Equation (15). The procedure is iterative reweighted least squares [26]. The keyword \"iterative\" reflects the nature of the algorithm that the current weight (16) is dependent on the previous solution. There are a number of numerical algorithms that have been designed to solve Equation (15). Among them, least squares with QR-factorization (LSQR) [27] and conjugate gradient for least squares (CGLS) problems [28] have both been proven to be efficient and convergent. Their implementations are also straightforward. We will use the CGLS algorithm in this paper.  
Three numerical procedures need to be implemented: CGLS for solving Equation (13), IRLS for solving the weighted system (15) of the \(L_{p}\)-norm optimization (9), and inversion conditioning. The CGLS procedure is used in the IRLS during each iteration. In each of the following workflows, we assume \(\mathbf{A}\) is a \(M\times N\) matrix, \(\mathbf{m}\) is a column vector with \(N\) elements, and \(\mathbf{d}\) is a column vector with \(M\) elements.  
The CGLS procedure for solving the linear equation \(\mathbf{Am}=\mathbf{d}\) is given as follows:  
* Choose \(\mathbf{m}_{0}\). Set \(\mathbf{s}_{0}=\mathbf{d}-\mathbf{Am}_{0}\), \(\mathbf{r}_{0}=\mathbf{p}_{0}=\mathbf{A}^{T}(\mathbf{d}-\mathbf{Am}_{0})= \mathbf{A}^{T}\mathbf{s}_{0}\), \(\mathbf{q}_{0}=\mathbf{Ap}_{0}\);
* Repeat the following steps until the solution \(\mathbf{m}_{k}\) does not change or \(k\) reaches an upper limit:
* \(\alpha_{k+1}=\frac{(\mathbf{r}_{k},\mathbf{r}_{k})}{(\mathbf{q}_{k},\mathbf{q }_{k})}\);
* \(\beta_{k+1}=\frac{(\mathbf{r}_{k+1},\mathbf{r}_{k+1})}{(\mathbf{r}_{k},\mathbf{ r}_{k})}\);  
Figure 4: Velocity inversion using travel time information under high-frequency assumption. The blue line indicates the source-receiver pair where the travel time value has the wrong measurement, the source at \((3,8)\) m, and the receiver is at \((9,0)\) m. (**a**) The model grids and straight rays between source-receiver pairs. (**b**) The inversion results using perfect data with \(p=2.0\). (**c–f**), Inversion results using the wrong data with \(p=1.3,1.2,1.1,\) and \(1.0\), respectively.  
The IRLS workflow for the \(L_{p}\) norm optimization problem for linear equation \(\mathbf{Am}=\mathbf{d}\) is illustrated below:  
* Choose \(\mathbf{m}_{0}\) and a small positive number \(\varepsilon\);
* Repeat the following steps until the solution converges:
* \(r_{i}=d_{i}-\sum\limits_{j=1}^{N}A_{ij}m_{j}\), \(i=1,\cdots,M\);
* Compute the weights for, \(W_{i}=\left\{\begin{array}{ll}|r_{i}|^{p-2}&,\quad|r_{i}|>\varepsilon\\ \varepsilon^{p-2}&,\quad|r_{i}|\leq\varepsilon\end{array},i=1,\cdots,M\);
* Apply the weights to the data, \(z_{i}=W_{i}d_{i}\), \(i=1,\cdots,M\);
* Apply the weights to the linear matrix, \(\mathbf{B}=\sqrt{W}\mathbf{A}\);
* Solve the new equation \(\mathbf{B}\mathbf{m}=\mathbf{z}\) using CGLS.  
There are two loops in the IRLS algorithm. The inner loop is inside the CGLS for solving linear equations as an \(L_{2}\) norm problem. The outer loop is for the computation of the new matrix \(\mathbf{B}\) and the new data vector \(\mathbf{z}\) from \(\mathbf{A}\) and \(\mathbf{d}\) based on the _current_ model estimation. It is this automatic updating for the equations that make IRLS immune to the impact of the data measurements with large errors.  
Lastly, the inversion conditioning is implemented as a non-stationary smoothing operator [15; 29] and applied after the IRLS workflow. The smoothing operation can also be implemented by adding the regularization term to the original cost function (9).