diff --git a/python/mxnet/ndarray/numpy/_op.py b/python/mxnet/ndarray/numpy/_op.py
index ff0e48d47664..31b2a87a17de 100644
--- a/python/mxnet/ndarray/numpy/_op.py
+++ b/python/mxnet/ndarray/numpy/_op.py
@@ -45,7 +45,7 @@
            'hypot', 'bitwise_and', 'bitwise_xor', 'bitwise_or', 'rad2deg', 'deg2rad', 'unique', 'lcm',
            'tril', 'identity', 'take', 'ldexp', 'vdot', 'inner', 'outer', 'kron',
            'equal', 'not_equal', 'greater', 'less', 'greater_equal', 'less_equal', 'roll', 'rot90', 'einsum',
-           'true_divide', 'nonzero', 'quantile', 'percentile', 'shares_memory', 'may_share_memory',
+           'true_divide', 'nonzero', 'quantile', 'percentile', 'shares_memory', 'may_share_memory', 'interp',
            'diff', 'ediff1d', 'resize', 'polyval', 'nan_to_num', 'isnan', 'isinf', 'isposinf', 'isneginf', 'isfinite',
            'where', 'bincount', 'pad', 'cumsum', 'diag', 'diagonal']
 
@@ -7081,6 +7081,86 @@ def may_share_memory(a, b, max_work=None):
     return _api_internal.share_memory(a, b).item()
 
 
+@set_module('mxnet.ndarray.numpy')
+def interp(x, xp, fp, left=None, right=None, period=None):  # pylint: disable=too-many-arguments
+    """
+    One-dimensional linear interpolation.
+    Returns the one-dimensional piecewise linear interpolant to a function
+    with given values at discrete data-points.
+
+    Parameters
+    ----------
+    x : ndarray
+        The x-coordinates of the interpolated values.
+    xp : 1-D array of floats
+        The x-coordinates of the data points, must be increasing if argument
+        `period` is not specified. Otherwise, `xp` is internally sorted after
+        normalizing the periodic boundaries with ``xp = xp % period``.
+    fp : 1-D array of floats
+        The y-coordinates of the data points, same length as `xp`.
+    left : optional float corresponding to fp
+        Value to return for `x < xp[0]`, default is `fp[0]`.
+    right : optional float corresponding to fp
+        Value to return for `x > xp[-1]`, default is `fp[-1]`.
+    period : None or float, optional
+        A period for the x-coordinates. This parameter allows the proper
+        interpolation of angular x-coordinates. Parameters `left` and `right`
+        are ignored if `period` is specified.
+        .. versionadded:: 1.10.0
+
+    Returns
+    -------
+    y : float (corresponding to fp) or ndarray
+        The interpolated values, same shape as `x`.
+    Raises
+    ------
+    ValueError
+        If `xp` and `fp` have different length
+        If `xp` or `fp` are not 1-D sequences
+        If `period == 0`
+
+    Notes
+    -----
+    Does not check that the x-coordinate sequence `xp` is increasing.
+    If `xp` is not increasing, the results are nonsense.
+    A simple check for increasing is::
+        np.all(np.diff(xp) > 0)
+
+    Examples
+    --------
+    >>> xp = [1, 2, 3]
+    >>> fp = [3, 2, 0]
+    >>> np.interp(2.5, xp, fp)
+    1.0
+    >>> np.interp([0, 1, 1.5, 2.72, 3.14], xp, fp)
+    array([ 3. ,  3. ,  2.5 ,  0.56,  0. ])
+    >>> UNDEF = -99.0
+    >>> np.interp(3.14, xp, fp, right=UNDEF)
+    -99.0
+    Plot an interpolant to the sine function:
+    >>> x = np.linspace(0, 2*np.pi, 10)
+    >>> y = np.sin(x)
+    >>> xvals = np.linspace(0, 2*np.pi, 50)
+    >>> yinterp = np.interp(xvals, x, y)
+    >>> import matplotlib.pyplot as plt
+    >>> plt.plot(x, y, 'o')
+    [<matplotlib.lines.Line2D object at 0x...>]
+    >>> plt.plot(xvals, yinterp, '-x')
+    [<matplotlib.lines.Line2D object at 0x...>]
+    >>> plt.show()
+    Interpolation with periodic x-coordinates:
+    >>> x = [-180, -170, -185, 185, -10, -5, 0, 365]
+    >>> xp = [190, -190, 350, -350]
+    >>> fp = [5, 10, 3, 4]
+    >>> np.interp(x, xp, fp, period=360)
+    array([7.5, 5., 8.75, 6.25, 3., 3.25, 3.5, 3.75])
+    """
+    if not isinstance(x, numeric_types):
+        x = x.astype(float)
+    return _api_internal.interp(xp.astype(float), fp.astype(float), x, left,
+                                right, period)
+
+
 @set_module('mxnet.ndarray.numpy')
 def diff(a, n=1, axis=-1, prepend=None, append=None):  # pylint: disable=redefined-outer-name
     r"""
diff --git a/python/mxnet/numpy/fallback.py b/python/mxnet/numpy/fallback.py
index 54df715b60c0..ca221874346f 100644
--- a/python/mxnet/numpy/fallback.py
+++ b/python/mxnet/numpy/fallback.py
@@ -50,7 +50,6 @@
     'histogramdd',
     'i0',
     'in1d',
-    'interp',
     'intersect1d',
     'isclose',
     'isin',
@@ -138,7 +137,6 @@
 histogramdd = onp.histogramdd
 i0 = onp.i0
 in1d = onp.in1d
-interp = onp.interp
 intersect1d = onp.intersect1d
 isclose = onp.isclose
 isin = onp.isin
diff --git a/python/mxnet/numpy/multiarray.py b/python/mxnet/numpy/multiarray.py
index 281a6f7cc3fc..73207da741e6 100644
--- a/python/mxnet/numpy/multiarray.py
+++ b/python/mxnet/numpy/multiarray.py
@@ -69,7 +69,7 @@
            'flip', 'flipud', 'fliplr', 'around', 'round', 'round_', 'arctan2', 'hypot',
            'bitwise_and', 'bitwise_xor', 'bitwise_or', 'rad2deg', 'deg2rad',
            'unique', 'lcm', 'tril', 'identity', 'take', 'ldexp', 'vdot', 'inner', 'outer', 'kron',
-           'equal', 'not_equal',
+           'equal', 'not_equal', 'interp',
            'greater', 'less', 'greater_equal', 'less_equal', 'roll', 'rot90', 'einsum', 'true_divide', 'nonzero',
            'quantile', 'percentile', 'shares_memory', 'may_share_memory', 'diff', 'ediff1d', 'resize', 'matmul',
            'nan_to_num', 'isnan', 'isinf', 'isposinf', 'isneginf', 'isfinite', 'polyval', 'where', 'bincount',
@@ -9275,6 +9275,83 @@ def resize(a, new_shape):
     return _mx_nd_np.resize(a, new_shape)
 
 
+@set_module('mxnet.numpy')
+def interp(x, xp, fp, left=None, right=None, period=None):  # pylint: disable=too-many-arguments
+    """
+    One-dimensional linear interpolation.
+    Returns the one-dimensional piecewise linear interpolant to a function
+    with given values at discrete data-points.
+
+    Parameters
+    ----------
+    x : ndarray
+        The x-coordinates of the interpolated values.
+    xp : 1-D array of floats
+        The x-coordinates of the data points, must be increasing if argument
+        `period` is not specified. Otherwise, `xp` is internally sorted after
+        normalizing the periodic boundaries with ``xp = xp % period``.
+    fp : 1-D array of floats
+        The y-coordinates of the data points, same length as `xp`.
+    left : optional float corresponding to fp
+        Value to return for `x < xp[0]`, default is `fp[0]`.
+    right : optional float corresponding to fp
+        Value to return for `x > xp[-1]`, default is `fp[-1]`.
+    period : None or float, optional
+        A period for the x-coordinates. This parameter allows the proper
+        interpolation of angular x-coordinates. Parameters `left` and `right`
+        are ignored if `period` is specified.
+        .. versionadded:: 1.10.0
+
+    Returns
+    -------
+    y : float (corresponding to fp) or ndarray
+        The interpolated values, same shape as `x`.
+    Raises
+    ------
+    ValueError
+        If `xp` and `fp` have different length
+        If `xp` or `fp` are not 1-D sequences
+        If `period == 0`
+
+    Notes
+    -----
+    Does not check that the x-coordinate sequence `xp` is increasing.
+    If `xp` is not increasing, the results are nonsense.
+    A simple check for increasing is::
+        np.all(np.diff(xp) > 0)
+
+    Examples
+    --------
+    >>> xp = [1, 2, 3]
+    >>> fp = [3, 2, 0]
+    >>> np.interp(2.5, xp, fp)
+    1.0
+    >>> np.interp([0, 1, 1.5, 2.72, 3.14], xp, fp)
+    array([ 3. ,  3. ,  2.5 ,  0.56,  0. ])
+    >>> UNDEF = -99.0
+    >>> np.interp(3.14, xp, fp, right=UNDEF)
+    -99.0
+    Plot an interpolant to the sine function:
+    >>> x = np.linspace(0, 2*np.pi, 10)
+    >>> y = np.sin(x)
+    >>> xvals = np.linspace(0, 2*np.pi, 50)
+    >>> yinterp = np.interp(xvals, x, y)
+    >>> import matplotlib.pyplot as plt
+    >>> plt.plot(x, y, 'o')
+    [<matplotlib.lines.Line2D object at 0x...>]
+    >>> plt.plot(xvals, yinterp, '-x')
+    [<matplotlib.lines.Line2D object at 0x...>]
+    >>> plt.show()
+    Interpolation with periodic x-coordinates:
+    >>> x = [-180, -170, -185, 185, -10, -5, 0, 365]
+    >>> xp = [190, -190, 350, -350]
+    >>> fp = [5, 10, 3, 4]
+    >>> np.interp(x, xp, fp, period=360)
+    array([7.5, 5., 8.75, 6.25, 3., 3.25, 3.5, 3.75])
+    """
+    return _mx_nd_np.interp(x, xp, fp, left=left, right=right, period=period)
+
+
 # pylint: disable=redefined-outer-name
 @set_module('mxnet.numpy')
 def full_like(a, fill_value, dtype=None, order='C', ctx=None, out=None): # pylint: disable=too-many-arguments
diff --git a/python/mxnet/numpy_dispatch_protocol.py b/python/mxnet/numpy_dispatch_protocol.py
index 110f2273a852..8a1beba8e87c 100644
--- a/python/mxnet/numpy_dispatch_protocol.py
+++ b/python/mxnet/numpy_dispatch_protocol.py
@@ -110,6 +110,7 @@ def _run_with_array_ufunc_proto(*args, **kwargs):
     'fliplr',
     'inner',
     'insert',
+    'interp',
     'max',
     'amax',
     'mean',
diff --git a/python/mxnet/symbol/numpy/_symbol.py b/python/mxnet/symbol/numpy/_symbol.py
index a2a4cd9d3584..d204c3e56967 100644
--- a/python/mxnet/symbol/numpy/_symbol.py
+++ b/python/mxnet/symbol/numpy/_symbol.py
@@ -48,7 +48,7 @@
            'flatnonzero',
            'swapaxes', 'clip', 'argmax', 'argmin', 'std', 'var', 'indices', 'copysign', 'ravel', 'unravel_index',
            'diag_indices_from', 'hanning', 'hamming', 'blackman', 'flip', 'flipud', 'fliplr',
-           'hypot', 'bitwise_and', 'bitwise_xor', 'bitwise_or', 'rad2deg', 'deg2rad', 'unique', 'lcm',
+           'hypot', 'bitwise_and', 'bitwise_xor', 'bitwise_or', 'rad2deg', 'deg2rad', 'unique', 'lcm', 'interp',
            'tril', 'identity', 'take', 'ldexp', 'vdot', 'inner', 'outer', 'kron',
            'equal', 'not_equal', 'greater', 'less', 'greater_equal', 'less_equal', 'roll', 'rot90', 'einsum',
            'true_divide', 'quantile', 'percentile', 'shares_memory', 'may_share_memory', 'diff', 'ediff1d',
@@ -6349,6 +6349,59 @@ def ediff1d(ary, to_end=None, to_begin=None):
                             to_begin_scalar=to_begin, to_end_scalar=to_end)
 
 
+@set_module('mxnet.symbol.numpy')
+def interp(x, xp, fp, left=None, right=None, period=None):  # pylint: disable=too-many-arguments
+    """
+    One-dimensional linear interpolation.
+    Returns the one-dimensional piecewise linear interpolant to a function
+    with given values at discrete data-points.
+
+    Parameters
+    ----------
+    x : _Symbol
+        The x-coordinates of the interpolated values.
+    xp : _Symbol
+        The x-coordinates of the data points, must be increasing if argument
+        `period` is not specified. Otherwise, `xp` is internally sorted after
+        normalizing the periodic boundaries with ``xp = xp % period``.
+    fp : _Symbol
+        The y-coordinates of the data points, same length as `xp`.
+    left : optional float corresponding to fp
+        Value to return for `x < xp[0]`, default is `fp[0]`.
+    right : optional float corresponding to fp
+        Value to return for `x > xp[-1]`, default is `fp[-1]`.
+    period : None or float, optional
+        A period for the x-coordinates. This parameter allows the proper
+        interpolation of angular x-coordinates. Parameters `left` and `right`
+        are ignored if `period` is specified.
+        .. versionadded:: 1.10.0
+
+    Returns
+    -------
+    y : _Symbol
+        The interpolated values, same shape as `x`.
+
+    Raises
+    ------
+    ValueError
+        If `xp` and `fp` have different length
+        If `xp` or `fp` are not 1-D sequences
+        If `period == 0`
+
+    Notes
+    -----
+    Does not check that the x-coordinate sequence `xp` is increasing.
+    If `xp` is not increasing, the results are nonsense.
+    A simple check for increasing is::
+        np.all(np.diff(xp) > 0)
+    """
+    if isinstance(x, numeric_types):
+        return _npi.interp(xp.astype(float), fp.astype(float), left=left,
+                           right=right, period=period, x_scalar=x, x_is_scalar=True)
+    return _npi.interp(xp.astype(float), fp.astype(float), x.astype(float), left=left,
+                       right=right, period=period, x_scalar=0.0, x_is_scalar=False)
+
+
 @set_module('mxnet.symbol.numpy')
 def resize(a, new_shape):
     """
diff --git a/src/api/operator/numpy/np_interp_op.cc b/src/api/operator/numpy/np_interp_op.cc
new file mode 100644
index 000000000000..7959383c9230
--- /dev/null
+++ b/src/api/operator/numpy/np_interp_op.cc
@@ -0,0 +1,78 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ * 
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ * 
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+/*!
+ * \file np_interp_op.cc
+ * \brief Implementation of the API of functions in src/operator/numpy/np_interp_op.cc
+ */
+#include <mxnet/api_registry.h>
+#include <mxnet/runtime/packed_func.h>
+#include "../utils.h"
+#include "../../../operator/numpy/np_interp_op-inl.h"
+
+namespace mxnet {
+
+MXNET_REGISTER_API("_npi.interp")
+.set_body([](runtime::MXNetArgs args, runtime::MXNetRetValue* ret) {
+  using namespace runtime;
+  static const nnvm::Op* op = Op::Get("_npi_interp");
+  nnvm::NodeAttrs attrs;
+  op::NumpyInterpParam param;
+  if (args[3].type_code() == kNull) {
+    param.left = dmlc::nullopt;
+  } else {
+    param.left = args[3].operator double();
+  }
+  if (args[4].type_code() == kNull) {
+    param.right = dmlc::nullopt;
+  } else {
+    param.right = args[4].operator double();
+  }
+  if (args[5].type_code() == kNull) {
+    param.period = dmlc::nullopt;
+  } else {
+    param.period = args[5].operator double();
+  }
+  if (args[2].type_code() == kDLInt || args[2].type_code() == kDLFloat) {
+    param.x_scalar = args[2].operator double();
+    param.x_is_scalar = true;
+    attrs.op = op;
+    attrs.parsed = std::move(param);
+    SetAttrDict<op::NumpyInterpParam>(&attrs);
+    NDArray* inputs[] = {args[0].operator mxnet::NDArray*(), args[1].operator mxnet::NDArray*()};
+    int num_inputs = 2;
+    int num_outputs = 0;
+    auto ndoutputs = Invoke(op, &attrs, num_inputs, inputs, &num_outputs, nullptr);
+    *ret = reinterpret_cast<mxnet::NDArray*>(ndoutputs[0]);
+  } else {
+    param.x_scalar = 0.0;
+    param.x_is_scalar = false;
+    attrs.op = op;
+    attrs.parsed = std::move(param);
+    SetAttrDict<op::NumpyInterpParam>(&attrs);
+    NDArray* inputs[] = {args[0].operator mxnet::NDArray*(), args[1].operator mxnet::NDArray*(),
+                         args[2].operator mxnet::NDArray*()};
+    int num_inputs = 3;
+    int num_outputs = 0;
+    auto ndoutputs = Invoke(op, &attrs, num_inputs, inputs, &num_outputs, nullptr);
+    *ret = reinterpret_cast<mxnet::NDArray*>(ndoutputs[0]);
+  }
+});
+
+}  // namespace mxnet
diff --git a/src/operator/numpy/np_interp_op-inl.h b/src/operator/numpy/np_interp_op-inl.h
new file mode 100644
index 000000000000..c0d595699a82
--- /dev/null
+++ b/src/operator/numpy/np_interp_op-inl.h
@@ -0,0 +1,286 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+/*!
+ * Copyright (c) 2020 by Contributors
+ * \file np_interp_op-inl.h
+*/
+
+#ifndef MXNET_OPERATOR_NUMPY_NP_INTERP_OP_INL_H_
+#define MXNET_OPERATOR_NUMPY_NP_INTERP_OP_INL_H_
+
+#include <vector>
+#include <string>
+#include <unordered_map>
+#include "../tensor/ordering_op-inl.h"
+#include "../tensor/matrix_op-inl.h"
+#include "../tensor/elemwise_binary_scalar_op.h"
+#include "../../common/utils.h"
+#include "../mshadow_op.h"
+#include "../operator_common.h"
+#include "../elemwise_op_common.h"
+#include "np_broadcast_reduce_op.h"
+
+namespace mxnet {
+namespace op {
+
+struct NumpyInterpParam : public dmlc::Parameter<NumpyInterpParam> {
+  dmlc::optional<double> left;
+  dmlc::optional<double> right;
+  dmlc::optional<double> period;
+  double x_scalar;
+  bool x_is_scalar;
+  DMLC_DECLARE_PARAMETER(NumpyInterpParam) {
+    DMLC_DECLARE_FIELD(left)
+      .set_default(dmlc::optional<double>())
+      .describe("Value to return for x < xp[0], default is fp[0].");
+    DMLC_DECLARE_FIELD(right)
+      .set_default(dmlc::optional<double>())
+      .describe("Value to return for x > xp[-1], default is fp[-1].");
+    DMLC_DECLARE_FIELD(period)
+      .set_default(dmlc::optional<double>())
+      .describe("A period for the x-coordinates. This parameter allows"
+                "the proper interpolation of angular x-coordinates. Parameters"
+                "left and right are ignored if period is specified.");
+    DMLC_DECLARE_FIELD(x_scalar).set_default(0.0)
+      .describe("x is a scalar input");
+    DMLC_DECLARE_FIELD(x_is_scalar).set_default(false)
+      .describe("Flag that determines whether input is a scalar");
+  }
+  void SetAttrDict(std::unordered_map<std::string, std::string>* dict) {
+    std::ostringstream left_s, right_s, period_s, x_scalar_s, x_is_scalar_s;
+    left_s << left;
+    right_s << right;
+    period_s << period;
+    x_scalar_s << x_scalar;
+    x_is_scalar_s << x_is_scalar;
+    (*dict)["left"] = left_s.str();
+    (*dict)["right"] = right_s.str();
+    (*dict)["period"] = period_s.str();
+    (*dict)["x_scalar"] = x_scalar_s.str();
+    (*dict)["x_is_scalar"] = x_is_scalar_s.str();
+  }
+};
+
+struct interp {
+  MSHADOW_XINLINE static void Map(int i,
+                                  double* out,
+                                  const double* x,
+                                  const double* xp,
+                                  const double* fp,
+                                  const int dsize,
+                                  const double left,
+                                  const double right,
+                                  const bool has_left,
+                                  const bool has_right) {
+    double x_value = x[i];
+    double xp_low = xp[0];
+    double xp_above = xp[dsize-1];
+    double lval = has_left ? left : fp[0];
+    double rval = has_right ? right : fp[dsize-1];
+
+    if (x_value > xp_above) {
+      out[i] = rval;
+    } else if (x_value < xp_low) {
+      out[i] = lval;
+    } else {
+      int imin = 0;
+      int imax = dsize;
+      int imid;
+      while (imin < imax) {
+        imid = static_cast<int>((imax + imin) / 2);
+        if (x_value >= xp[imid]) {
+          imin = imid + 1;
+        } else {
+          imax = imid;
+        }
+      }  // biserction search
+
+      int j = imin;
+      if (j == dsize) {
+        out[i] = fp[dsize-1];
+      } else if (x_value == xp[j-1]) {
+        out[i] = fp[j-1];  // void potential non-finite interpolation
+      } else {
+        double xp_below = xp[j-1];
+        double xp_above = xp[j];
+        double weight_above = (x_value - xp_below) / (xp_above - xp_below);
+        double weigth_below = 1 - weight_above;
+        double x1 = fp[j-1] * weigth_below;
+        double x2 = fp[j] * weight_above;
+        out[i] = x1 + x2;
+      }
+    }
+  }
+};
+
+struct interp_period {
+  MSHADOW_XINLINE static void Map(int i,
+                                  double* out,
+                                  const double* x,
+                                  const double* xp,
+                                  const double* fp,
+                                  const index_t* idx,
+                                  const int dsize,
+                                  const double period) {
+    double x_value = x[i];
+    int imin = 0;
+    int imax = dsize;
+    int imid;
+    while (imin < imax) {
+      imid = static_cast<int>((imax + imin) / 2);
+      if (x_value >= xp[idx[imid]]) {
+        imin = imid + 1;
+      } else {
+        imax = imid;
+      }
+    }  // biserction search
+
+    int j = imin;
+    double xp_below, xp_above;
+    double fp1, fp2;
+    if (j == 0) {
+      xp_below = xp[idx[dsize-1]] - period;
+      xp_above = xp[idx[0]];
+      fp1 = fp[idx[dsize-1]];
+      fp2 = fp[idx[0]];
+    } else if (j == dsize) {
+      xp_below = xp[idx[dsize-1]];
+      xp_above = xp[idx[0]] + period;
+      fp1 = fp[idx[dsize-1]];
+      fp2 = fp[idx[0]];
+    } else {
+      xp_below = xp[idx[j-1]];
+      xp_above = xp[idx[j]];
+      fp1 = fp[idx[j-1]];
+      fp2 = fp[idx[j]];
+    }
+    double weight_above = (x_value - xp_below) / (xp_above - xp_below);
+    double weigth_below = 1 - weight_above;
+    double x1 = fp1 * weigth_below;
+    double x2 = fp2 * weight_above;
+    out[i] = x1 + x2;
+  }
+};
+
+template<typename xpu, typename OP>
+void NumpyInterpForward(const nnvm::NodeAttrs& attrs,
+                        const OpContext &ctx,
+                        const std::vector<TBlob> &inputs,
+                        const std::vector<OpReqType> &req,
+                        const std::vector<TBlob> &outputs) {
+  if (req[0] == kNullOp) return;
+  using namespace mxnet;
+  using namespace mxnet_op;
+  using namespace mshadow;
+  using namespace mshadow::expr;
+  CHECK_GE(inputs.size(), 2U);
+  CHECK_EQ(outputs.size(), 1U);
+
+  Stream<xpu> *s = ctx.get_stream<xpu>();
+  const NumpyInterpParam& param = nnvm::get<NumpyInterpParam>(attrs.parsed);
+  dmlc::optional<double> left = param.left;
+  dmlc::optional<double> right = param.right;
+  dmlc::optional<double> period = param.period;
+  bool x_is_scalar = param.x_is_scalar;
+
+  TBlob xp = inputs[0];
+  const TBlob &fp = inputs[1];
+  const TBlob &out = outputs[0];
+  bool has_left = left.has_value() ? true : false;
+  bool has_right = right.has_value() ? true : false;
+  double left_value = left.has_value() ? left.value() : 0.0;
+  double right_value = right.has_value() ? right.value() : 0.0;
+
+  CHECK_GE(xp.Size(), 1U) <<"ValueError: array of sample points is empty";
+
+  TopKParam topk_param = TopKParam();
+  topk_param.axis = dmlc::optional<int>(-1);
+  topk_param.is_ascend = true;
+  topk_param.k = 0;
+  topk_param.ret_typ = topk_enum::kReturnIndices;
+
+  size_t topk_temp_size;  // Used by Sort
+  size_t topk_workspace_size = TopKWorkspaceSize<xpu, double>(xp, topk_param, &topk_temp_size);
+  size_t size_x = x_is_scalar ? 8 : 0;
+  size_t size_norm_x = x_is_scalar ? 8 : inputs[2].Size() * sizeof(double);
+  size_t size_norm_xp = xp.Size() * sizeof(double);
+  size_t size_norm = period.has_value()? size_norm_x + size_norm_xp : 0;
+  size_t size_idx = period.has_value()? xp.Size() * sizeof(index_t) : 0;
+  size_t workspace_size =
+    topk_workspace_size + size_x + size_norm + size_idx;
+
+  Tensor<xpu, 1, char> temp_mem =
+    ctx.requested[0].get_space_typed<xpu, 1, char>(Shape1(workspace_size), s);
+
+  char* workspace_curr_ptr = temp_mem.dptr_;
+
+  TBlob x, idx;
+  if (x_is_scalar) {
+    double x_scalar = param.x_scalar;
+    Tensor<cpu, 1, double> host_x(&x_scalar, Shape1(1), ctx.get_stream<cpu>());
+    Tensor<xpu, 1, double> device_x(reinterpret_cast<double*>(workspace_curr_ptr),
+                                    Shape1(1), ctx.get_stream<xpu>());
+    Copy(device_x, host_x, ctx.get_stream<xpu>());
+    x = TBlob(device_x.dptr_, TShape(0, 1), xpu::kDevMask);
+    workspace_curr_ptr += 8;
+  } else {
+    x = inputs[2];
+  }   // handle input x is a scalar
+
+  // normalize the input data by periodic boundaries.
+  if (period.has_value()) {
+    double* norm_xp_ptr;
+    double* norm_x_ptr;
+    double period_value = period.value();
+    index_t* idx_ptr;
+    CHECK_NE(period_value, 0.0)<< "period must be a non-zero value";
+
+    norm_xp_ptr = reinterpret_cast<double*>(workspace_curr_ptr);
+    norm_x_ptr = reinterpret_cast<double*>(workspace_curr_ptr + size_norm_xp);
+    idx_ptr = reinterpret_cast<index_t*>(workspace_curr_ptr + size_norm_xp + size_norm_x);
+
+    TBlob norm_x = TBlob(norm_x_ptr, x.shape_, xpu::kDevMask);
+    TBlob norm_xp = TBlob(norm_xp_ptr, xp.shape_, xpu::kDevMask);
+    const OpReqType ReqType = kWriteTo;
+    Kernel<op_with_req<OP, ReqType>, xpu>::Launch(
+      s, x.Size(), norm_x.dptr<double>(), x.dptr<double>(), period_value);
+    Kernel<op_with_req<OP, ReqType>, xpu>::Launch(
+      s, xp.Size(), norm_xp.dptr<double>(), xp.dptr<double>(), period_value);
+
+    workspace_curr_ptr += size_x + size_norm + size_idx;
+    idx = TBlob(idx_ptr, xp.shape_, xpu::kDevMask);
+    std::vector<OpReqType> req_TopK = {kWriteTo};
+    std::vector<TBlob> ret = {idx};
+
+    TopKImplwithWorkspace<xpu, double, index_t>(ctx.run_ctx, req_TopK, norm_xp, ret, topk_param,
+                                                workspace_curr_ptr, topk_temp_size, s);
+    Kernel<interp_period, xpu>::Launch(
+      s, norm_x.Size(), out.dptr<double>(), norm_x.dptr<double>(), norm_xp.dptr<double>(),
+      fp.dptr<double>(), idx.dptr<index_t>(), norm_xp.Size(), period_value);
+  } else {
+    Kernel<interp, xpu>::Launch(
+      s, x.Size(), out.dptr<double>(), x.dptr<double>(), xp.dptr<double>(), fp.dptr<double>(),
+      xp.Size(), left_value, right_value, has_left, has_right);
+  }
+}
+
+}  // namespace op
+}  // namespace mxnet
+
+#endif  // MXNET_OPERATOR_NUMPY_NP_INTERP_OP_INL_H_
diff --git a/src/operator/numpy/np_interp_op.cc b/src/operator/numpy/np_interp_op.cc
new file mode 100644
index 000000000000..af4bb1cd0bb3
--- /dev/null
+++ b/src/operator/numpy/np_interp_op.cc
@@ -0,0 +1,92 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+/*!
+ * Copyright (c) 2020 by Contributors
+ * \file np_interp_op.cc
+ * \brief CPU Implementation of Numpy-compatible interp
+*/
+
+#include "np_interp_op-inl.h"
+
+namespace mxnet {
+namespace op {
+
+inline bool NumpyInterpShape(const nnvm::NodeAttrs& attrs,
+                             std::vector<TShape> *in_attrs,
+                             std::vector<TShape> *out_attrs) {
+  CHECK_GE(in_attrs->size(), 2U);
+  CHECK_EQ(out_attrs->size(), 1U);
+  const NumpyInterpParam& param = nnvm::get<NumpyInterpParam>(attrs.parsed);
+
+  TShape oshape;
+  CHECK_EQ(in_attrs->at(0).ndim(), 1U)
+    << "ValueError: Data points must be 1-D array";
+  CHECK_EQ(in_attrs->at(1).ndim(), 1U)
+    << "ValueError: Data points must be 1-D array";
+  CHECK_EQ(in_attrs->at(0)[0], in_attrs->at(1)[0])
+    << "ValueError: fp and xp are not of the same length";
+  oshape = param.x_is_scalar ? TShape(0, 1) : in_attrs->at(2);
+  SHAPE_ASSIGN_CHECK(*out_attrs, 0, oshape);
+  return shape_is_known(out_attrs->at(0));
+}
+
+inline bool NumpyInterpType(const nnvm::NodeAttrs& attrs,
+                            std::vector<int> *in_attrs,
+                            std::vector<int> *out_attrs) {
+  CHECK_GE(in_attrs->size(), 2U);
+  CHECK_EQ(out_attrs->size(), 1U);
+
+  TYPE_ASSIGN_CHECK(*out_attrs, 0, mshadow::kFloat64);
+  return out_attrs->at(0) != -1;
+}
+
+DMLC_REGISTER_PARAMETER(NumpyInterpParam);
+
+NNVM_REGISTER_OP(_npi_interp)
+.set_num_inputs([](const NodeAttrs& attrs) {
+  const NumpyInterpParam& param =
+    nnvm::get<NumpyInterpParam>(attrs.parsed);
+  return param.x_is_scalar ? 2 : 3;
+  })
+.set_num_outputs(1)
+.set_attr_parser(ParamParser<NumpyInterpParam>)
+.set_attr<mxnet::FInferShape>("FInferShape", NumpyInterpShape)
+.set_attr<nnvm::FInferType>("FInferType", NumpyInterpType)
+.set_attr<nnvm::FListInputNames>("FListInputNames",
+  [](const NodeAttrs& attrs) {
+    const NumpyInterpParam& param =
+      nnvm::get<NumpyInterpParam>(attrs.parsed);
+    return param.x_is_scalar ?
+           std::vector<std::string>{"xp", "fp"} :
+           std::vector<std::string>{"xp", "fp", "x"};
+  })
+.set_attr<FCompute>("FCompute<cpu>", NumpyInterpForward<cpu, mshadow_op::mod>)
+.set_attr<FResourceRequest>("FResourceRequest",
+  [](const NodeAttrs& attrs) {
+    return std::vector<ResourceRequest>{ResourceRequest::kTempSpace};
+  })
+.set_attr<THasDeterministicOutput>("THasDeterministicOutput", true)
+.set_attr<nnvm::FGradient>("FGradient", MakeZeroGradNodes)
+.add_argument("xp", "NDArray-or-Symbol", "Input x-coordinates")
+.add_argument("fp", "NDArray-or-Symbol", "Input y-coordinates")
+.add_argument("x", "NDArray-or-Symbol", "Input data")
+.add_arguments(NumpyInterpParam::__FIELDS__());
+
+}  // namespace op
+}  // namespace mxnet
diff --git a/src/operator/numpy/np_interp_op.cu b/src/operator/numpy/np_interp_op.cu
new file mode 100644
index 000000000000..d3753e3eab0f
--- /dev/null
+++ b/src/operator/numpy/np_interp_op.cu
@@ -0,0 +1,34 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+/*!
+ * Copyright (c) 2020 by Contributors
+ * \file np_interp_op.cu
+ * \brief GPU Implementation of Numpy-compatible interp
+*/
+
+#include "np_interp_op-inl.h"
+
+namespace mxnet {
+namespace op {
+
+NNVM_REGISTER_OP(_npi_interp)
+.set_attr<FCompute>("FCompute<gpu>", NumpyInterpForward<gpu, mshadow_op::mod>);
+
+}  // namespace op
+}  // namespace mxnet
diff --git a/tests/python/unittest/test_numpy_interoperability.py b/tests/python/unittest/test_numpy_interoperability.py
index a4492a3beab1..6a1a832e985a 100644
--- a/tests/python/unittest/test_numpy_interoperability.py
+++ b/tests/python/unittest/test_numpy_interoperability.py
@@ -1325,6 +1325,38 @@ def _add_workload_insert():
     OpArgMngr.add_workload('insert', np.array(1), 0, np.array([0]))
 
 
+def _add_workload_interp():
+    xp0 = np.linspace(0, 1, 5)
+    fp0 = np.linspace(0, 1, 5)
+    x0 = np.linspace(0, 1, 50)
+    xp1 = np.array([1, 2, 3, 4])
+    fp1 = np.array([1, 2, np.inf, 4])
+    x1 = np.array([1, 2, 2.5, 3, 4])
+    xp2 = np.arange(0, 10, 0.0001)
+    fp2 = np.sin(xp2)
+    xp3 = np.array([190, -190, 350, -350])
+    fp3 = np.array([5, 10, 3, 4])
+    x3 = np.array([-180, -170, -185, 185, -10, -5, 0, 365])
+
+    OpArgMngr.add_workload('interp', x0, xp0, fp0)
+    OpArgMngr.add_workload('interp', x1, xp1, fp1)
+    OpArgMngr.add_workload('interp', np.pi, xp2, fp2)
+    OpArgMngr.add_workload('interp', x3, xp3, fp3, period=360)
+    for size in range(1, 10):
+        xp = np.arange(size, dtype=np.float64)
+        fp = np.ones(size, dtype=np.float64)
+        incpts = np.array([-1, 0, size - 1, size], dtype=np.float64)
+        decpts = incpts[::-1]
+        OpArgMngr.add_workload('interp', incpts, xp, fp)
+        OpArgMngr.add_workload('interp', decpts, xp, fp)
+        OpArgMngr.add_workload('interp', incpts, xp, fp, left=0)
+        OpArgMngr.add_workload('interp', decpts, xp, fp, left=0)
+        OpArgMngr.add_workload('interp', incpts, xp, fp, right=2)
+        OpArgMngr.add_workload('interp', decpts, xp, fp, right=2)
+        OpArgMngr.add_workload('interp', incpts, xp, fp, left=0, right=2)
+        OpArgMngr.add_workload('interp', decpts, xp, fp, left=0, right=2)
+
+
 def _add_workload_hypot():
     OpArgMngr.add_workload('hypot', np.array(1), np.array(1))
     OpArgMngr.add_workload('hypot', np.array(0), np.array(0))
@@ -2852,6 +2884,7 @@ def _prepare_workloads():
     _add_workload_true_divide()
     _add_workload_inner()
     _add_workload_insert()
+    _add_workload_interp()
     _add_workload_hypot()
     _add_workload_lcm()
     _add_workload_bitwise_and()
diff --git a/tests/python/unittest/test_numpy_op.py b/tests/python/unittest/test_numpy_op.py
index d3964af0bab7..3ccdc3afc972 100644
--- a/tests/python/unittest/test_numpy_op.py
+++ b/tests/python/unittest/test_numpy_op.py
@@ -8358,6 +8358,73 @@ def hybrid_forward(self, F, a):
             assert_almost_equal(elem_mx.asnumpy(), elem_np, rtol=rtol, atol=atol)
 
 
+@with_seed()
+@use_np
+def test_np_interp():
+    class TestInterp(HybridBlock):
+        def __init__(self, left=None, right=None, period=None):
+            super(TestInterp, self).__init__()
+            self._left = left
+            self._right = right
+            self._period = period
+        
+        def hybrid_forward(self, F, x, xp, fp):
+            return F.np.interp(x, xp, fp, left=self._left, right=self._right, period=self._period)
+    
+    class TestInterpScalar(HybridBlock):
+        def __init__(self, x=None, left=None, right=None, period=None):
+            super(TestInterpScalar, self).__init__()
+            self._x = x
+            self._left = left
+            self._right = right
+            self._period = period
+        
+        def hybrid_forward(self, F, xp, fp):
+            return F.np.interp(self._x, xp, fp, left=self._left, right=self._right, period=self._period)
+
+    xtypes = [np.int64, np.float32, np.float64]
+    dtypes = [np.int32, np.int64, np.float32, np.float64]
+    xshapes = [
+        (), (3,), (5,), (20,),
+        (2, 2), (4, 4), (8, 8),
+        (5, 5, 5), (8, 0, 8)
+    ]
+    dsizes = [10, 30]
+    periods = [None, 2*np.pi]
+    lefts = [None, -10, 0]
+    rights= [None, 20, 50]
+    flags = [True, False]
+    combinations = itertools.product(flags, flags, xshapes, dsizes, xtypes, dtypes, lefts, rights, periods)
+    for hybridize, x_scalar, xshape, dsize, xtype, dtype, left, right, period in combinations:
+        rtol = 1e-3
+        atol = 1e-5
+        if period is not None:
+            x = np.random.uniform(-np.pi, np.pi, size=xshape).astype(xtype)
+            xp = np.random.uniform(0, 2*np.pi, size=dsize)
+            fp = np.sin(xp)
+        else:
+            x = np.random.uniform(0, 100, size=xshape).astype(xtype)
+            xp = np.sort(np.random.choice(100, dsize, replace=False).astype(dtype))
+            fp = np.random.uniform(-50, 50, size=dsize).astype(dtype) 
+        np_x = x.asnumpy()
+        if x_scalar and xshape == ():
+            x = x.item()
+            np_x = x
+            test_interp = TestInterpScalar(x=x, left=left, right=right, period=period)
+        else: 
+            test_interp = TestInterp(left=left, right=right, period=period)
+        if hybridize:
+            test_interp.hybridize()
+        mx_out = test_interp(xp, fp) if (x_scalar and xshape == ()) else test_interp(x, xp, fp)
+        np_out = _np.interp(np_x, xp.asnumpy(), fp.asnumpy(), left=left, right=right, period=period)
+        assert mx_out.shape == np_out.shape
+        assert_almost_equal(mx_out.asnumpy(), np_out, atol=atol, rtol=rtol)
+
+        mx_out = np.interp(x, xp, fp, left=left, right=right, period=period)
+        np_out = _np.interp(np_x ,xp.asnumpy(), fp.asnumpy(), left=left, right=right, period=period)
+        assert_almost_equal(mx_out.asnumpy(), np_out, atol=atol, rtol=rtol)
+
+
 @with_seed()
 @use_np
 def test_np_bincount():
