diff --git a/.github/workflows/array-api-skips.txt b/.github/workflows/array-api-skips.txt
index 9589c591f5cb..1da2bc1ee5d7 100644
--- a/.github/workflows/array-api-skips.txt
+++ b/.github/workflows/array-api-skips.txt
@@ -17,13 +17,7 @@ array_api_tests/test_signatures.py::test_func_signature[unique_counts]
 array_api_tests/test_signatures.py::test_func_signature[unique_inverse]
 array_api_tests/test_signatures.py::test_func_signature[unique_values]
 
-# do not return a namedtuple
-array_api_tests/test_linalg.py::test_eigh
-array_api_tests/test_linalg.py::test_slogdet
-array_api_tests/test_linalg.py::test_svd
-
 # hypothesis found failures
-array_api_tests/test_linalg.py::test_qr
 array_api_tests/test_operators_and_elementwise_functions.py::test_clip
 
 # unexpected result is returned - unmute when dpctl-1986 is resolved
diff --git a/doc/conf.py b/doc/conf.py
index c2a475e09da9..03cd711c30f4 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -9,6 +9,7 @@
 from datetime import datetime
 
 from sphinx.ext.autodoc import FunctionDocumenter
+from sphinx.ext.napoleon import NumpyDocstring, docstring
 
 from dpnp.dpnp_algo.dpnp_elementwise_common import DPNPBinaryFunc, DPNPUnaryFunc
 
@@ -231,3 +232,59 @@ def _can_document_member(member, *args, **kwargs):
 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
+
+        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
diff --git a/dpnp/dpnp_array.py b/dpnp/dpnp_array.py
index 4e3625efa6c2..af6d20d627f5 100644
--- a/dpnp/dpnp_array.py
+++ b/dpnp/dpnp_array.py
@@ -760,12 +760,14 @@ def argsort(
 
     def asnumpy(self):
         """
-        Copy content of the array into :class:`numpy.ndarray` instance of the same shape and data type.
+        Copy content of the array into :class:`numpy.ndarray` instance of
+        the same shape and data type.
 
         Returns
         -------
-        numpy.ndarray
-            An instance of :class:`numpy.ndarray` populated with the array content.
+        out : numpy.ndarray
+            An instance of :class:`numpy.ndarray` populated with the array
+            content.
 
         """
 
diff --git a/dpnp/dpnp_iface_manipulation.py b/dpnp/dpnp_iface_manipulation.py
index da5720ff2152..0c796c61961c 100644
--- a/dpnp/dpnp_iface_manipulation.py
+++ b/dpnp/dpnp_iface_manipulation.py
@@ -1124,7 +1124,7 @@ def broadcast_shapes(*args):
 
     Returns
     -------
-    tuple
+    out : tuple
         Broadcasted shape.
 
     See Also
diff --git a/dpnp/dpnp_iface_statistics.py b/dpnp/dpnp_iface_statistics.py
index 7f320a166cab..4c93f634b260 100644
--- a/dpnp/dpnp_iface_statistics.py
+++ b/dpnp/dpnp_iface_statistics.py
@@ -843,7 +843,7 @@ def max(a, axis=None, out=None, keepdims=False, initial=None, where=True):
         dimension ``a.ndim - len(axis)``.
 
     Limitations
-    -----------.
+    -----------
     Parameters `where`, and `initial` are only supported with their default
     values. Otherwise ``NotImplementedError`` exception will be raised.
 
diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py
index 377bbcdb37c2..364fb4de2c0a 100644
--- a/dpnp/linalg/dpnp_iface_linalg.py
+++ b/dpnp/linalg/dpnp_iface_linalg.py
@@ -39,6 +39,8 @@
 # pylint: disable=invalid-name
 # pylint: disable=no-member
 
+from typing import NamedTuple
+
 import numpy
 from dpctl.tensor._numpy_helper import normalize_axis_tuple
 
@@ -100,6 +102,12 @@
 ]
 
 
+# pylint:disable=missing-class-docstring
+class EigResult(NamedTuple):
+    eigenvalues: dpnp.ndarray
+    eigenvectors: dpnp.ndarray
+
+
 def cholesky(a, /, *, upper=False):
     """
     Cholesky decomposition.
@@ -123,6 +131,7 @@ def cholesky(a, /, *, upper=False):
     upper : {bool}, optional
         If ``True``, the result must be the upper-triangular Cholesky factor.
         If ``False``, the result must be the lower-triangular Cholesky factor.
+
         Default: ``False``.
 
     Returns
@@ -174,6 +183,7 @@ def cond(x, p=None):
         Order of the norm used in the condition number computation.
         ``inf`` means the `dpnp.inf` object, and the Frobenius norm is
         the root-of-sum-of-squares norm.
+
         Default: ``None``.
 
     Returns
@@ -244,6 +254,7 @@ def cross(x1, x2, /, *, axis=-1):
     axis : int, optional
         The axis (dimension) of `x1` and `x2` containing the vectors for
         which to compute the cross-product.
+
         Default: ``-1``.
 
     Returns
@@ -451,17 +462,18 @@ 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 will be of complex type, unless the imaginary part is
-        zero in which case it will be cast to a real type. When `a`
-        is real the resulting eigenvalues will be real (0 imaginary
-        part) or occur in conjugate pairs
+        The eigenvalues are not necessarily ordered. The resulting array will
+        be of complex type, unless the imaginary part is zero in which case it
+        will be cast to a real type. When `a` is real the resulting eigenvalues
+        will be real (zero imaginary part) or occur in conjugate pairs.
     eigenvectors : (..., M, M) dpnp.ndarray
-        The normalized (unit "length") eigenvectors, such that the
-        column ``v[:,i]`` is the eigenvector corresponding to the
-        eigenvalue ``w[i]``.
+        The normalized (unit "length") eigenvectors, such that the column
+        ``eigenvectors[:,i]`` is the eigenvector corresponding to the
+        eigenvalue ``eigenvalues[i]``.
 
     Note
     ----
@@ -532,7 +544,7 @@ def eig(a):
     # Since geev function from OneMKL LAPACK is not implemented yet,
     # use NumPy for this calculation.
     w_np, v_np = numpy.linalg.eig(dpnp.asnumpy(a))
-    return (
+    return EigResult(
         dpnp.array(w_np, sycl_queue=a_sycl_queue, usm_type=a_usm_type),
         dpnp.array(v_np, sycl_queue=a_sycl_queue, usm_type=a_usm_type),
     )
@@ -561,16 +573,19 @@ def eigh(a, UPLO="L"):
         considered to preserve the Hermite matrix property.
         It therefore follows that the imaginary part of the diagonal
         will always be treated as zero.
+
         Default: ``"L"``.
 
     Returns
     -------
-    w : (..., M) dpnp.ndarray
-        The eigenvalues in ascending order, each repeated according to
-        its multiplicity.
-    v : (..., M, M) dpnp.ndarray
-        The column ``v[:, i]`` is the normalized eigenvector corresponding
-        to the eigenvalue ``w[i]``.
+    A namedtuple with the following attributes:
+
+    eigenvalues : (..., M) dpnp.ndarray
+        The eigenvalues in ascending order, each repeated according to its
+        multiplicity.
+    eigenvectors : (..., M, M) dpnp.ndarray
+        The column ``eigenvectors[:, i]`` is the normalized eigenvector
+        corresponding to the eigenvalue ``eigenvalues[i]``.
 
     See Also
     --------
@@ -644,7 +659,7 @@ def eigvals(a):
     Illustration, using the fact that the eigenvalues of a diagonal matrix
     are its diagonal elements, that multiplying a matrix on the left
     by an orthogonal matrix, `Q`, and on the right by `Q.T` (the transpose
-    of `Q`), preserves the eigenvalues of the "middle" matrix.  In other words,
+    of `Q`), preserves the eigenvalues of the "middle" matrix. In other words,
     if `Q` is orthogonal, then ``Q * A * Q.T`` has the same eigenvalues as
     ``A``:
 
@@ -698,6 +713,7 @@ def eigvalsh(a, UPLO="L"):
         considered to preserve the Hermite matrix property.
         It therefore follows that the imaginary part of the diagonal
         will always be treated as zero.
+
         Default: ``"L"``.
 
     Returns
@@ -809,6 +825,7 @@ def lstsq(a, b, rcond=None):
         of `a`.
         The default uses the machine precision times ``max(M, N)``. Passing
         ``-1`` will use machine precision.
+
         Default: ``None``.
 
     Returns
@@ -839,7 +856,7 @@ def lstsq(a, b, rcond=None):
     gradient of roughly 1 and cut the y-axis at, more or less, -1.
 
     We can rewrite the line equation as ``y = Ap``, where ``A = [[x 1]]``
-    and ``p = [[m], [c]]``.  Now use `lstsq` to solve for `p`:
+    and ``p = [[m], [c]]``. Now use `lstsq` to solve for `p`:
 
     >>> A = np.vstack([x, np.ones(len(x))]).T
     >>> A
@@ -949,10 +966,12 @@ def matrix_norm(x, /, *, keepdims=False, ord="fro"):
         If this is set to ``True``, the axes which are normed over are left in
         the result as dimensions with size one. With this option the result
         will broadcast correctly against the original `x`.
+
         Default: ``False``.
     ord : {None, 1, -1, 2, -2, dpnp.inf, -dpnp.inf, 'fro', 'nuc'}, optional
         The order of the norm. For details see the table under ``Notes``
         section in :obj:`dpnp.linalg.norm`.
+
         Default: ``"fro"``.
 
     Returns
@@ -1079,16 +1098,19 @@ def matrix_rank(A, tol=None, hermitian=False, *, rtol=None):
         `rtol` can be set at a time. If none of them are provided, defaults
         to ``S.max() * max(M, N) * eps`` where `S` is an array with singular
         values for `A`, and `eps` is the epsilon value for datatype of `S`.
+
         Default: ``None``.
     hermitian : bool, optional
         If ``True``, `A` is assumed to be Hermitian (symmetric if real-valued),
         enabling a more efficient method for finding singular values.
+
         Default: ``False``.
     rtol : (...) {None, float, dpnp.ndarray, usm_ndarray}, optional
         Parameter for the relative tolerance component. Only `tol` or `rtol`
         can be set at a time. If none of them are provided, defaults to
         ``max(M, N) * eps`` where `eps` is the epsilon value for datatype
         of `S` (an array with singular values for `A`).
+
         Default: ``None``.
 
     Returns
@@ -1191,6 +1213,7 @@ def multi_dot(arrays, *, out=None):
         C-contiguous, and its dtype must be the dtype that would be returned
         for `dot(a, b)`. If these conditions are not met, an exception is
         raised, instead of attempting to be flexible.
+
         Default: ``None``.
 
     Returns
@@ -1252,11 +1275,12 @@ def norm(x, ord=None, axis=None, keepdims=False):
     Parameters
     ----------
     x : {dpnp.ndarray, usm_ndarray}
-        Input array.  If `axis` is ``None``, `x` must be 1-D or 2-D, unless
+        Input array. If `axis` is ``None``, `x` must be 1-D or 2-D, unless
         `ord` is ``None``. If both `axis` and `ord` are ``None``, the 2-norm
         of ``x.ravel`` will be returned.
     ord : {int, float, inf, -inf, "fro", "nuc"}, optional
         Norm type. inf means dpnp's `inf` object.
+
         Default: ``None``.
     axis : {None, int, 2-tuple of ints}, optional
         If `axis` is an integer, it specifies the axis of `x` along which to
@@ -1264,11 +1288,13 @@ def norm(x, ord=None, axis=None, keepdims=False):
         axes that hold 2-D matrices, and the matrix norms of these matrices
         are computed. If `axis` is ``None`` then either a vector norm (when
         `x` is 1-D) or a matrix norm (when `x` is 2-D) is returned.
+
         Default: ``None``.
     keepdims : bool, optional
         If this is set to ``True``, the axes which are normed over are left in
         the result as dimensions with size one. With this option the result
         will broadcast correctly against the original `x`.
+
         Default: ``False``.
 
     Returns
@@ -1486,15 +1512,18 @@ def pinv(a, rcond=None, hermitian=False, *, rtol=None):
         are set to zero. Broadcasts against the stack of matrices.
         Only `rcond` or `rtol` can be set at a time. If none of them are
         provided, defaults to ``max(M, N) * dpnp.finfo(a.dtype).eps``.
+
         Default: ``None``.
     hermitian : bool, optional
         If ``True``, a is assumed to be Hermitian (symmetric if real-valued),
         enabling a more efficient method for finding singular values.
+
         Default: ``False``.
     rtol : (...) {None, float, dpnp.ndarray, usm_ndarray}, optional
         Same as `rcond`, but it's an Array API compatible parameter name.
         Only `rcond` or `rtol` can be set at a time. If none of them are
         provided, defaults to ``max(M, N) * dpnp.finfo(a.dtype).eps``.
+
         Default: ``None``.
 
     Returns
@@ -1557,20 +1586,22 @@ 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
+    the attributes `Q` and `R`:
+
+    Q : dpnp.ndarray of float or complex, optional
         A matrix with orthonormal columns.
-        When mode = "complete" the result is an orthogonal/unitary matrix
-        depending on whether or not a is real/complex.
-        The determinant may be either +/- 1 in that case.
-        In case the number of dimensions in the input array is greater
-        than 2 then a stack of the matrices with above properties is returned.
-    R : dpnp.ndarray
-        The upper-triangular matrix or a stack of upper-triangular matrices
-        if the number of dimensions in the input array is greater than 2.
-    (h, tau) : tuple of dpnp.ndarray
-        The `h` array contains the Householder reflectors that generate Q along
-        with R. The `tau` array contains scaling factors for the reflectors.
+        When mode is ``"complete"`` the result is an orthogonal/unitary matrix
+        depending on whether or not `a` is real/complex. The determinant may be
+        either ``+/- 1`` in that case. In case the number of dimensions in the
+        input array is greater than 2 then a stack of the matrices with above
+        properties is returned.
+    R : dpnp.ndarray of float or complex, optional
+        The upper-triangular matrix or a stack of upper-triangular matrices if
+        the number of dimensions in the input array is greater than 2.
+    (h, tau) : tuple of dpnp.ndarray of float or complex, optional
+        The array `h` contains the Householder reflectors that generate `Q`
+        along with `R`. The `tau` array contains scaling factors for the
+        reflectors.
 
     Examples
     --------
@@ -1698,33 +1729,39 @@ def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
     full_matrices : {bool}, optional
         If ``True``, it returns `u` and `Vh` with full-sized matrices.
         If ``False``, the matrices are reduced in size.
+
         Default: ``True``.
     compute_uv : {bool}, optional
         If ``False``, it only returns singular values.
+
         Default: ``True``.
     hermitian : {bool}, optional
         If True, a is assumed to be Hermitian (symmetric if real-valued),
         enabling a more efficient method for finding singular values.
+
         Default: ``False``.
 
     Returns
     -------
-    u : { (…, M, M), (…, M, K) } dpnp.ndarray
+    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`.
-        If `full_matrices` is ``True``, `u` has the shape (…, M, M).
-        If `full_matrices` is ``False``, `u` has the shape (…, M, K), where
-        K = min(M, N), and N is the number of columns of the input array `a`.
-        If `compute_uv` is ``False``, neither `u` or `Vh` are computed.
-    s : (…, K) dpnp.ndarray
+        The shape of the matrix `U` depends on the value of `full_matrices`.
+        If `full_matrices` is ``True``, `U` has the shape (…, M, M). If
+        `full_matrices` is ``False``, `U` has the shape (…, M, K), where
+        ``K = min(M, N)``, and N is the number of columns of the input array
+        `a`. If `compute_uv` is ``False``, neither `U` or `Vh` are computed.
+    S : (…, K) dpnp.ndarray
         Vector containing the singular values of `a`, sorted in descending
-        order. The length of `s` is min(M, N).
+        order. The length of `S` is min(M, N).
     Vh : { (…, N, N), (…, K, N) } dpnp.ndarray
         Unitary matrix, where N is the number of columns of the input array `a`.
         The shape of the matrix `Vh` depends on the value of `full_matrices`.
         If `full_matrices` is ``True``, `Vh` has the shape (…, N, N).
         If `full_matrices` is ``False``, `Vh` has the shape (…, K, N).
-        If `compute_uv` is ``False``, neither `u` or `Vh` are computed.
+        If `compute_uv` is ``False``, neither `U` or `Vh` are computed.
 
     Examples
     --------
@@ -1852,6 +1889,8 @@ 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
@@ -1924,6 +1963,8 @@ def tensordot(a, b, /, *, axes=2):
           applying to `a`, second to `b`. Both elements array_like must be of
           the same length.
 
+          Default: ``2``.
+
     Returns
     -------
     out : dpnp.ndarray
@@ -2009,6 +2050,7 @@ def tensorinv(a, ind=2):
     ind : int, optional
         Number of first indices that are involved in the inverse sum.
         Must be a positive integer.
+
         Default: ``2``.
 
     Returns
@@ -2074,6 +2116,7 @@ def tensorsolve(a, b, axes=None):
     axes : {None, tuple of ints}, optional
         Axes in `a` to reorder to the right, before inversion.
         If ``None`` , no reordering is done.
+
         Default: ``None``.
 
     Returns
@@ -2158,6 +2201,7 @@ def trace(x, /, *, offset=0, dtype=None):
         `a` is of integer type of precision less than the default integer
         precision, then the default integer precision is used. Otherwise, the
         precision is the same as that of `a`.
+
         Default: ``None``.
 
     Returns
@@ -2243,6 +2287,7 @@ def vecdot(x1, x2, /, *, axis=-1):
         Second input array.
     axis : int, optional
         Axis over which to compute the dot product.
+
         Default: ``-1``.
 
     Returns
@@ -2289,15 +2334,18 @@ def vector_norm(x, /, *, axis=None, keepdims=False, ord=2):
         (dimensions) along which to compute batched vector norms. If ``None``,
         the vector norm must be computed over all array values (i.e.,
         equivalent to computing the vector norm of a flattened array).
+
         Default: ``None``.
     keepdims : bool, optional
         If this is set to ``True``, the axes which are normed over are left in
         the result as dimensions with size one. With this option the result
         will broadcast correctly against the original `x`.
+
         Default: ``False``.
     ord : {int, float, inf, -inf, 'fro', 'nuc'}, optional
         The order of the norm. For details see the table under ``Notes``
         section in :obj:`dpnp.linalg.norm`.
+
         Default: ``2``.
 
     Returns
diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py
index 97ff26a70432..f105cff5a6c2 100644
--- a/dpnp/linalg/dpnp_utils_linalg.py
+++ b/dpnp/linalg/dpnp_utils_linalg.py
@@ -38,6 +38,8 @@
 # pylint: disable=protected-access
 # pylint: disable=useless-import-alias
 
+from typing import NamedTuple
+
 import dpctl.tensor._tensor_impl as ti
 import dpctl.utils as dpu
 import numpy
@@ -70,6 +72,29 @@
     "dpnp_svd",
 ]
 
+
+# pylint:disable=missing-class-docstring
+class EighResult(NamedTuple):
+    eigenvalues: dpnp.ndarray
+    eigenvectors: dpnp.ndarray
+
+
+class QRResult(NamedTuple):
+    Q: dpnp.ndarray
+    R: dpnp.ndarray
+
+
+class SlogdetResult(NamedTuple):
+    sign: dpnp.ndarray
+    logabsdet: dpnp.ndarray
+
+
+class SVDResult(NamedTuple):
+    U: dpnp.ndarray
+    S: dpnp.ndarray
+    Vh: dpnp.ndarray
+
+
 _jobz = {"N": 0, "V": 1}
 _upper_lower = {"U": 0, "L": 1}
 
@@ -162,7 +187,7 @@ def _batched_eigh(a, UPLO, eigen_mode, w_type, v_type):
         # Convert to contiguous to align with NumPy
         if a_orig_order == "C":
             v = dpnp.ascontiguousarray(v)
-        return w, v
+        return EighResult(w, v)
     return w
 
 
@@ -476,7 +501,7 @@ def _batched_qr(a, mode="reduced"):
 
     r = _triu_inplace(r)
 
-    return (
+    return QRResult(
         q.reshape(batch_shape + q.shape[-2:]),
         r.reshape(batch_shape + r.shape[-2:]),
     )
@@ -632,7 +657,7 @@ def _batched_svd(
         u = dpnp.ascontiguousarray(u)
         vt = dpnp.ascontiguousarray(vt)
         # Swap `u` and `vt` for transposed input to restore correct order
-        return (vt, s, u) if trans_flag else (u, s, vt)
+        return SVDResult(vt, s, u) if trans_flag else SVDResult(u, s, vt)
     return s
 
 
@@ -819,9 +844,9 @@ def _hermitian_svd(a, compute_uv):
     # but dpnp.linalg.eigh returns s sorted ascending so we re-order
     # the eigenvalues and related arrays to have the correct order
     if compute_uv:
-        s, u = dpnp.linalg.eigh(a)
+        s, u = dpnp_eigh(a, eigen_mode="V")
         sgn = dpnp.sign(s)
-        s = dpnp.absolute(s)
+        s = dpnp.abs(s, out=s)
         sidx = dpnp.argsort(s)[..., ::-1]
         # Rearrange the signs according to sorted indices
         sgn = dpnp.take_along_axis(sgn, sidx, axis=-1)
@@ -832,11 +857,10 @@ def _hermitian_svd(a, compute_uv):
         # Singular values are unsigned, move the sign into v
         # Compute V^T adjusting for the sign and conjugating
         vt = dpnp.transpose(u * sgn[..., None, :]).conjugate()
-        return u, s, vt
+        return SVDResult(u, s, vt)
 
-    # TODO: use dpnp.linalg.eighvals when it is updated
-    s, _ = dpnp.linalg.eigh(a)
-    s = dpnp.abs(s)
+    s = dpnp_eigh(a, eigen_mode="N")
+    s = dpnp.abs(s, out=s)
     return dpnp.sort(s)[..., ::-1]
 
 
@@ -1423,7 +1447,7 @@ def _zero_batched_qr(a, mode, m, n, k, res_type):
     batch_shape = a.shape[:-2]
 
     if mode == "reduced":
-        return (
+        return QRResult(
             dpnp.empty_like(
                 a,
                 shape=batch_shape + (m, k),
@@ -1443,7 +1467,7 @@ def _zero_batched_qr(a, mode, m, n, k, res_type):
             usm_type=a_usm_type,
             sycl_queue=a_sycl_queue,
         )
-        return (
+        return QRResult(
             q,
             dpnp.empty_like(
                 a,
@@ -1530,7 +1554,7 @@ def _zero_batched_svd(
             usm_type=usm_type,
             sycl_queue=exec_q,
         )
-        return u, s, vt
+        return SVDResult(u, s, vt)
     return s
 
 
@@ -1548,22 +1572,28 @@ def _zero_k_qr(a, mode, m, n, res_type):
     m, n = a.shape
 
     if mode == "reduced":
-        return dpnp.empty_like(
-            a,
-            shape=(m, 0),
-            dtype=res_type,
-        ), dpnp.empty_like(
-            a,
-            shape=(0, n),
-            dtype=res_type,
+        return QRResult(
+            dpnp.empty_like(
+                a,
+                shape=(m, 0),
+                dtype=res_type,
+            ),
+            dpnp.empty_like(
+                a,
+                shape=(0, n),
+                dtype=res_type,
+            ),
         )
     if mode == "complete":
-        return dpnp.identity(
-            m, dtype=res_type, sycl_queue=a_sycl_queue, usm_type=a_usm_type
-        ), dpnp.empty_like(
-            a,
-            shape=(m, n),
-            dtype=res_type,
+        return QRResult(
+            dpnp.identity(
+                m, dtype=res_type, sycl_queue=a_sycl_queue, usm_type=a_usm_type
+            ),
+            dpnp.empty_like(
+                a,
+                shape=(m, n),
+                dtype=res_type,
+            ),
         )
     if mode == "r":
         return dpnp.empty_like(
@@ -1648,7 +1678,7 @@ def _zero_m_n_batched_svd(
                 usm_type=usm_type,
                 sycl_queue=exec_q,
             )
-        return u, s, vt
+        return SVDResult(u, s, vt)
     return s
 
 
@@ -1692,7 +1722,7 @@ def _zero_m_n_svd(
             usm_type=usm_type,
             sycl_queue=exec_q,
         )
-        return u, s, vt
+        return SVDResult(u, s, vt)
     return s
 
 
@@ -1993,7 +2023,7 @@ def dpnp_det(a):
     return det.reshape(shape)
 
 
-def dpnp_eigh(a, UPLO, eigen_mode="V"):
+def dpnp_eigh(a, UPLO="L", eigen_mode="V"):
     """
     dpnp_eigh(a, UPLO, eigen_mode="V")
 
@@ -2016,7 +2046,7 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"):
         w = dpnp.empty_like(a, shape=a.shape[:-1], dtype=w_type)
         if eigen_mode == "V":
             v = dpnp.empty_like(a, dtype=v_type)
-            return w, v
+            return EighResult(w, v)
         return w
 
     if a.ndim > 2:
@@ -2097,7 +2127,7 @@ def dpnp_eigh(a, UPLO, eigen_mode="V"):
     else:
         out_v = v
 
-    return (w, out_v) if eigen_mode == "V" else w
+    return EighResult(w, out_v) if eigen_mode == "V" else w
 
 
 def dpnp_inv(a):
@@ -2546,7 +2576,7 @@ def dpnp_qr(a, mode="reduced"):
     r = a_t[:, :mc].transpose()
 
     r = _triu_inplace(r)
-    return (q, r)
+    return QRResult(q, r)
 
 
 def dpnp_solve(a, b):
@@ -2675,7 +2705,7 @@ def dpnp_slogdet(a):
             usm_type=a_usm_type,
             sycl_queue=a_sycl_queue,
         )
-        return sign, logdet
+        return SlogdetResult(sign, logdet)
 
     lu, ipiv, dev_info = _lu_factor(a, res_type)
 
@@ -2687,7 +2717,7 @@ def dpnp_slogdet(a):
 
     logdet = logdet.astype(logdet_dtype, copy=False)
     singular = dev_info > 0
-    return (
+    return SlogdetResult(
         dpnp.where(singular, res_type.type(0), sign).reshape(shape),
         dpnp.where(singular, logdet_dtype.type("-inf"), logdet).reshape(shape),
     )
@@ -2815,10 +2845,10 @@ def dpnp_svd(
         # For A^T = V S^T U^T, `u_h` becomes V and `vt_h` becomes U^T.
         # Transpose and swap them back to restore correct order for A.
         if trans_flag:
-            return vt_h.T, s_h, u_h.T
+            return SVDResult(vt_h.T, s_h, u_h.T)
         # gesvd call writes `u_h` and `vt_h` in Fortran order;
         # Convert to contiguous to align with NumPy
         u_h = dpnp.ascontiguousarray(u_h)
         vt_h = dpnp.ascontiguousarray(vt_h)
-        return u_h, s_h, vt_h
+        return SVDResult(u_h, s_h, vt_h)
     return s_h
diff --git a/dpnp/tests/test_linalg.py b/dpnp/tests/test_linalg.py
index 111c7845bdd8..d67e859c2675 100644
--- a/dpnp/tests/test_linalg.py
+++ b/dpnp/tests/test_linalg.py
@@ -521,7 +521,8 @@ def test_eigenvalues(self, func, shape, dtype, order):
         # we verify them through the eigen equation A*v=w*v.
         if func in ("eig", "eigh"):
             w, _ = getattr(numpy.linalg, func)(a)
-            w_dp, v_dp = getattr(dpnp.linalg, func)(a_dp)
+            result = getattr(dpnp.linalg, func)(a_dp)
+            w_dp, v_dp = result.eigenvalues, result.eigenvectors
 
             self.assert_eigen_decomposition(a_dp, w_dp, v_dp)
 
@@ -545,7 +546,8 @@ def test_eigenvalue_empty(self, func, shape, dtype):
 
         if func == "eig":
             w, v = getattr(numpy.linalg, func)(a_np)
-            w_dp, v_dp = getattr(dpnp.linalg, func)(a_dp)
+            result = getattr(dpnp.linalg, func)(a_dp)
+            w_dp, v_dp = result.eigenvalues, result.eigenvectors
 
             assert_dtype_allclose(v_dp, v)
 
@@ -2388,16 +2390,18 @@ def test_qr(self, dtype, shape, mode):
             dpnp_r = dpnp.linalg.qr(ia, mode)
         else:
             np_q, np_r = numpy.linalg.qr(a, mode)
-            dpnp_q, dpnp_r = dpnp.linalg.qr(ia, mode)
 
             # check decomposition
             if mode in ("complete", "reduced"):
+                result = dpnp.linalg.qr(ia, mode)
+                dpnp_q, dpnp_r = result.Q, result.R
                 assert_almost_equal(
                     dpnp.matmul(dpnp_q, dpnp_r),
                     a,
                     decimal=5,
                 )
             else:  # mode=="raw"
+                dpnp_q, dpnp_r = dpnp.linalg.qr(ia, mode)
                 assert_dtype_allclose(dpnp_q, np_q)
 
         if mode in ("raw", "r"):
@@ -2421,15 +2425,18 @@ def test_qr_large(self, dtype, shape, mode):
             dpnp_r = dpnp.linalg.qr(ia, mode)
         else:
             np_q, np_r = numpy.linalg.qr(a, mode)
-            dpnp_q, dpnp_r = dpnp.linalg.qr(ia, mode)
+
             # check decomposition
             if mode in ("complete", "reduced"):
+                result = dpnp.linalg.qr(ia, mode)
+                dpnp_q, dpnp_r = result.Q, result.R
                 assert_almost_equal(
                     dpnp.matmul(dpnp_q, dpnp_r),
                     a,
                     decimal=5,
                 )
             else:  # mode=="raw"
+                dpnp_q, dpnp_r = dpnp.linalg.qr(ia, mode)
                 assert_allclose(np_q, dpnp_q, atol=1e-4)
         if mode in ("raw", "r"):
             assert_allclose(np_r, dpnp_r, atol=1e-4)
@@ -2457,7 +2464,12 @@ def test_qr_empty(self, dtype, shape, mode):
             dpnp_r = dpnp.linalg.qr(ia, mode)
         else:
             np_q, np_r = numpy.linalg.qr(a, mode)
-            dpnp_q, dpnp_r = dpnp.linalg.qr(ia, mode)
+
+            if mode in ("complete", "reduced"):
+                result = dpnp.linalg.qr(ia, mode)
+                dpnp_q, dpnp_r = result.Q, result.R
+            else:
+                dpnp_q, dpnp_r = dpnp.linalg.qr(ia, mode)
 
             assert_dtype_allclose(dpnp_q, np_q)
 
@@ -2474,7 +2486,12 @@ def test_qr_strides(self, mode):
             dpnp_r = dpnp.linalg.qr(ia[::2, ::2], mode)
         else:
             np_q, np_r = numpy.linalg.qr(a[::2, ::2], mode)
-            dpnp_q, dpnp_r = dpnp.linalg.qr(ia[::2, ::2], mode)
+
+            if mode in ("complete", "reduced"):
+                result = dpnp.linalg.qr(ia[::2, ::2], mode)
+                dpnp_q, dpnp_r = result.Q, result.R
+            else:
+                dpnp_q, dpnp_r = dpnp.linalg.qr(ia[::2, ::2], mode)
 
             assert_dtype_allclose(dpnp_q, np_q)
 
@@ -2486,7 +2503,12 @@ def test_qr_strides(self, mode):
             dpnp_r = dpnp.linalg.qr(ia[::-2, ::-2], mode)
         else:
             np_q, np_r = numpy.linalg.qr(a[::-2, ::-2], mode)
-            dpnp_q, dpnp_r = dpnp.linalg.qr(ia[::-2, ::-2], mode)
+
+            if mode in ("complete", "reduced"):
+                result = dpnp.linalg.qr(ia[::-2, ::-2], mode)
+                dpnp_q, dpnp_r = result.Q, result.R
+            else:
+                dpnp_q, dpnp_r = dpnp.linalg.qr(ia[::-2, ::-2], mode)
 
             assert_dtype_allclose(dpnp_q, np_q)
 
@@ -2660,7 +2682,8 @@ def test_slogdet_2d(self, dtype):
         a_dp = dpnp.array(a_np)
 
         sign_expected, logdet_expected = numpy.linalg.slogdet(a_np)
-        sign_result, logdet_result = dpnp.linalg.slogdet(a_dp)
+        result = dpnp.linalg.slogdet(a_dp)
+        sign_result, logdet_result = result.sign, result.logabsdet
 
         assert_allclose(sign_expected, sign_result)
         assert_allclose(logdet_expected, logdet_result, rtol=1e-3, atol=1e-4)
@@ -2678,7 +2701,8 @@ def test_slogdet_3d(self, dtype):
         a_dp = dpnp.array(a_np)
 
         sign_expected, logdet_expected = numpy.linalg.slogdet(a_np)
-        sign_result, logdet_result = dpnp.linalg.slogdet(a_dp)
+        result = dpnp.linalg.slogdet(a_dp)
+        sign_result, logdet_result = result.sign, result.logabsdet
 
         assert_allclose(sign_expected, sign_result)
         assert_allclose(logdet_expected, logdet_result, rtol=1e-3, atol=1e-4)
@@ -2698,13 +2722,15 @@ def test_slogdet_strides(self):
 
         # positive strides
         sign_expected, logdet_expected = numpy.linalg.slogdet(a_np[::2, ::2])
-        sign_result, logdet_result = dpnp.linalg.slogdet(a_dp[::2, ::2])
+        result = dpnp.linalg.slogdet(a_dp[::2, ::2])
+        sign_result, logdet_result = result.sign, result.logabsdet
         assert_allclose(sign_expected, sign_result)
         assert_allclose(logdet_expected, logdet_result, rtol=1e-3, atol=1e-4)
 
         # negative strides
         sign_expected, logdet_expected = numpy.linalg.slogdet(a_np[::-2, ::-2])
-        sign_result, logdet_result = dpnp.linalg.slogdet(a_dp[::-2, ::-2])
+        result = dpnp.linalg.slogdet(a_dp[::-2, ::-2])
+        sign_result, logdet_result = result.sign, result.logabsdet
         assert_allclose(sign_expected, sign_result)
         assert_allclose(logdet_expected, logdet_result, rtol=1e-3, atol=1e-4)
 
@@ -2732,7 +2758,8 @@ def test_slogdet_singular_matrix(self, matrix):
         a_dp = dpnp.array(a_np)
 
         sign_expected, logdet_expected = numpy.linalg.slogdet(a_np)
-        sign_result, logdet_result = dpnp.linalg.slogdet(a_dp)
+        result = dpnp.linalg.slogdet(a_dp)
+        sign_result, logdet_result = result.sign, result.logabsdet
 
         assert_allclose(sign_expected, sign_result)
         assert_allclose(logdet_expected, logdet_result, rtol=1e-3, atol=1e-4)
@@ -2748,7 +2775,8 @@ def test_slogdet_singular_matrix_3D(self):
         a_dp = dpnp.array(a_np)
 
         sign_expected, logdet_expected = numpy.linalg.slogdet(a_np)
-        sign_result, logdet_result = dpnp.linalg.slogdet(a_dp)
+        result = dpnp.linalg.slogdet(a_dp)
+        sign_result, logdet_result = result.sign, result.logabsdet
 
         assert_allclose(sign_expected, sign_result)
         assert_allclose(logdet_expected, logdet_result, rtol=1e-3, atol=1e-4)
@@ -2841,13 +2869,14 @@ def test_svd(self, dtype, shape):
         a = numpy.arange(shape[0] * shape[1], dtype=dtype).reshape(shape)
         dp_a = dpnp.array(a)
 
-        np_u, np_s, np_vt = numpy.linalg.svd(a)
-        dp_u, dp_s, dp_vt = dpnp.linalg.svd(dp_a)
+        np_u, np_s, np_vh = numpy.linalg.svd(a)
+        result = dpnp.linalg.svd(dp_a)
+        dp_u, dp_s, dp_vh = result.U, result.S, result.Vh
 
-        self.check_types_shapes(dp_u, dp_s, dp_vt, np_u, np_s, np_vt)
+        self.check_types_shapes(dp_u, dp_s, dp_vh, np_u, np_s, np_vh)
         self.get_tol(dtype)
         self.check_decomposition(
-            dp_a, dp_u, dp_s, dp_vt, np_u, np_s, np_vt, True
+            dp_a, dp_u, dp_s, dp_vh, np_u, np_s, np_vh, True
         )
 
     @pytest.mark.parametrize("dtype", get_float_complex_dtypes())
@@ -2860,25 +2889,26 @@ def test_svd_hermitian(self, dtype, compute_vt, shape):
         dp_a = dpnp.array(a)
 
         if compute_vt:
-            np_u, np_s, np_vt = numpy.linalg.svd(
+            np_u, np_s, np_vh = numpy.linalg.svd(
                 a, compute_uv=compute_vt, hermitian=True
             )
-            dp_u, dp_s, dp_vt = dpnp.linalg.svd(
+            result = dpnp.linalg.svd(
                 dp_a, compute_uv=compute_vt, hermitian=True
             )
+            dp_u, dp_s, dp_vh = result.U, result.S, result.Vh
         else:
             np_s = numpy.linalg.svd(a, compute_uv=compute_vt, hermitian=True)
             dp_s = dpnp.linalg.svd(dp_a, compute_uv=compute_vt, hermitian=True)
-            np_u = np_vt = dp_u = dp_vt = None
+            np_u = np_vh = dp_u = dp_vh = None
 
         self.check_types_shapes(
-            dp_u, dp_s, dp_vt, np_u, np_s, np_vt, compute_vt
+            dp_u, dp_s, dp_vh, np_u, np_s, np_vh, compute_vt
         )
 
         self.get_tol(dtype)
 
         self.check_decomposition(
-            dp_a, dp_u, dp_s, dp_vt, np_u, np_s, np_vt, compute_vt
+            dp_a, dp_u, dp_s, dp_vh, np_u, np_s, np_vh, compute_vt
         )
 
     def test_svd_errors(self):