# Fast and Efficient Lunar Finite Element Gravity Model
## 2 Orthogonal Approximation  
The following are some examples of treatments of discrete approximation using Chebyshev polynomials: [24; 38; 39; 40; 41; 42; 43; 44; 45]. Although [41; 42] hold more significance in the field of orthogonal approximation, in [44], the orthogonal approximation is used to show that multi-resolution approximation of input and output maps can be linear and nonlinear. Typically, a linear least-squares approximation is used for an orthogonal approximation when dealing with experimental data or an arbitrary function generating polynomials. In this paper, a least-squares approximation with classical Chebyshev polynomials is used to analyze discrete samples of the \"to-be\" approximated function. Furthermore, arranging the regression matrix to be \"Kronecker factorable\" (a technique used in [21]) allows the use of array algebra identities to generate the multidimensional orthogonal least-squares operators. This approach can be seamlessly adapted to higher-dimensional approximations. This section summarizes how using the Kronecker product technique extends a single classical Chebyshev polynomial approximation result to higher dimensions. This method avoids dimensional challenges related to how many attributes a dataset has and paves the way for efficient computation. The importance of this method allows considering the advantages of orthogonal approximation in gravity modeling. Most mathematical derivations are based on the work presented in [21]. First, let us briefly consider the following approximation of a single-valued function, \(f(x)\), with \(x\) as one independent variable:  
\[g(x),\{x_{min}\leq x\leq x_{max}\} \tag{3}\]  
For simplification in multi-dimensional analysis, a new independent variable is introduced \(\xi\), where \(\{-1\leq\xi\leq 1\}\). Thus, through the process of normalization, the forward and inverse transformation for \(\xi\) can be produced:  
\[ \begin{split} Forward&\big{\{}\,\xi(x)=2(x-x_{ min})/(x_{max}-x_{min})-1\\ Inverse&\big{\{}\,x(\xi)=x_{min}+\,(\xi+1)\,(x_{ max}-x_{min})/2\end{split}  \tag{4}\]  
The final function to be approximated can be created by substituting the forward and inverse functions into Equation (3):  
\[f(\xi)=g(x_{min}+(\xi+1)\,(x_{max}-x_{min})/2\,) \tag{5}\]  
In a general case, for basis functions, \(\phi_{n}(\xi)\), we seek to approximate \(f(\xi)\) as a linear combination of a set of \(N+1\) linearly independent basis functions \(\{\phi_{0}(\xi),\phi_{1}(\xi),...,\phi_{N}(\xi)\}\) as the following:  
\[f(\xi)\approx\,\sum_{n=0}^{N}a_{n}\phi_{n}(\xi) \tag{6}\]  
For discrete cases, a set of sample points is introduced in which nodes is set as \(\{\xi_{0},\xi_{i},...,\xi_{M};M\geq N\}\). With this, the residual approximation error at each node is considered to be:  
\[r_{j}=f(\xi_{j})-\sum_{n=0}^{N}a_{n}\phi_{n}(\xi);\,j=0,1,...,M \tag{7}\]  
In general, the decision to use a least-squares approximation routine is motivated by the advantages of applying the one-dimensional analysis to higher dimensions. Using least-squares approximation, the objective is to find the coefficient vector \(a\) that minimizes the weighted sum square of the residuals:  
\[J=\frac{1}{2}(f-\Phi a)^{T}W(f-\Phi a);W=W^{T} \tag{8}\]  
where \(\Phi=\sum_{n=0}^{N}a_{n}\phi_{n}(\xi)\), and \(W\) is a positive definite weighted matrix determined by classical orthogonal conditions for Chebyshev polynomials. Using a weighted least-squares solution and orthogonal conditions, this result can be displayed as normal equations:  
\[a=(\Phi^{T}W\Phi)^{-1}\Phi^{T}W\mathbf{f}=C_{x}\mathbf{f} \tag{9}\]  
This finalized relationship can be displayed in matrix form as the following:\[f= \begin{bmatrix}f(\xi_{0})\\ f(\xi_{i})\\ \vdots\\ f(\xi_{N})\end{bmatrix} ,C_{x}=\frac{1}{M_{x}}