## 3 Our Proposed Approach
### Line Search Algorithm for the L-BFGS-B  
After selecting the best CUDA optimization configuration described in Section 3.2, we integrate our CUDA implementation of the Okada model into the Monte Carlo method estimating the earthquake source parameters with the L-BFGS-B optimization. Then, we analyze the computational cost of subroutines in the L-BFGS-B to find the most time-consuming part and alleviate the bottleneck. The L-BFGS-B proceeds in the sequence of generalized Cauchy point computation (algorithm CP), subspace minimization, line search, objective function and gradient up  
Figure 6: Comparison of the L and H options for CSE. date, and update of a new limited-memory Hessian approximation, as described in Algorithm 1[30]. We define the objective function \(f\) as the misfit of geodetic measurement \(\mathbf{d}\) and model displacement \(\mathbf{G}(x)\) and use the central finite difference equation for gradient approximation, as formulated by Equations (7) and (8), respectively. ```
1:Result: optimized output \(x_{out}\)
2:Initialization: random starting point \(x_{0}\), the number of BFGS corrections stored \(m\), bound constraint \(l,\ u\), termination condition \(maxIter,\ ftol,\ gtol,k=0\)
3:repeat
4: Algorithm CP: compute generalized Cauchy point \(x_{k}^{c}\)
5: Subspace minimization: compute search direction \(d_{k}\)
6: Line search: select step length \(\alpha_{k}\)
7: Update point: \(x_{k+1}=x_{k}+\alpha_{k}d_{k}\)
8: Compute objective function \(f_{k+1}=\frac{1}{2}\|\mathbf{G}(x_{k+1}-\mathbf{d})\|_{2}^{2}\) and gradient \(g_{k+1}\)
9: Update limited memory BFGS matrix \(H_{k+1}\)
10:\(k=k+1\)
11: Compute projected gradient \(P(x_{k}-g_{k},l,u)\)
12:until\(\frac{f_{k-1}-f_{k}}{\max\left(\|f_{k}\|,\|f_{k-1}\|,1\right)}<ftol\) and \(\|P(x_{k}-g_{k},l,u)-x_{k}\|_{\infty}<gtol\) and \(k>maxIter\)
13:\(x_{out}=x_{k}\)
```  
**Algorithm 1** L-BFGS-B algorithm. \[f=\frac{1}{2}\|\mathbf{G}(x)-\mathbf{d}\|_{2}^{2} \tag{7}\]  
\[g\approx\frac{f(x+h)-f(x-h)}{h} \tag{8}\]  
As summarized in Table 4, we analyzed the computational cost of L-BFGS-B subroutines [30; 50], where \(m\) is the number of BFGS corrections stored, \(n\) is the number of variables to be optimized, \(t\) is the number of free variables, \(n_{int}\) is the number of segment explorations of the algorithm CP, \(liter\) is the number of iterations of the line search algorithm, \(w\) is the mesh width of the Okada model, and \(h\) is the mesh height of the model, respectively. We optimized nine parameters, all input parameters of the Okada model except _Open_, in this paper and usually used small values of \(m\) (e.g., \(3\leq m\leq 20\)) [50]. However, the mesh size of the Okada model varies from hundreds to thousands of squares depending on the scale of geodetic measurement, and thus the subroutine that includes the most dislocation model calculations takes the highest computational burden. Therefore, the line search algorithm can be a bottleneck of the L-BFGS-B because the Okada model is repeatedly calculated. To alleviate this bottleneck, we compared four different candidates of the line search algorithm: Armijo, MT, BWW, and cuLBFGSB. We performed the five thousand iterations of the Monte Carlo method for each candidate in the real earthquake data mentioned in \begin{table}
\begin{tabular}{l c} \hline \hline
**Subroutine** & **Computational Cost** \\ \hline Algorithm CP & \((2m+2)n+O(m^{2})\cdot n_{int}\) \\ Subspace minimization (direct primal method) & \(2m^{2}t+6 mt+4t+O(m^{3})\) \\ Line search (Armijo) & \(liter\cdot O(wh)\) \\ Line search (MT, BWW) & \(liter\cdot(2+2n)O(wh)\) \\ Update limited memory BFGS matrix & \(2n+O(m^{3})\) \\ Compute objective function \(f\) & \(O(wh)\) \\ Compute gradient \(g\) & \((1+2n)O(wh)\) \\ \hline \hline \end{tabular}
\end{table} 
Table 4: Computational cost of the L-BFGS-B subroutines. Section 3.1 and evaluated the performance and correctness of the results.