Using least squares for computing gradients

Question

I am developing a code where I am using the least squares method to compute gradients.

Generally, we use least squares to obtain some model based on a set of data (${q_1 \cdots q_N}$) at locations (${x_1 \cdots x_N}$), and use this model to predict $q_c$ at $x_c$.

In my case, I am looking for $\nabla{q}$ at $x_c$ using least squares. Should the $q$ value at $x_c$ be included in a least squares model for its gradient at the location, or only the neighboring values?

I would like use the weighted least squares method. It seems the weighting matrix is defined as a diagonal matrix with the inverse of the variance squared along the main diagonal. The popular choice for the variance seems to be the distance between the neighboring data point and the point in question. So if I include the $q_c$ value at $x_c$ in my least squares model, then the variance would be zero for this data point, which would be problematic unless I add some $\epsilon$ value to the variance before taking it's inverse.

The least squares equation is

\begin{align} \nabla q = q^TZM^T\nabla\phi(x) \\ \\ M=(Z^TZ)^{-1}\\ \\ \phi = \begin{bmatrix} 1 & x & y & z & xy & xz & yz & 0.5x^2 & 0.5y^2 & 0.5z^2 \\ \end{bmatrix}\\ \\ Z_{unweight} = \begin{bmatrix} 1 & x_1 & y_1 & z_1 & x_1y_1 & x_1z_1 & y_1z_1 & 0.5x_1^2 & 0.5y_1^2 & 0.5z_1^2 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots &\vdots & \vdots \\ 1 & x_N & y_N & z_N & x_Ny_N & x_Nz_N & y_Nz_N & 0.5x_N^2 & 0.5y_N^2 & 0.5z_N^2 \\ \end{bmatrix} \\ \\ q = \begin{bmatrix} q_1 & \cdots & q_N \end{bmatrix} \end{align}

$q_c$ is known, but I am not sure if this should be included in the above model.

Answer 1

You are asking the following question. Given $(q_0,q_1,q_2,\ldots,q_N)$ find $\nabla q$ at $\vec{r}=\vec{r}_0$. You can do this two ways.

I. Assume

$$ q(x,y,z) = a + \vec{b} \cdot (\vec{r} - \vec{r}_0) $$ and find $a$ and $\vec{b}$ by $$ \min_{a,\vec{b}}\sum_{i=0}^N(q(x_i,y_i,z_i) - q_i)^2 $$

II. Assume

$$ q(x,y,z) = q_0 + \vec{b} \cdot (\vec{r} - \vec{r}_0) $$ and find $\vec{b}$ by $$ \min_{\vec{b}}\sum_{i=1}^N(q(x_i,y_i,z_i) - q_i)^2 $$

In either case, $\vec{b}$ will give an approximation to $\nabla q$ at $\vec{r}_0$.

In (II), you have $q(x_0,y_0,z_0)=q_0$ but in (I) you only have $q(x_0,y_0,z_0) \approx q_0$.

If you want to put in weights based on some inverse distance in above norm, then even in (I) you are forcing $q(x_0,y_0,z_0)=q_0$ to hold. In this case, you can just use approach (II) after putting in the weights appropriately in the norm.