This is a lengthy answer. It will, in order: (1) explain what vectorization is (2) what the documentation wants you to do (3) give a modified function to match that documentation (4) how to confirm we generalized the function without changing it.
(1) What vectorization means in the context of such python-functions is easiest to understand if we consider a more simple function, taking a real number and returning a real number:
# ℝ → ℝ
def f_normal(x: float) -> float:
# take a number, double it
return x*2
# ℝⁿ → ℝⁿ
def f_vec(x: ndarray) -> ndarray:
# take an array of numbers, return an equal-sized array where every number is doubled
return x*2
These functions calculate the same thing conceptually, but one works on single values, the other works on a lot of values at once. Why would you do that? Loops in python are slow. Operating on many values at once at native-coded-library-level is a huge speed gain. In this simple case, the functions are literally the same. This is no coincidence. Numpy is designed around pythons duck-typing in such a way that this works. But we can take this even further: The function also works for ℝⁿˣᵐ → ℝⁿˣᵐ, for ℝⁿˣᵐˣᵏ → ℝⁿˣᵐˣᵏ. We say it works for all shapes. scipy.differentiate.hessian relies on this property of the function! Plug this toy function into the hessian and print the argument x. It is never just a number, it is multi-dimensional matrices.
Now consider your function which takes an array and returns a number:
# ℝᵐ -> ℝ
def func(x: ndarray) -> float:
# takes an array, returns a number
return np.sum(np.cross(x[0:3], x[3:6]))
(2) Now we can understand the critical part of the doc:
argument f must be vectorized to accept an array of shape (m, ...). The first axis represents the m inputs of f; the remaining axes indicated by ellipses are for evaluating the function at several abscissae in a single call.
argument f must return an array of shape (...).
We must modify your function so that it can not just doℝᵐ → ℝ but also ℝᵐˣⁿ → ℝⁿ and ℝᵐˣⁿˣᵏ → ℝⁿˣᵏ (and more).
So now that we understand the goal, how can we modify the function to work the way scipy expects it to?
The trouble is this expression: x[0:3], will not be returning a 3D-vector, it will return any multidimensional matrix, we can only describe the shape as (3, ...), and np.cross (and np.sum) will not know what to do with such a thing. We can tell it to operate on this first dimension and essentially loop over the other numbers with the axis argument of np.cross.
(3)
def vec_func(x: np.ndarray) -> np.ndarray:
return np.sum(np.cross(x[0:3], x[3:6], axis=0), axis=0)
(4)
This works with the hessian-function from scipy. We can informally confirm this function does produce the same results as the non-vectorized code with this snippet:
inp = np.array([[1, 2, 3, 4, 5, 6], [2, 3, 4, 5, 9, 0], [8, 5, 8, 1, -9, 10]])
print(func(inp[0]), func(inp[1]), func(inp[2])) # 0 -13 -27
print(vec_func(inp[0]), vec_func(inp[1]), vec_func(inp[2])) # 0 -13 -27
# note the T for transpose
print(vec_func(inp.T)) # [ 0 -13 -27 ]
