ENH: Add functions for producing 2D and 3D pseudo Vandermonde matrices that are useful for least squares fits to data depending on two or three variables using the various polynomial basis.

The new functions have names polyvander2d, and polyvander3d,
where 'poly' can be replaced by any of 'leg', 'cheb', 'lag',
'herm', or 'herme'.

6 files changed, 589 insertions, 18 deletions

diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py
index 6212f2bc5..d6ccf25ca 100644
--- a/numpy/polynomial/chebyshev.py
+++ b/numpy/polynomial/chebyshev.py
@@ -40,6 +40,8 @@ Misc Functions
 - `chebfromroots` -- create a Chebyshev series with specified roots.
 - `chebroots` -- find the roots of a Chebyshev series.
 - `chebvander` -- Vandermonde-like matrix for Chebyshev polynomials.
+- `chebvander2d` -- Vandermonde-like matrix for 2D power series.
+- `chebvander3d` -- Vandermonde-like matrix for 3D power series.
 - `chebfit` -- least-squares fit returning a Chebyshev series.
 - `chebpts1` -- Chebyshev points of the first kind.
 - `chebpts2` -- Chebyshev points of the second kind.
@@ -90,10 +92,10 @@ from polytemplate import polytemplate
 
 __all__ = ['chebzero', 'chebone', 'chebx', 'chebdomain', 'chebline',
     'chebadd', 'chebsub', 'chebmulx', 'chebmul', 'chebdiv', 'chebpow',
-    'chebval', 'chebval2d', 'chebval3d', 'chebgrid2d', 'chebgrid3d',
-    'chebder', 'chebint', 'cheb2poly', 'poly2cheb', 'chebfromroots',
-    'chebvander', 'chebfit', 'chebtrim', 'chebroots', 'chebpts1',
-    'chebpts2', 'Chebyshev']
+    'chebval', 'chebder', 'chebint', 'cheb2poly', 'poly2cheb',
+    'chebfromroots', 'chebvander', 'chebfit', 'chebtrim', 'chebroots',
+    'chebpts1', 'chebpts2', 'Chebyshev', 'chebval2d', 'chebval3d',
+    'chebgrid2d', 'chebgrid3d', 'chebvander2d','chebvander3d']
 
 chebtrim = pu.trimcoef
 
@@ -1315,6 +1317,98 @@ def chebvander(x, deg) :
     return np.rollaxis(v, 0, v.ndim)
 
 
+def chebvander2d(x, y, deg) :
+    """Psuedo Vandermonde matrix of given degree.
+
+    Returns the pseudo Vandermonde matrix for 2D Chebyshev series in `x`
+    and `y`. The sample point coordinates must all have the same shape
+    after conversion to arrays and the dtype will be converted to either
+    float64 or complex128 depending on whether any of `x` or 'y' are
+    complex.  The maximum degrees of the 2D Chebyshev series in each
+    variable are specified in the list `deg` in the form ``[xdeg, ydeg]``.
+    The return array has the shape ``x.shape + (order,)`` if `x`, and `y`
+    are arrays or ``(1, order) if they are scalars. Here order is the
+    number of elements in a flattened coefficient array of original shape
+    ``(xdeg + 1, ydeg + 1)``.  The flattening is done so that the resulting
+    pseudo Vandermonde array can be easily used in least squares fits.
+
+    Parameters
+    ----------
+    x,y : array_like
+        Arrays of point coordinates, each of the same shape.
+    deg : list
+        List of maximum degrees of the form [x_deg, y_deg].
+
+    Returns
+    -------
+    vander2d : ndarray
+        The shape of the returned matrix is described above.
+
+    See Also
+    --------
+    chebvander, chebvander3d. chebval2d, chebval3d
+
+    """
+    ideg = [int(d) for d in deg]
+    is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
+    if is_valid != [1, 1]:
+        raise ValueError("degrees must be non-negative integers")
+    degx, degy = deg
+    x, y = np.array((x, y), copy=0) + 0.0
+
+    vx = chebvander(x, degx)
+    vy = chebvander(y, degy)
+    v = np.einsum("...i,...j->...ij", vx, vy)
+    return v.reshape(v.shape[:-2] + (-1,))
+
+
+def chebvander3d(x, y, z, deg) :
+    """Psuedo Vandermonde matrix of given degree.
+
+    Returns the pseudo Vandermonde matrix for 3D Chebyshev series in `x`,
+    `y`, or `z`. The sample point coordinates must all have the same shape
+    after conversion to arrays and the dtype will be converted to either
+    float64 or complex128 depending on whether any of `x`, `y`, or 'z' are
+    complex.  The maximum degrees of the 3D Chebeshev series in each
+    variable are specified in the list `deg` in the form ``[xdeg, ydeg,
+    zdeg]``. The return array has the shape ``x.shape + (order,)`` if `x`,
+    `y`, and `z` are arrays or ``(1, order) if they are scalars. Here order
+    is the number of elements in a flattened coefficient array of original
+    shape ``(xdeg + 1, ydeg + 1, zdeg + 1)``.  The flattening is done so
+    that the resulting pseudo Vandermonde array can be easily used in least
+    squares fits.
+
+    Parameters
+    ----------
+    x,y,z : array_like
+        Arrays of point coordinates, each of the same shape.
+    deg : list
+        List of maximum degrees of the form [x_deg, y_deg, z_deg].
+
+    Returns
+    -------
+    vander3d : ndarray
+        The shape of the returned matrix is described above.
+
+    See Also
+    --------
+    chebvander, chebvander3d. chebval2d, chebval3d
+
+    """
+    ideg = [int(d) for d in deg]
+    is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
+    if is_valid != [1, 1, 1]:
+        raise ValueError("degrees must be non-negative integers")
+    degx, degy, degz = deg
+    x, y, z = np.array((x, y, z), copy=0) + 0.0
+
+    vx = chebvander(x, deg_x)
+    vy = chebvander(y, deg_y)
+    vz = chebvander(z, deg_z)
+    v = np.einsum("...i,...j,...k->...ijk", vx, vy, vz)
+    return v.reshape(v.shape[:-3] + (-1,))
+
+
 def chebfit(x, y, deg, rcond=None, full=False, w=None):
     """
     Least squares fit of Chebyshev series to data.
diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py
index 2fd28a3ff..801e8f8d9 100644
--- a/numpy/polynomial/hermite.py
+++ b/numpy/polynomial/hermite.py
@@ -37,6 +37,8 @@ Misc Functions
 - `hermfromroots` -- create a Hermite series with specified roots.
 - `hermroots` -- find the roots of a Hermite series.
 - `hermvander` -- Vandermonde-like matrix for Hermite polynomials.
+- `hermvander2d` -- Vandermonde-like matrix for 2D power series.
+- `hermvander3d` -- Vandermonde-like matrix for 3D power series.
 - `hermfit` -- least-squares fit returning a Hermite series.
 - `hermtrim` -- trim leading coefficients from a Hermite series.
 - `hermline` -- Hermite series of given straight line.
@@ -62,9 +64,10 @@ from polytemplate import polytemplate
 
 __all__ = ['hermzero', 'hermone', 'hermx', 'hermdomain', 'hermline',
     'hermadd', 'hermsub', 'hermmulx', 'hermmul', 'hermdiv', 'hermpow',
-    'hermval', 'hermval2d', 'hermval3d', 'hermgrid2d', 'hermgrid3d',
-    'hermder', 'hermint', 'herm2poly', 'poly2herm', 'hermfromroots',
-    'hermvander', 'hermfit', 'hermtrim', 'hermroots', 'Hermite']
+    'hermval', 'hermder', 'hermint', 'herm2poly', 'poly2herm',
+    'hermfromroots', 'hermvander', 'hermfit', 'hermtrim', 'hermroots',
+    'Hermite', 'hermval2d', 'hermval3d', 'hermgrid2d', 'hermgrid3d',
+    'hermvander2d', 'hermvander3d']
 
 hermtrim = pu.trimcoef
 
@@ -1102,6 +1105,98 @@ def hermvander(x, deg) :
     return np.rollaxis(v, 0, v.ndim)
 
 
+def hermvander2d(x, y, deg) :
+    """Pseudo Vandermonde matrix of given degree.
+
+    Returns the pseudo Vandermonde matrix for 2D Hermite series in `x` and
+    `y`. The sample point coordinates must all have the same shape after
+    conversion to arrays and the dtype will be converted to either float64
+    or complex128 depending on whether any of `x` or 'y' are complex.  The
+    maximum degrees of the 2D Hermite series in each variable are specified
+    in the list `deg` in the form ``[xdeg, ydeg]``. The return array has
+    the shape ``x.shape + (order,)`` if `x`, and `y` are arrays or
+    ``(1, order) if they are scalars. Here order is the number of elements
+    in a flattened coefficient array of original shape ``(xdeg + 1, ydeg +
+    1)``.  The flattening is done so that the resulting pseudo Vandermonde
+    array can be easily used in least squares fits.
+
+    Parameters
+    ----------
+    x,y : array_like
+        Arrays of point coordinates, each of the same shape.
+    deg : list
+        List of maximum degrees of the form [x_deg, y_deg].
+
+    Returns
+    -------
+    vander2d : ndarray
+        The shape of the returned matrix is described above.
+
+    See Also
+    --------
+    hermvander, hermvander3d. hermval2d, hermval3d
+
+    """
+    ideg = [int(d) for d in deg]
+    is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
+    if is_valid != [1, 1]:
+        raise ValueError("degrees must be non-negative integers")
+    degx, degy = deg
+    x, y = np.array((x, y), copy=0) + 0.0
+
+    vx = hermvander(x, degx)
+    vy = hermvander(y, degy)
+    v = np.einsum("...i,...j->...ij", vx, vy)
+    return v.reshape(v.shape[:-2] + (-1,))
+
+
+def hermvander3d(x, y, z, deg) :
+    """Psuedo Vandermonde matrix of given degree.
+
+    Returns the pseudo Vandermonde matrix for 3D Hermite series in `x`,
+    `y`, or `z`. The sample point coordinates must all have the same shape
+    after conversion to arrays and the dtype will be converted to either
+    float64 or complex128 depending on whether any of `x`, `y`, or 'z' are
+    complex.  The maximum degrees of the 3D Hermite series in each variable
+    are specified in the list `deg` in the form ``[xdeg, ydeg, zdeg]``. The
+    return array has the shape ``x.shape + (order,)`` if `x`, `y`, and `z`
+    are arrays or ``(1, order) if they are scalars. Here order is the
+    number of elements in a flattened coefficient array of original shape
+    ``(xdeg + 1, ydeg + 1, zdeg + 1)``.  The flattening is done so that the
+    resulting pseudo Vandermonde array can be easily used in least squares
+    fits.
+
+    Parameters
+    ----------
+    x,y,z : array_like
+        Arrays of point coordinates, each of the same shape.
+    deg : list
+        List of maximum degrees of the form [x_deg, y_deg, z_deg].
+
+    Returns
+    -------
+    vander3d : ndarray
+        The shape of the returned matrix is described above.
+
+    See Also
+    --------
+    hermvander, hermvander3d. hermval2d, hermval3d
+
+    """
+    ideg = [int(d) for d in deg]
+    is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
+    if is_valid != [1, 1, 1]:
+        raise ValueError("degrees must be non-negative integers")
+    degx, degy, degz = deg
+    x, y, z = np.array((x, y, z), copy=0) + 0.0
+
+    vx = hermvander(x, deg_x)
+    vy = hermvander(y, deg_y)
+    vz = hermvander(z, deg_z)
+    v = np.einsum("...i,...j,...k->...ijk", vx, vy, vz)
+    return v.reshape(v.shape[:-3] + (-1,))
+
+
 def hermfit(x, y, deg, rcond=None, full=False, w=None):
     """
     Least squares fit of Hermite series to data.
diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py
index 17e285e15..8307cf80c 100644
--- a/numpy/polynomial/hermite_e.py
+++ b/numpy/polynomial/hermite_e.py
@@ -37,6 +37,8 @@ Misc Functions
 - `hermefromroots` -- create a Hermite_e series with specified roots.
 - `hermeroots` -- find the roots of a Hermite_e series.
 - `hermevander` -- Vandermonde-like matrix for Hermite_e polynomials.
+- `hermevander2d` -- Vandermonde-like matrix for 2D power series.
+- `hermevander3d` -- Vandermonde-like matrix for 3D power series.
 - `hermefit` -- least-squares fit returning a Hermite_e series.
 - `hermetrim` -- trim leading coefficients from a Hermite_e series.
 - `hermeline` -- Hermite_e series of given straight line.
@@ -62,9 +64,11 @@ from polytemplate import polytemplate
 
 __all__ = ['hermezero', 'hermeone', 'hermex', 'hermedomain', 'hermeline',
     'hermeadd', 'hermesub', 'hermemulx', 'hermemul', 'hermediv', 'hermpow',
-    'legval', 'legval2d', 'legval3d', 'leggrid2d', 'leggrid3d',
+    'hermeval',
     'hermeder', 'hermeint', 'herme2poly', 'poly2herme', 'hermefromroots',
-    'hermevander', 'hermefit', 'hermetrim', 'hermeroots', 'HermiteE']
+    'hermevander', 'hermefit', 'hermetrim', 'hermeroots', 'HermiteE',
+    'hermeval2d', 'hermeval3d', 'hermegrid2d', 'hermegrid3d', 'hermevander2d',
+    'hermevander3d']
 
 hermetrim = pu.trimcoef
 
@@ -1097,6 +1101,98 @@ def hermevander(x, deg) :
     return np.rollaxis(v, 0, v.ndim)
 
 
+def hermevander2d(x, y, deg) :
+    """Pseudo Vandermonde matrix of given degree.
+
+    Returns the pseudo Vandermonde matrix for 2D Hermite_e series in `x`
+    and `y`. The sample point coordinates must all have the same shape
+    after conversion to arrays and the dtype will be converted to either
+    float64 or complex128 depending on whether any of `x` or 'y' are
+    complex.  The maximum degrees of the 2D Hermite_e series in each
+    variable are specified in the list `deg` in the form ``[xdeg, ydeg]``.
+    The return array has the shape ``x.shape + (order,)`` if `x`, and `y`
+    are arrays or ``(1, order) if they are scalars. Here order is the
+    number of elements in a flattened coefficient array of original shape
+    ``(xdeg + 1, ydeg + 1)``.  The flattening is done so that the resulting
+    pseudo Vandermonde array can be easily used in least squares fits.
+
+    Parameters
+    ----------
+    x,y : array_like
+        Arrays of point coordinates, each of the same shape.
+    deg : list
+        List of maximum degrees of the form [x_deg, y_deg].
+
+    Returns
+    -------
+    vander2d : ndarray
+        The shape of the returned matrix is described above.
+
+    See Also
+    --------
+    hermevander, hermevander3d. hermeval2d, hermeval3d
+
+    """
+    ideg = [int(d) for d in deg]
+    is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
+    if is_valid != [1, 1]:
+        raise ValueError("degrees must be non-negative integers")
+    degx, degy = deg
+    x, y = np.array((x, y), copy=0) + 0.0
+
+    vx = hermevander(x, degx)
+    vy = hermevander(y, degy)
+    v = np.einsum("...i,...j->...ij", vx, vy)
+    return v.reshape(v.shape[:-2] + (-1,))
+
+
+def hermevander3d(x, y, z, deg) :
+    """Psuedo Vandermonde matrix of given degree.
+
+    Returns the pseudo Vandermonde matrix for 3D Hermite_e series in `x`,
+    `y`, or `z`. The sample point coordinates must all have the same shape
+    after conversion to arrays and the dtype will be converted to either
+    float64 or complex128 depending on whether any of `x`, `y`, or 'z' are
+    complex.  The maximum degrees of the 3D Hermite_e series in each
+    variable are specified in the list `deg` in the form ``[xdeg, ydeg,
+    zdeg]``. The return array has the shape ``x.shape + (order,)`` if `x`,
+    `y`, and `z` are arrays or ``(1, order) if they are scalars. Here order
+    is the number of elements in a flattened coefficient array of original
+    shape ``(xdeg + 1, ydeg + 1, zdeg + 1)``.  The flattening is done so
+    that the resulting pseudo Vandermonde array can be easily used in least
+    squares fits.
+
+    Parameters
+    ----------
+    x,y,z : array_like
+        Arrays of point coordinates, each of the same shape.
+    deg : list
+        List of maximum degrees of the form [x_deg, y_deg, z_deg].
+
+    Returns
+    -------
+    vander3d : ndarray
+        The shape of the returned matrix is described above.
+
+    See Also
+    --------
+    hermevander, hermevander3d. hermeval2d, hermeval3d
+
+    """
+    ideg = [int(d) for d in deg]
+    is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
+    if is_valid != [1, 1, 1]:
+        raise ValueError("degrees must be non-negative integers")
+    degx, degy, degz = deg
+    x, y, z = np.array((x, y, z), copy=0) + 0.0
+
+    vx = hermevander(x, deg_x)
+    vy = hermevander(y, deg_y)
+    vz = hermevander(z, deg_z)
+    v = np.einsum("...i,...j,...k->...ijk", vx, vy, vz)
+    return v.reshape(v.shape[:-3] + (-1,))
+
+
 def hermefit(x, y, deg, rcond=None, full=False, w=None):
     """
     Least squares fit of Hermite series to data.
diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py
index dbbb2c401..39463b0cb 100644
--- a/numpy/polynomial/laguerre.py
+++ b/numpy/polynomial/laguerre.py
@@ -37,6 +37,8 @@ Misc Functions
 - `lagfromroots` -- create a Laguerre series with specified roots.
 - `lagroots` -- find the roots of a Laguerre series.
 - `lagvander` -- Vandermonde-like matrix for Laguerre polynomials.
+- `lagvander2d` -- Vandermonde-like matrix for 2D power series.
+- `lagvander3d` -- Vandermonde-like matrix for 3D power series.
 - `lagfit` -- least-squares fit returning a Laguerre series.
 - `lagtrim` -- trim leading coefficients from a Laguerre series.
 - `lagline` -- Laguerre series of given straight line.
@@ -62,9 +64,9 @@ from polytemplate import polytemplate
 
 __all__ = ['lagzero', 'lagone', 'lagx', 'lagdomain', 'lagline',
     'lagadd', 'lagsub', 'lagmulx', 'lagmul', 'lagdiv', 'lagpow',
-    'lagval', 'lagval2d', 'lagval3d', 'laggrid2d', 'laggrid3d',
-    'lagder', 'lagint', 'lag2poly', 'poly2lag', 'lagfromroots',
-    'lagvander', 'lagfit', 'lagtrim', 'lagroots', 'Laguerre']
+    'lagval', 'lagder', 'lagint', 'lag2poly', 'poly2lag', 'lagfromroots',
+    'lagvander', 'lagfit', 'lagtrim', 'lagroots', 'Laguerre', 'lagval2d',
+    'lagval3d', 'laggrid2d', 'laggrid3d', 'lagvander2d', 'lagvander3d']
 
 lagtrim = pu.trimcoef
 
@@ -1100,6 +1102,98 @@ def lagvander(x, deg) :
     return np.rollaxis(v, 0, v.ndim)
 
 
+def lagvander2d(x, y, deg) :
+    """Pseudo Vandermonde matrix of given degree.
+
+    Returns the pseudo Vandermonde matrix for 2D Laguerre series in `x` and
+    `y`. The sample point coordinates must all have the same shape after
+    conversion to arrays and the dtype will be converted to either float64
+    or complex128 depending on whether any of `x` or 'y' are complex.  The
+    maximum degrees of the 2D Laguerre series in each variable are specified in
+    the list `deg` in the form ``[xdeg, ydeg]``. The return array has the
+    shape ``x.shape + (order,)`` if `x`, and `y` are arrays or ``(1, order)
+    if they are scalars. Here order is the number of elements in a
+    flattened coefficient array of original shape ``(xdeg + 1, ydeg + 1)``.
+    The flattening is done so that the resulting pseudo Vandermonde array
+    can be easily used in least squares fits.
+
+    Parameters
+    ----------
+    x,y : array_like
+        Arrays of point coordinates, each of the same shape.
+    deg : list
+        List of maximum degrees of the form [x_deg, y_deg].
+
+    Returns
+    -------
+    vander2d : ndarray
+        The shape of the returned matrix is described above.
+
+    See Also
+    --------
+    lagvander, lagvander3d. lagval2d, lagval3d
+
+    """
+    ideg = [int(d) for d in deg]
+    is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
+    if is_valid != [1, 1]:
+        raise ValueError("degrees must be non-negative integers")
+    degx, degy = deg
+    x, y = np.array((x, y), copy=0) + 0.0
+
+    vx = lagvander(x, degx)
+    vy = lagvander(y, degy)
+    v = np.einsum("...i,...j->...ij", vx, vy)
+    return v.reshape(v.shape[:-2] + (-1,))
+
+
+def lagvander3d(x, y, z, deg) :
+    """Psuedo Vandermonde matrix of given degree.
+
+    Returns the pseudo Vandermonde matrix for 3D Laguerre series in `x`,
+    `y`, or `z`. The sample point coordinates must all have the same shape
+    after conversion to arrays and the dtype will be converted to either
+    float64 or complex128 depending on whether any of `x`, `y`, or 'z' are
+    complex.  The maximum degrees of the 3D Laguerre series in each
+    variable are specified in the list `deg` in the form ``[xdeg, ydeg,
+    zdeg]``. The return array has the shape ``x.shape + (order,)`` if `x`,
+    `y`, and `z` are arrays or ``(1, order) if they are scalars. Here order
+    is the number of elements in a flattened coefficient array of original
+    shape ``(xdeg + 1, ydeg + 1, zdeg + 1)``.  The flattening is done so
+    that the resulting pseudo Vandermonde array can be easily used in least
+    squares fits.
+
+    Parameters
+    ----------
+    x,y,z : array_like
+        Arrays of point coordinates, each of the same shape.
+    deg : list
+        List of maximum degrees of the form [x_deg, y_deg, z_deg].
+
+    Returns
+    -------
+    vander3d : ndarray
+        The shape of the returned matrix is described above.
+
+    See Also
+    --------
+    lagvander, lagvander3d. lagval2d, lagval3d
+
+    """
+    ideg = [int(d) for d in deg]
+    is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
+    if is_valid != [1, 1, 1]:
+        raise ValueError("degrees must be non-negative integers")
+    degx, degy, degz = deg
+    x, y, z = np.array((x, y, z), copy=0) + 0.0
+
+    vx = lagvander(x, deg_x)
+    vy = lagvander(y, deg_y)
+    vz = lagvander(z, deg_z)
+    v = np.einsum("...i,...j,...k->...ijk", vx, vy, vz)
+    return v.reshape(v.shape[:-3] + (-1,))
+
+
 def lagfit(x, y, deg, rcond=None, full=False, w=None):
     """
     Least squares fit of Laguerre series to data.
diff --git a/numpy/polynomial/legendre.py b/numpy/polynomial/legendre.py
index 5a72217d0..00dbacebe 100644
--- a/numpy/polynomial/legendre.py
+++ b/numpy/polynomial/legendre.py
@@ -38,6 +38,8 @@ Misc Functions
 - `legfromroots` -- create a Legendre series with specified roots.
 - `legroots` -- find the roots of a Legendre series.
 - `legvander` -- Vandermonde-like matrix for Legendre polynomials.
+- `legvander2d` -- Vandermonde-like matrix for 2D power series.
+- `legvander3d` -- Vandermonde-like matrix for 3D power series.
 - `legfit` -- least-squares fit returning a Legendre series.
 - `legtrim` -- trim leading coefficients from a Legendre series.
 - `legline` -- Legendre series representing given straight line.
@@ -62,10 +64,10 @@ import warnings
 from polytemplate import polytemplate
 
 __all__ = ['legzero', 'legone', 'legx', 'legdomain', 'legline',
-    'legadd', 'legsub', 'legmulx', 'legmul', 'legdiv', 'legpow',
-    'legval', 'legval2d', 'legval3d', 'leggrid2d', 'leggrid3d',
+    'legadd', 'legsub', 'legmulx', 'legmul', 'legdiv', 'legpow', 'legval',
     'legder', 'legint', 'leg2poly', 'poly2leg', 'legfromroots',
-    'legvander', 'legfit', 'legtrim', 'legroots', 'Legendre']
+    'legvander', 'legfit', 'legtrim', 'legroots', 'Legendre','legval2d',
+    'legval3d', 'leggrid2d', 'leggrid3d', 'legvander2d', 'legvander3d']
 
 legtrim = pu.trimcoef
 
@@ -1103,6 +1105,98 @@ def legvander(x, deg) :
     return np.rollaxis(v, 0, v.ndim)
 
 
+def legvander2d(x, y, deg) :
+    """Pseudo Vandermonde matrix of given degree.
+
+    Returns the pseudo Vandermonde matrix for 2D Legendre series in `x` and
+    `y`. The sample point coordinates must all have the same shape after
+    conversion to arrays and the dtype will be converted to either float64
+    or complex128 depending on whether any of `x` or 'y' are complex.  The
+    maximum degrees of the 2D Legendre series in each variable are specified in
+    the list `deg` in the form ``[xdeg, ydeg]``. The return array has the
+    shape ``x.shape + (order,)`` if `x`, and `y` are arrays or ``(1, order)
+    if they are scalars. Here order is the number of elements in a
+    flattened coefficient array of original shape ``(xdeg + 1, ydeg + 1)``.
+    The flattening is done so that the resulting pseudo Vandermonde array
+    can be easily used in least squares fits.
+
+    Parameters
+    ----------
+    x,y : array_like
+        Arrays of point coordinates, each of the same shape.
+    deg : list
+        List of maximum degrees of the form [x_deg, y_deg].
+
+    Returns
+    -------
+    vander2d : ndarray
+        The shape of the returned matrix is described above.
+
+    See Also
+    --------
+    legvander, legvander3d. legval2d, legval3d
+
+    """
+    ideg = [int(d) for d in deg]
+    is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
+    if is_valid != [1, 1]:
+        raise ValueError("degrees must be non-negative integers")
+    degx, degy = deg
+    x, y = np.array((x, y), copy=0) + 0.0
+
+    vx = legvander(x, degx)
+    vy = legvander(y, degy)
+    v = np.einsum("...i,...j->...ij", vx, vy)
+    return v.reshape(v.shape[:-2] + (-1,))
+
+
+def legvander3d(x, y, z, deg) :
+    """Psuedo Vandermonde matrix of given degree.
+
+    Returns the pseudo Vandermonde matrix for 3D Legendre series in `x`, `y`,
+    or `z`. The sample point coordinates must all have the same shape after
+    conversion to arrays and the dtype will be converted to either float64
+    or complex128 depending on whether any of `x`, `y`, or 'z' are complex.
+    The maximum degrees of the 3D Legendre series in each variable are
+    specified in the list `deg` in the form ``[xdeg, ydeg, zdeg]``. The
+    return array has the shape ``x.shape + (order,)`` if `x`, `y`, and `z`
+    are arrays or ``(1, order) if they are scalars. Here order is the
+    number of elements in a flattened coefficient array of original shape
+    ``(xdeg + 1, ydeg + 1, zdeg + 1)``.  The flattening is done so that the
+    resulting pseudo Vandermonde array can be easily used in least squares
+    fits.
+
+    Parameters
+    ----------
+    x,y,z : array_like
+        Arrays of point coordinates, each of the same shape.
+    deg : list
+        List of maximum degrees of the form [x_deg, y_deg, z_deg].
+
+    Returns
+    -------
+    vander3d : ndarray
+        The shape of the returned matrix is described above.
+
+    See Also
+    --------
+    legvander, legvander3d. legval2d, legval3d
+
+    """
+    ideg = [int(d) for d in deg]
+    is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
+    if is_valid != [1, 1, 1]:
+        raise ValueError("degrees must be non-negative integers")
+    degx, degy, degz = deg
+    x, y, z = np.array((x, y, z), copy=0) + 0.0
+
+    vx = legvander(x, deg_x)
+    vy = legvander(y, deg_y)
+    vz = legvander(z, deg_z)
+    v = np.einsum("...i,...j,...k->...ijk", vx, vy, vz)
+    return v.reshape(v.shape[:-3] + (-1,))
+
+
 def legfit(x, y, deg, rcond=None, full=False, w=None):
     """
     Least squares fit of Legendre series to data.
diff --git a/numpy/polynomial/polynomial.py b/numpy/polynomial/polynomial.py
index a0c9f903a..bf9c1a216 100644
--- a/numpy/polynomial/polynomial.py
+++ b/numpy/polynomial/polynomial.py
@@ -37,6 +37,8 @@ Misc Functions
 - `polyfromroots` -- create a polynomial with specified roots.
 - `polyroots` -- find the roots of a polynomial.
 - `polyvander` -- Vandermonde-like matrix for powers.
+- `polyvander2d` -- Vandermonde-like matrix for 2D power series.
+- `polyvander3d` -- Vandermonde-like matrix for 3D power series.
 - `polyfit` -- least-squares fit returning a polynomial.
 - `polytrim` -- trim leading coefficients from a polynomial.
 - `polyline` -- polynomial representing given straight line.
@@ -54,9 +56,9 @@ from __future__ import division
 
 __all__ = ['polyzero', 'polyone', 'polyx', 'polydomain', 'polyline',
     'polyadd', 'polysub', 'polymulx', 'polymul', 'polydiv', 'polypow',
-    'polyval', 'polyval2d', 'polyval3d', 'polygrid2d', 'polygrid3d',
-    'polyder', 'polyint', 'polyfromroots', 'polyvander', 'polyfit',
-    'polytrim', 'polyroots', 'Polynomial']
+    'polyval', 'polyder', 'polyint', 'polyfromroots', 'polyvander',
+    'polyfit', 'polytrim', 'polyroots', 'Polynomial','polyval2d',
+    'polyval3d', 'polygrid2d', 'polygrid3d', 'polyvander2d','polyvander3d']
 
 import numpy as np
 import numpy.linalg as la
@@ -891,6 +893,10 @@ def polyvander(x, deg) :
         The shape of the returned matrix is ``x.shape + (deg+1,)``. The last
         index is the degree.
 
+    See Also
+    --------
+    polyvander2d, polyvander3d
+
     """
     ideg = int(deg)
     if ideg != deg:
@@ -908,6 +914,98 @@ def polyvander(x, deg) :
     return np.rollaxis(v, 0, v.ndim)
 
 
+def polyvander2d(x, y, deg) :
+    """Pseudo Vandermonde matrix of given degree.
+
+    Returns the pseudo Vandermonde matrix for 2D polynomials in `x` and
+    `y`. The sample point coordinates must all have the same shape after
+    conversion to arrays and the dtype will be converted to either float64
+    or complex128 depending on whether any of `x` or 'y' are complex.  The
+    maximum degrees of the 2D polynomials in each variable are specified in
+    the list `deg` in the form ``[xdeg, ydeg]``. The return array has the
+    shape ``x.shape + (order,)`` if `x`, and `y` are arrays or ``(1, order)
+    if they are scalars. Here order is the number of elements in a
+    flattened coefficient array of original shape ``(xdeg + 1, ydeg + 1)``.
+    The flattening is done so that the resulting pseudo Vandermonde array
+    can be easily used in least squares fits.
+
+    Parameters
+    ----------
+    x,y : array_like
+        Arrays of point coordinates, each of the same shape.
+    deg : list
+        List of maximum degrees of the form [x_deg, y_deg].
+
+    Returns
+    -------
+    vander2d : ndarray
+        The shape of the returned matrix is described above.
+
+    See Also
+    --------
+    polyvander, polyvander3d. polyval2d, polyval3d
+
+    """
+    ideg = [int(d) for d in deg]
+    is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
+    if is_valid != [1, 1]:
+        raise ValueError("degrees must be non-negative integers")
+    degx, degy = deg
+    x, y = np.array((x, y), copy=0) + 0.0
+
+    vx = polyvander(x, degx)
+    vy = polyvander(y, degy)
+    v = np.einsum("...i,...j->...ij", vx, vy)
+    return v.reshape(v.shape[:-2] + (-1,))
+
+
+def polyvander3d(x, y, z, deg) :
+    """Psuedo Vandermonde matrix of given degree.
+
+    Returns the pseudo Vandermonde matrix for 3D polynomials in `x`, `y`,
+    or `z`. The sample point coordinates must all have the same shape after
+    conversion to arrays and the dtype will be converted to either float64
+    or complex128 depending on whether any of `x`, `y`, or 'z' are complex.
+    The maximum degrees of the 3D polynomials in each variable are
+    specified in the list `deg` in the form ``[xdeg, ydeg, zdeg]``. The
+    return array has the shape ``x.shape + (order,)`` if `x`, `y`, and `z`
+    are arrays or ``(1, order) if they are scalars. Here order is the
+    number of elements in a flattened coefficient array of original shape
+    ``(xdeg + 1, ydeg + 1, zdeg + 1)``.  The flattening is done so that the
+    resulting pseudo Vandermonde array can be easily used in least squares
+    fits.
+
+    Parameters
+    ----------
+    x,y,z : array_like
+        Arrays of point coordinates, each of the same shape.
+    deg : list
+        List of maximum degrees of the form [x_deg, y_deg, z_deg].
+
+    Returns
+    -------
+    vander3d : ndarray
+        The shape of the returned matrix is described above.
+
+    See Also
+    --------
+    polyvander, polyvander3d. polyval2d, polyval3d
+
+    """
+    ideg = [int(d) for d in deg]
+    is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
+    if is_valid != [1, 1, 1]:
+        raise ValueError("degrees must be non-negative integers")
+    degx, degy, degz = deg
+    x, y, z = np.array((x, y, z), copy=0) + 0.0
+
+    vx = polyvander(x, deg_x)
+    vy = polyvander(y, deg_y)
+    vz = polyvander(z, deg_z)
+    v = np.einsum("...i,...j,...k->...ijk", vx, vy, vz)
+    return v.reshape(v.shape[:-3] + (-1,))
+
+
 def polyfit(x, y, deg, rcond=None, full=False, w=None):
     """
     Least-squares fit of a polynomial to data.