diff --git a/doc/_static/dpnp-custom.css b/doc/_static/dpnp-custom.css
index 6c4495a213ba..c39eb9f307d3 100644
--- a/doc/_static/dpnp-custom.css
+++ b/doc/_static/dpnp-custom.css
@@ -57,14 +57,8 @@ dt.sig.sig-object .sig-param > .n * {
     font-weight: 700 !important;
 }
 
-/* Parameter/return descriptions: indented block on new line (via custom.js) */
-dl.field-list dd .param-desc {
-    display: block;
-    padding-left: 1.5em;
-}
-
-/* Parameter lists: no bullets, keep indentation */
-dl.field-list dd ul.simple {
-    list-style: none !important;
-    padding-left: 1.2em !important;
+/* numpydoc param/return types: italic like NumPy, not bold */
+.classifier {
+    font-style: italic;
+    font-weight: 400;
 }
diff --git a/doc/_static/dpnp-custom.js b/doc/_static/dpnp-custom.js
deleted file mode 100644
index b40ef0365036..000000000000
--- a/doc/_static/dpnp-custom.js
+++ /dev/null
@@ -1,67 +0,0 @@
-(function() {
-var separators = [ " – ", " -- " ];
-
-function findSeparator(container)
-{
-    var walker = document.createTreeWalker(container, NodeFilter.SHOW_TEXT);
-    var node;
-    while ((node = walker.nextNode())) {
-        for (var i = 0; i < separators.length; i++) {
-            var idx = node.nodeValue.indexOf(separators[i]);
-            if (idx !== -1)
-                return {node : node, offset : idx, sep : separators[i]};
-        }
-    }
-    return null;
-}
-
-// Splits <p> at the separator; wraps everything after it in .param-desc.
-function reformatP(p)
-{
-    var found = findSeparator(p);
-    if (!found)
-        return null;
-
-    var afterNode = found.node.splitText(found.offset);
-    afterNode.nodeValue = afterNode.nodeValue.slice(found.sep.length);
-
-    var range = document.createRange();
-    range.setStartBefore(afterNode);
-    range.setEndAfter(p.lastChild);
-
-    var desc = document.createElement("span");
-    desc.className = "param-desc";
-    desc.appendChild(range.extractContents());
-    p.appendChild(desc);
-    return desc;
-}
-
-// Browsers auto-close nested <p> tags, so multi-paragraph descriptions
-// arrive as sibling <p> elements inside <li>. Fold them into the same desc.
-function reformatEntry(container)
-{
-    var firstP = container.querySelector(":scope > p");
-    if (!firstP)
-        return;
-
-    var desc = reformatP(firstP);
-    if (!desc)
-        return;
-
-    var sibling;
-    while ((sibling = firstP.nextElementSibling) && sibling.tagName === "P") {
-        if (desc.textContent.trim())
-            desc.appendChild(document.createElement("br"));
-        while (sibling.firstChild)
-            desc.appendChild(sibling.firstChild);
-        sibling.remove();
-    }
-}
-
-document.querySelectorAll("dl.field-list dd ul.simple li")
-    .forEach(reformatEntry);
-document.querySelectorAll("dl.field-list dd").forEach(function(dd) {
-    if (!dd.querySelector("ul.simple"))
-        reformatEntry(dd);
-});
-}());
diff --git a/doc/conf.py b/doc/conf.py
index 0c31bc413caa..eaa6f7ce4070 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -8,8 +8,10 @@
 from datetime import datetime
 from urllib.parse import urljoin
 
+from jinja2.sandbox import SandboxedEnvironment
+from numpydoc.docscrape import NumpyDocString
+from numpydoc.docscrape_sphinx import SphinxDocString
 from sphinx.ext.autodoc import FunctionDocumenter
-from sphinx.ext.napoleon import NumpyDocstring, docstring
 
 from dpnp.dpnp_algo.dpnp_elementwise_common import (
     DPNPBinaryFunc,
@@ -64,7 +66,7 @@
     "sphinx.ext.viewcode",
     "sphinx.ext.githubpages",
     "sphinx.ext.intersphinx",
-    "sphinx.ext.napoleon",
+    "numpydoc",
     "sphinx.ext.autodoc",
     "sphinx.ext.autosummary",
     "sphinx_copybutton",
@@ -238,66 +240,64 @@ def _can_document_member(member, *args, **kwargs):
 # If true, `todo` and `todoList` produce output, else they produce nothing.
 todo_include_todos = True
 
-# Napoleon settings
-napoleon_use_ivar = True
-napoleon_include_special_with_doc = True
-napoleon_custom_sections = ["limitations"]
-
-
-# Napoleon extension can't properly render "Returns" section in case of
-# namedtuple as a return type. That patch proposes to extend the parse logic
-# which allows text in a header of "Returns" section.
-def _parse_returns_section_patched(self, section: str) -> list[str]:
-    fields = self._consume_returns_section()
-    multi = len(fields) > 1
-    use_rtype = False if multi else self._config.napoleon_use_rtype
-    lines: list[str] = []
-    header: list[str] = []
-    is_logged_header = False
-
-    for _name, _type, _desc in fields:
-        # self._consume_returns_section() stores the header block
-        # into `_type` argument, while `_name` has to be empty string and
-        # `_desc` has to be empty list of strings
-        if _name == "" and (not _desc or len(_desc) == 1 and _desc[0] == ""):
-            if not is_logged_header:
-                docstring.logger.info(
-                    "parse a header block of 'Returns' section",
-                    location=self._get_location(),
-                )
-                is_logged_header = True
-
-            # build a list with lines of the header block
-            header.extend([_type])
-            continue
+# Members come from autosummary; don't let numpydoc duplicate them
+numpydoc_show_class_members = False
+
+# Keep the dpnp-only "Limitations" section (numpydoc drops unknown sections):
+# register it and give it a slot in the template below
+NumpyDocString.sections.setdefault("Limitations", [])
+
+_NUMPYDOC_TEMPLATE = """\
+{{index}}
+{{summary}}
+{{extended_summary}}
+{{parameters}}
+{{attributes}}
+{{methods}}
+{{returns}}
+{{yields}}
+{{receives}}
+{{other_parameters}}
+{{raises}}
+{{warns}}
+{{warnings}}
+{{limitations}}
+{{see_also}}
+{{notes}}
+{{references}}
+{{examples}}
+"""
+
+_orig_load_config = SphinxDocString.load_config
+
+
+def _load_config_with_limitations(self, config):
+    _orig_load_config(self, config)
+    # Use our template with the "limitations" slot
+    self.template = SandboxedEnvironment().from_string(_NUMPYDOC_TEMPLATE)
+
+
+SphinxDocString.load_config = _load_config_with_limitations
+
+_orig_str = SphinxDocString.__str__
+
+
+def _str_with_limitations(self, indent=0, func_role="obj"):
+    # Wrap render() to fill the "limitations" slot (a rubric, like "Notes")
+    orig_render = self.template.render
+
+    def render(**ns):
+        ns["limitations"] = "\n".join(self._str_section("Limitations"))
+        return orig_render(**ns)
+
+    self.template.render = render
+    try:
+        return _orig_str(self, indent=indent, func_role=func_role)
+    finally:
+        self.template.render = orig_render
+
 
-        if use_rtype:
-            field = self._format_field(_name, "", _desc)
-        else:
-            field = self._format_field(_name, _type, _desc)
-
-        if multi:
-            if lines:
-                lines.extend(self._format_block("          * ", field))
-            else:
-                if header:
-                    # add the header block + the 1st parameter stored in `field`
-                    lines.extend([":returns:", ""])
-                    lines.extend(self._format_block(" " * 4, header))
-                    lines.extend(self._format_block("          * ", field))
-                else:
-                    lines.extend(self._format_block(":returns: * ", field))
-        else:
-            if any(field):  # only add :returns: if there's something to say
-                lines.extend(self._format_block(":returns: ", field))
-            if _type and use_rtype:
-                lines.extend([f":rtype: {_type}", ""])
-    if lines and lines[-1]:
-        lines.append("")
-    return lines
-
-
-NumpyDocstring._parse_returns_section = _parse_returns_section_patched
+SphinxDocString.__str__ = _str_with_limitations
 
 
 # TODO: Remove once dpnp.tensor docs are generated in dpnp
diff --git a/dpnp/dpnp_array.py b/dpnp/dpnp_array.py
index b225fb2c7329..5a93c8018a82 100644
--- a/dpnp/dpnp_array.py
+++ b/dpnp/dpnp_array.py
@@ -1986,8 +1986,8 @@ def sort(
         :obj:`dpnp.searchsorted` : Find elements in a sorted array.
         :obj:`dpnp.partition` : Partial sort.
 
-        Note
-        ----
+        Notes
+        -----
         `axis` in :obj:`dpnp.sort` could be integer or ``None``. If ``None``,
         the array is flattened before sorting. However, `axis` in
         :obj:`dpnp.ndarray.sort` can only be integer since it sorts an array
diff --git a/dpnp/dpnp_iface_linearalgebra.py b/dpnp/dpnp_iface_linearalgebra.py
index b8c01bdc854b..1d0843847631 100644
--- a/dpnp/dpnp_iface_linearalgebra.py
+++ b/dpnp/dpnp_iface_linearalgebra.py
@@ -257,7 +257,7 @@ def einsum(
         The calculation based on the Einstein summation convention.
 
     See Also
-    -------
+    --------
     :obj:`dpnp.einsum_path` : Evaluates the lowest cost contraction order
                               for an einsum expression.
     :obj:`dpnp.dot` : Returns the dot product of two arrays.
diff --git a/dpnp/dpnp_iface_manipulation.py b/dpnp/dpnp_iface_manipulation.py
index 8d3050da4e0e..dd53033db279 100644
--- a/dpnp/dpnp_iface_manipulation.py
+++ b/dpnp/dpnp_iface_manipulation.py
@@ -749,7 +749,7 @@ def asarray_chkfinite(
         already an ndarray.
 
     Raises
-    -------
+    ------
     ValueError
         Raises ``ValueError`` if `a` contains NaN (Not a Number) or
         Inf (Infinity).
@@ -4537,8 +4537,6 @@ def unique_all(x, /):
 
     Returns
     -------
-    A namedtuple with the following attributes:
-
     values : dpnp.ndarray
         The unique elements of an input array.
     indices : dpnp.ndarray
@@ -4591,8 +4589,6 @@ def unique_counts(x, /):
 
     Returns
     -------
-    A namedtuple with the following attributes:
-
     values : dpnp.ndarray
         The unique elements of an input array.
     counts : dpnp.ndarray
@@ -4637,8 +4633,6 @@ def unique_inverse(x, /):
 
     Returns
     -------
-    A namedtuple with the following attributes:
-
     values : dpnp.ndarray
         The unique elements of an input array.
     inverse_indices : dpnp.ndarray
diff --git a/dpnp/dpnp_iface_mathematical.py b/dpnp/dpnp_iface_mathematical.py
index c66d827faff6..7b7ffba9dbfa 100644
--- a/dpnp/dpnp_iface_mathematical.py
+++ b/dpnp/dpnp_iface_mathematical.py
@@ -2329,7 +2329,7 @@ def ediff1d(ary, to_end=None, to_begin=None):
     returned array is determined by the Type Promotion Rules.
 
 Limitations
-----------
+-----------
 Parameters `where` and `subok` are supported with their default values.
 Keyword argument `kwargs` is currently unsupported.
 Otherwise ``NotImplementedError`` exception will be raised.
@@ -3929,8 +3929,8 @@ def _check_nan_inf(val, val_dt):
 :obj:`dpnp.negative` : Return the numerical negative of each element of `x`.
 :obj:`dpnp.copysign` : Change the sign of `x1` to that of `x2`, element-wise.
 
-Note
-----
+Notes
+-----
 Equivalent to `x.copy()`, but only defined for types that support arithmetic.
 
 Examples
diff --git a/dpnp/dpnp_iface_nanfunctions.py b/dpnp/dpnp_iface_nanfunctions.py
index 10fffb342305..b528316b485f 100644
--- a/dpnp/dpnp_iface_nanfunctions.py
+++ b/dpnp/dpnp_iface_nanfunctions.py
@@ -129,8 +129,6 @@ def nanargmax(a, axis=None, out=None, *, keepdims=False):
     the user is recommended to filter NaNs themselves and use `dpnp.argmax`
     on the filtered array.
 
-    Warnings
-    --------
     The results cannot be trusted if a slice contains only NaNs
     and -Infs.
 
@@ -213,8 +211,6 @@ def nanargmin(a, axis=None, out=None, *, keepdims=False):
     the user is recommended to filter NaNs themselves and use `dpnp.argmax`
     on the filtered array.
 
-    Warnings
-    --------
     The results cannot be trusted if a slice contains only NaNs
     and -Infs.
 
diff --git a/dpnp/dpnp_iface_trigonometric.py b/dpnp/dpnp_iface_trigonometric.py
index ee3462ab610e..b82f090de476 100644
--- a/dpnp/dpnp_iface_trigonometric.py
+++ b/dpnp/dpnp_iface_trigonometric.py
@@ -907,8 +907,8 @@ def cumlogsumexp(
     :obj:`dpnp.logsumexp` : Logarithm of the sum of elements of the inputs,
                             element-wise.
 
-    Note
-    ----
+    Notes
+    -----
     This function is equivalent of `numpy.logaddexp.accumulate`.
 
     Examples
@@ -1889,8 +1889,8 @@ def logsumexp(x, /, *, axis=None, dtype=None, keepdims=False, out=None):
     :obj:`dpnp.cumlogsumexp` : Cumulative the natural logarithm of the sum of
         elements in the input array.
 
-    Note
-    ----
+    Notes
+    -----
     This function is equivalent of `numpy.logaddexp.reduce`.
 
     Examples
@@ -2171,8 +2171,8 @@ def reduce_hypot(x, /, *, axis=None, dtype=None, keepdims=False, out=None):
     --------
     :obj:`dpnp.hypot` : Calculates :math:`\sqrt{x1^2 + x2^2}`, element-wise.
 
-    Note
-    ----
+    Notes
+    -----
     This function is equivalent of `numpy.hypot.reduce`.
 
     Examples
diff --git a/dpnp/fft/dpnp_iface_fft.py b/dpnp/fft/dpnp_iface_fft.py
index 90e1a112bdaf..0a5c3cc4351d 100644
--- a/dpnp/fft/dpnp_iface_fft.py
+++ b/dpnp/fft/dpnp_iface_fft.py
@@ -559,7 +559,7 @@ def hfft(a, n=None, axis=-1, norm=None, out=None):
     --------
     :obj:`dpnp.fft` : For definition of the DFT and conventions used.
     :obj:`dpnp.fft.rfft` : The one-dimensional FFT of real input.
-    :obj:`dpnp.fft.ihfft` :The inverse of :obj:`dpnp.fft.hfft`.
+    :obj:`dpnp.fft.ihfft` : The inverse of :obj:`dpnp.fft.hfft`.
 
 
     Notes
@@ -1101,7 +1101,7 @@ def irfft(a, n=None, axis=-1, norm=None, out=None):
     :obj:`dpnp.fft.rfft` : The one-dimensional FFT of real input, of which
                         :obj:`dpnp.fft.irfft` is inverse.
     :obj:`dpnp.fft.fft` : The one-dimensional FFT of general (complex) input.
-    :obj:`dpnp.fft.irfft2` :The inverse of the two-dimensional FFT of
+    :obj:`dpnp.fft.irfft2` : The inverse of the two-dimensional FFT of
                         real input.
     :obj:`dpnp.fft.irfftn` : The inverse of the *N*-dimensional FFT of
                         real input.
diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py
index 76910692ea0c..74fe3b65febc 100644
--- a/dpnp/linalg/dpnp_iface_linalg.py
+++ b/dpnp/linalg/dpnp_iface_linalg.py
@@ -468,8 +468,6 @@ def eig(a):
 
     Returns
     -------
-    A namedtuple with the following attributes:
-
     eigenvalues : (..., M) dpnp.ndarray
         The eigenvalues, each repeated according to its multiplicity.
         The eigenvalues are not necessarily ordered. The resulting array is
@@ -481,8 +479,8 @@ def eig(a):
         ``eigenvectors[:,i]`` is the eigenvector corresponding to the
         eigenvalue ``eigenvalues[i]``.
 
-    Note
-    ----
+    Notes
+    -----
     Since there is no proper OneMKL LAPACK function, DPNP will calculate
     through a fallback on NumPy call.
 
@@ -584,8 +582,6 @@ def eigh(a, UPLO="L"):
 
     Returns
     -------
-    A namedtuple with the following attributes:
-
     eigenvalues : (..., M) dpnp.ndarray
         The eigenvalues in ascending order, each repeated according to its
         multiplicity.
@@ -645,8 +641,8 @@ def eigvals(a):
         They are not necessarily ordered, nor are they necessarily
         real for real matrices.
 
-    Note
-    ----
+    Notes
+    -----
     Since there is no proper OneMKL LAPACK function, DPNP will calculate
     through a fallback on NumPy call.
 
@@ -1591,9 +1587,6 @@ def qr(a, mode="reduced"):
 
     Returns
     -------
-    When mode is "reduced" or "complete", the result will be a namedtuple with
-    the attributes `Q` and `R`:
-
     Q : dpnp.ndarray of float or complex, optional
         A matrix with orthonormal columns.
         When mode is ``"complete"`` the result is an orthogonal/unitary matrix
@@ -1609,6 +1602,12 @@ def qr(a, mode="reduced"):
         along with `R`. The `tau` array contains scaling factors for the
         reflectors.
 
+    Notes
+    -----
+    When `mode` is ``"reduced"`` or ``"complete"``, the result is a namedtuple
+    with the attributes `Q` and `R`. When `mode` is ``"raw"``, it returns
+    ``(h, tau)``.
+
     Examples
     --------
     >>> import dpnp as np
@@ -1749,9 +1748,6 @@ def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
 
     Returns
     -------
-    When `compute_uv` is ``True``, the result is a namedtuple with the
-    following attribute names:
-
     U : { (…, M, M), (…, M, K) } dpnp.ndarray
         Unitary matrix, where M is the number of rows of the input array `a`.
         The shape of the matrix `U` depends on the value of `full_matrices`.
@@ -1769,6 +1765,12 @@ def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
         If `full_matrices` is ``False``, `Vh` has the shape (…, K, N).
         If `compute_uv` is ``False``, neither `U` or `Vh` are computed.
 
+    Notes
+    -----
+    When `compute_uv` is ``True``, the result is a namedtuple with the
+    attribute names `U`, `S`, and `Vh`. When `compute_uv` is ``False``, only
+    the singular values `S` are returned.
+
     Examples
     --------
     >>> import dpnp as np
@@ -1895,8 +1897,6 @@ def slogdet(a):
 
     Returns
     -------
-    A namedtuple with the following attributes:
-
     sign : (...) dpnp.ndarray
         A number representing the sign of the determinant. For a real matrix,
         this is 1, 0, or -1. For a complex matrix, this is a complex number
diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py
index ac1c87a61b5b..5639de12ae5a 100644
--- a/dpnp/linalg/dpnp_utils_linalg.py
+++ b/dpnp/linalg/dpnp_utils_linalg.py
@@ -736,7 +736,7 @@ def _calculate_determinant_sign(ipiv, diag, res_type, n):
     values.
 
     Parameters
-    -----------
+    ----------
     ipiv : {dpnp.ndarray, usm_ndarray}
         The pivot indices from LU decomposition.
     diag : {dpnp.ndarray, usm_ndarray}
diff --git a/dpnp/scipy/linalg/_decomp_lu.py b/dpnp/scipy/linalg/_decomp_lu.py
index f96d56b0e423..5a57e962937e 100644
--- a/dpnp/scipy/linalg/_decomp_lu.py
+++ b/dpnp/scipy/linalg/_decomp_lu.py
@@ -95,9 +95,6 @@ def lu(
 
     Returns
     -------
-    The tuple ``(p, l, u)`` is returned if ``permute_l`` is ``False``
-    (default), else the tuple ``(pl, u)`` is returned, where:
-
     p : (..., M, M) dpnp.ndarray or (..., M) dpnp.ndarray
         Permutation matrix or permutation indices.
         If `p_indices` is ``False`` (default), a permutation matrix.
@@ -116,6 +113,9 @@ def lu(
 
     Notes
     -----
+    The tuple ``(p, l, u)`` is returned if `permute_l` is ``False`` (default),
+    else the tuple ``(pl, u)`` is returned.
+
     Permutation matrices are costly since they are nothing but row reorder of
     ``L`` and hence indices are strongly recommended to be used instead if the
     permutation is required. The relation in the 2D case then becomes simply
diff --git a/environments/building_docs.yml b/environments/building_docs.yml
index 6afaf65c9514..79407e4fdbcd 100644
--- a/environments/building_docs.yml
+++ b/environments/building_docs.yml
@@ -6,5 +6,6 @@ dependencies:
   - cupy
   - sphinx
   - furo
+  - numpydoc
   - pip:
     - -r base_build_docs.txt
diff --git a/pyproject.toml b/pyproject.toml
index 0cf3e9721872..a9519d1744aa 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -89,6 +89,7 @@ docs = [
   "Cython",
   "cupy",
   "furo",
+  "numpydoc",
   "sphinx",
   "sphinx-copybutton",
   "sphinx-design",