For computing the gradient of a scalar field, one can use the weighted least squares method as described here:

least squares gradient reconstruction

in the appendix (.pdf page 23). How to go about reconstructing a gradient of a vector field? I've read here:

least squares gradients from geophysical experimental data

  that one could consider the vector field components as scalar fields and feed this to the grad calculation, under some assumptions. I'm writing currently a generic C++ method that is supposed to compute this, where the result rank is determined by outer product trait classes defined for all combinations available (scalar vector, vector-tensor, tensor-vector, vector-vector, etc). Any thoughts on how to approach this?

For computing the gradient of a scalar field, one can use the weighted least squares method as described in the paper Revisiting the Least-squares Procedure for Gradient Reconstruction on Unstructured Meshes by Dimitri Mavriplis (pg. 23).

My question is: How can I reconstruct a gradient of a vector field?

In the paper Least-squares gradient calculation from multi-point observations of scalar and vector fields: methodology and applications with Cluster in the plasmasphere by J. De Keyser, et al., it seems that one could consider the vector field components as scalar fields and feed this to the gradient calculation, under some assumptions. 

I'm writing currently a generic C++ method that is supposed to compute this, where the result rank is determined by outer product trait classes defined for all combinations available (scalar vector, vector-tensor, tensor-vector, vector-vector, etc). Any thoughts on how to approach this?