CODEWITHSUNDEEP
numpyscipypolynomials
M.F.
M.F.

Reputation: 31

Coefficients of 2D Chebyshev series in numpy.polynomial.chebyshev

I understand that chebvander2d and chebval2d return the Vandermonde matrix and fitted values for 2D inputs, and chebfit returns the coefficients for 1D-input series, but how do I get the coefficients for 2D-input series?

Upvotes: 3

Views: 1082

Answers (2)

Georges
Georges

Reputation: 39

As already mentioned, numpy.polynomial.chebyshev.chebfit2d() does not seems to exist. Fortunately, numpy.polynomial.chebyshev.chebvander2d() does exist and works as expected. The following worked for me:

import numpy as np

def chebfit2d(x0, x1, y, deg, rcond=None, full=False, w=None):
    a = np.polynomial.chebyshev.chebvander2d(x0, x1, deg)
    if w is not None:
        a = np.matmul(np.diag(w), a)
        y = np.matmul(np.diag(w), y)
    x, residuals, rank, s = np.linalg.lstsq(a, y, rcond=rcond)
    if full:
        return(x.reshape(deg[0]+1, deg[1]+1), [residuals, rank, s, rcond])
    return(x.reshape(deg[0]+1, deg[1]+1))

The output can directly be passed to numpy.polynomial.chebyshev.chebval2d()

As numpy.polynomial.chebyshev.chebvander3d() also exists, the 3D extension should be straightforward. As a reminder, for an optimal use, x0 and x1 should be normalized so that they fill the [-1, +1] range.

Regarding Chris' comments, yes, the total degree of the polynomial is deg[0]+deg[1] and not all the possible terms for a polynomial of this degree are present (or independent). I see no issue there.

Upvotes: 1

Chris K
Chris K

Reputation: 1576

Short answer: It looks to me like this is not yet implemented. The whole of 2D polynomials seems more like a draft with some stub functions (as of June 2020).

Long answer (I came looking for the same thing, so I dug a little deeper):

First of all, this applies to all of the polynomial classes, not only chebyshev, so you also cannot fit an "ordinary" polynomial (power series). In fact, you cannot even construct one.

To understand the programming problem, let me recapture what a 2D polynomial looks like as a math formula, at an example polynomial of degree 2:
p(x, y) = c_00 + c_10 x + c_01 y + c_20 x^2 + c11 xy + c02 y^2
here the indices of c refer to the powers of x and y (the sum of the exponents must be <= degree).

First thing to notice is that, for degree d, there are (d+1)(d+2)/2 coefficients. They could be stored in the upper left part of a matrix or in a 1D array, e.g. aranged as in the formula above.
The documentation of functions like numpy.polynomial.polynomial.polyval2d implies that numpy expects the matrix variant: p(x, y) = sum_i,j c_i,j * x^i * y^j.
Side note: it may be confusing that the row index i ("y-coordinate") of the matrix is used as exponent of x, not y; maybe the role of i and j should be switched if this is eventually implementd, or at least there should be a note in the documentation.

This leads to the core problem: the data structure for the 2D coefficients is not defined anywhere; only indirectly, like above, it can be guessed that a matrix should be used. But compared to a 1D array this is a waste of space, and evaluation of the polynomial takes two nested loops instead of just one. Also: does the matrix have to be initialized with np.zeros or do the implemented functions make sure that the lower right part is never touched so that np.empty can be used?

If the whole (d+1)^2 matrix were used, as the polyval2d function doc suggests, the degree of the polynomial would actually be d*2 (if c_d,d != 0)
To test this, I wanted to construct a numpy.polynomial.polynomial.Polynomial (yes, three times polynomial) and check the degree attribute:

import numpy as np
import numpy.polynomial.polynomial as poly

coef = np.array([
  [5.00, 5.01, 5.02],
  [5.10, 5.11, 0.  ],
  [5.20, 0.  , 0.  ]
])

polyObj = poly.Polynomial(coef)
print(polyObj.degree)

This gave a ValueError: Coefficient array is not 1-d before the print statement was reached. So while polyval2d expects a 2D coefficient array, it is not (yet) possible to construct such a polynomial - not manually like this at least. With this insight, it is not surprising that there is no function (yet) that computes a fit for 2D polynomials.

Upvotes: 1

Related Questions