The np.polynomial.Polynomial does not seem to handle array-valued polynomials, although this is advertised at some places in its documentation.

E.g. help(np.polynomial.Polynomial.fit) says:

"yarray_like, shape (M,) or (M, K): y-coordinates of the sample points. Several data sets of sample points sharing the same x-coordinates can be fitted at once by passing in a 2D-array that contains one dataset per column."

However, it doesn't seem to work. Maybe this is a left-over "intended feature" that didn't make it in the final implementation, or the documentation is just wrong. Even if so, this feature would be useful to have.

Reproducing code example:

import numpy as np
np.polynomial.Polynomial.fit([0, 1, 2], np.array([[1, 2, 3], [4, 5, 6]]).T, deg=2)

Error message:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib64/python3.8/site-packages/numpy/polynomial/_polybase.py", line 898, in fit
    return cls(coef, domain=domain, window=window)
  File "/usr/lib64/python3.8/site-packages/numpy/polynomial/_polybase.py", line 266, in __init__
    [coef] = pu.as_series([coef], trim=False)
  File "/usr/lib64/python3.8/site-packages/numpy/polynomial/polyutils.py", line 182, in as_series
    raise ValueError("Coefficient array is not 1-d")
ValueError: Coefficient array is not 1-d

Numpy/Python version information:

1.18.4 3.8.2 (default, Feb 28 2020, 00:00:00)