# Parameter uncertainty in Least Squares fit to Vector Valued Function

Suppose I have to fit a vector field $\vec{y}_i$ with a vector valued function $\vec{f}_i = \vec{f}(\vec{x}_i; \vec{p})$, so both $\vec{y}_i, \vec{f}_i \in \mathbb{R}^m$, where $\vec{x}_i \in \mathbb{R}^n$, and the index $i$ indicates the i-th datapoint. We are interested in the values of the parameters $\vec{p} = (p_1, \ldots, p_k)$; these can be obtained by minimizing the least squares

$\chi^2=\sum_{i=1}^N \Big\vert \frac{\vec{f}(\vec{x}_i) - \vec{y}_i}{\vec{\sigma}_i} \Big\vert^2 = \sum_{i=1}^N \sum_{j=1}^m \Big( \frac{f_j(\vec{x}_i) - y_{ji}}{\sigma_{ji}} \Big)^2$

My question: what is the correct way to estimate the uncertainties in the parameters $\vec{p}$ of this model?

What I have tried: For a scalar function $g_i = g(\vec{x}_i;\vec{p})$, we typically approximate the covariance matrix as

$C = (J_{\vec{p}}^{T} W J_{\vec{p}})^{-1}$,

where $J_{\vec{p}}=\nabla_{\vec{p}} g(\vec{x}_i;\vec{p})$ is the jacobian of $g_i$ with respect to the parameters in the model, and $W$ is a diagonal matrix with the values of $1/\sigma_i^2$ on its diagonal. This way, $J_{\vec{p}}$ is a $k \times N$ matrix, and $W$ is $N \times N$ ($N$ is the number of datapoints).

I therefore assumed the natural way to generalize this to the vector model $\vec{f}_i$, is to make $J$ and $W$ higher dimensional tensors, such that $J_{\vec{p}}$ is $k \times m \times N$, and W is $m\times N \times N \times m$. This way, if I dot product them together as $J_{\vec{p}}^{T} W J_{\vec{p}}$ appropriately, the result will still be a $k \times k$ matrix as it should be, which after inversion should be the covariance matrix. But the results do not quite match my expectation (though they are close) so I want to make sure this is the right way to deal with this. Is my generalization correct, or if not, what should it be?