diff --git a/docarray/typing/tensor/__init__.py b/docarray/typing/tensor/__init__.py
index 58568a08c36..1a7e1f190b4 100644
--- a/docarray/typing/tensor/__init__.py
+++ b/docarray/typing/tensor/__init__.py
@@ -1,3 +1,5 @@
+from typing import Dict, Set, Type
+
 from docarray.typing.tensor.embedding import Embedding, NdArrayEmbedding
 from docarray.typing.tensor.ndarray import NdArray
 from docarray.typing.tensor.tensor import Tensor
@@ -7,8 +9,12 @@
     'Tensor',
     'Embedding',
     'NdArrayEmbedding',
+    'framework_types',
+    'type_to_framework',
 ]
 
+framework_types: Dict[str, Set] = {'numpy': {NdArray, NdArrayEmbedding}, 'torch': set()}
+
 try:
     import torch  # noqa: F401
 except ImportError:
@@ -18,3 +24,10 @@
     from docarray.typing.tensor.torch_tensor import TorchTensor  # noqa: F401
 
     __all__.extend(['TorchEmbedding', 'TorchTensor'])
+    framework_types['torch'] = {TorchTensor, TorchEmbedding}
+
+type_to_framework: Dict[Type, str] = {
+    type_: framework
+    for framework, type_set in framework_types.items()
+    for type_ in type_set
+}
diff --git a/docarray/utility/__init__.py b/docarray/utility/__init__.py
new file mode 100644
index 00000000000..7c9bada7262
--- /dev/null
+++ b/docarray/utility/__init__.py
@@ -0,0 +1,3 @@
+from docarray.utility.find import find, find_batched
+
+__all__ = ['find', 'find_batched']
diff --git a/docarray/utility/find.py b/docarray/utility/find.py
new file mode 100644
index 00000000000..dcf57395d4d
--- /dev/null
+++ b/docarray/utility/find.py
@@ -0,0 +1,321 @@
+import importlib
+from typing import Callable, List, NamedTuple, Optional, Type, Union
+
+import torch  # TODO(johannes) this breaks the optional import of torch
+
+from docarray import Document, DocumentArray
+from docarray.typing import Tensor
+from docarray.typing.tensor import type_to_framework
+
+# but will be fixed once we have a computational backend
+
+
+class FindResult(NamedTuple):
+    documents: DocumentArray
+    scores: Tensor
+
+
+def find(
+    index: DocumentArray,
+    query: Union[Tensor, Document],
+    embedding_field: str = 'embedding',
+    metric: str = 'cosine_sim',
+    limit: int = 10,
+    device: Optional[str] = None,
+    descending: Optional[bool] = None,
+) -> FindResult:
+    """
+    Find the closest Documents in the index to the query.
+    Supports PyTorch and NumPy embeddings.
+
+    .. note::
+        This utility function is likely to be removed once
+        Document Stores are available.
+        At that point, and in-memory Document Store will serve the same purpose
+        by exposing a .find() method.
+
+    .. note::
+        This is a simple implementation that assumes the same embedding field name for
+        both query and index, does not support nested search, and does not support
+        hybrid (multi-vector) search. These shortcoming will be addressed in future
+        versions.
+
+    EXAMPLE USAGE
+
+    .. code-block:: python
+
+        from docarray import DocumentArray, Document
+        from docarray.typing import TorchTensor
+        from docarray.utility.find import find
+
+
+        class MyDocument(Document):
+            embedding: TorchTensor
+
+
+        index = DocumentArray[MyDocument](
+            [MyDocument(embedding=torch.rand(128)) for _ in range(100)]
+        )
+
+        # use Document as query
+        query = MyDocument(embedding=torch.rand(128))
+        top_matches, scores = find(
+            index=index,
+            query=query,
+            embedding_field='tensor',
+            metric='cosine_sim',
+        )
+
+        # use tensor as query
+        query = torch.rand(128)
+        top_matches, scores = find(
+            index=index,
+            query=query,
+            embedding_field='tensor',
+            metric='cosine_sim',
+        )
+
+    :param index: the index of Documents to search in
+    :param query: the query to search for
+    :param embedding_field: the tensor-like field in the index to use
+        for the similarity computation
+    :param metric: the distance metric to use for the similarity computation.
+        Can be one of the following strings:
+        'cosine_sim' for cosine similarity, 'euclidean_dist' for euclidean distance,
+        'sqeuclidean_dist' for squared euclidean distance
+    :param limit: return the top `limit` results
+    :param device: the computational device to use,
+        can be either `cpu` or a `cuda` device.
+    :param descending: sort the results in descending order.
+        Per default, this is chosen based on the `metric` argument.
+    :return: A named tuple of the form (DocumentArray, Tensor),
+        where the first element contains the closes matches for the query,
+        and the second element contains the corresponding scores.
+    """
+    query = _extract_embedding_single(query, embedding_field)
+    return find_batched(
+        index=index,
+        query=query,
+        embedding_field=embedding_field,
+        metric=metric,
+        limit=limit,
+        device=device,
+        descending=descending,
+    )[0]
+
+
+def find_batched(
+    index: DocumentArray,
+    query: Union[Tensor, DocumentArray],
+    embedding_field: str = 'embedding',
+    metric: str = 'cosine_sim',
+    limit: int = 10,
+    device: Optional[str] = None,
+    descending: Optional[bool] = None,
+) -> List[FindResult]:
+    """
+    Find the closest Documents in the index to the queries.
+    Supports PyTorch and NumPy embeddings.
+
+    .. note::
+        This utility function is likely to be removed once
+        Document Stores are available.
+        At that point, and in-memory Document Store will serve the same purpose
+        by exposing a .find() method.
+
+    .. note::
+        This is a simple implementation that assumes the same embedding field name for
+        both query and index, does not support nested search, and does not support
+        hybrid (multi-vector) search. These shortcoming will be addressed in future
+        versions.
+
+        EXAMPLE USAGE
+
+    .. code-block:: python
+
+        from docarray import DocumentArray, Document
+        from docarray.typing import TorchTensor
+        from docarray.utility.find import find
+
+
+        class MyDocument(Document):
+            embedding: TorchTensor
+
+
+        index = DocumentArray[MyDocument](
+            [MyDocument(embedding=torch.rand(128)) for _ in range(100)]
+        )
+
+        # use DocumentArray as query
+        query = DocumentArray[MyDocument](
+            [MyDocument(embedding=torch.rand(128)) for _ in range(3)]
+        )
+        results = find(
+            index=index,
+            query=query,
+            embedding_field='tensor',
+            metric='cosine_sim',
+        )
+        top_matches, scores = results[0]
+
+        # use tensor as query
+        query = torch.rand(3, 128)
+        results, scores = find(
+            index=index,
+            query=query,
+            embedding_field='tensor',
+            metric='cosine_sim',
+        )
+        top_matches, scores = results[0]
+
+    :param index: the index of Documents to search in
+    :param query: the query to search for
+    :param embedding_field: the tensor-like field in the index to use
+        for the similarity computation
+    :param metric: the distance metric to use for the similarity computation.
+        Can be one of the following strings:
+        'cosine_sim' for cosine similarity, 'euclidean_dist' for euclidean distance,
+        'sqeuclidean_dist' for squared euclidean distance
+    :param limit: return the top `limit` results
+    :param device: the computational device to use,
+        can be either `cpu` or a `cuda` device.
+    :param descending: sort the results in descending order.
+        Per default, this is chosen based on the `metric` argument.
+    :return: a list of named tuples of the form (DocumentArray, Tensor),
+        where the first element contains the closes matches for each query,
+        and the second element contains the corresponding scores.
+    """
+    if descending is None:
+        descending = metric.endswith('_sim')  # similarity metrics are descending
+
+    embedding_type = _da_attr_type(index, embedding_field)
+
+    # get framework-specific distance and top_k function
+    metric_fn = _get_metric_fn(embedding_type, metric)
+    top_k_fn = _get_topk_fn(embedding_type)
+
+    # extract embeddings from query and index
+    index_embeddings = _extraxt_embeddings(index, embedding_field, embedding_type)
+    query_embeddings = _extraxt_embeddings(query, embedding_field, embedding_type)
+
+    # compute distances and return top results
+    dists = metric_fn(query_embeddings, index_embeddings, device=device)
+    top_scores, top_indices = top_k_fn(
+        dists, k=limit, device=device, descending=descending
+    )
+
+    index_doc_type = index.document_type
+    results = []
+    for indices_per_query, scores_per_query in zip(top_indices, top_scores):
+        docs_per_query = DocumentArray[index_doc_type]([])  # type: ignore
+        for idx in indices_per_query:  # workaround until #930 is fixed
+            docs_per_query.append(index[idx])
+        results.append(FindResult(scores=scores_per_query, documents=docs_per_query))
+    return results
+
+
+def _extract_embedding_single(
+    data: Union[DocumentArray, Document, Tensor],
+    embedding_field: str,
+) -> Tensor:
+    """Extract the embeddings from a single query,
+    and return it in a batched representation.
+
+    :param data: the data
+    :param embedding_field: the embedding field
+    :param embedding_type: type of the embedding: torch.Tensor, numpy.ndarray etc.
+    :return: the embeddings
+    """
+    if isinstance(data, Document):
+        emb = getattr(data, embedding_field)
+    else:  # treat data as tensor
+        emb = data
+    if len(emb.shape) == 1:
+        # TODO(johannes) solve this with computational backend,
+        #  this is ugly hack for now
+        if isinstance(emb, torch.Tensor):
+            emb = emb.unsqueeze(0)
+        else:
+            import numpy as np
+
+            if isinstance(emb, np.ndarray):
+                emb = np.expand_dims(emb, axis=0)
+    return emb
+
+
+def _extraxt_embeddings(
+    data: Union[DocumentArray, Document, Tensor],
+    embedding_field: str,
+    embedding_type: Type,
+) -> Tensor:
+    """Extract the embeddings from the data.
+
+    :param data: the data
+    :param embedding_field: the embedding field
+    :param embedding_type: type of the embedding: torch.Tensor, numpy.ndarray etc.
+    :return: the embeddings
+    """
+    # TODO(johannes) put docarray stack in the computational backend
+    if isinstance(data, DocumentArray):
+        emb = getattr(data, embedding_field)
+        if not data.is_stacked():
+            emb = embedding_type.__docarray_stack__(emb)
+    elif isinstance(data, Document):
+        emb = getattr(data, embedding_field)
+    else:  # treat data as tensor
+        emb = data
+
+    if len(emb.shape) == 1:
+        # TODO(johannes) solve this with computational backend,
+        #  this is ugly hack for now
+        if isinstance(emb, torch.Tensor):
+            emb = emb.unsqueeze(0)
+        else:
+            import numpy as np
+
+            if isinstance(emb, np.ndarray):
+                emb = np.expand_dims(emb, axis=0)
+    return emb
+
+
+def _da_attr_type(da: DocumentArray, attr: str) -> Type:
+    """Get the type of the attribute according to the Document type
+    (schema) of the DocumentArray.
+
+    :param da: the DocumentArray
+    :param attr: the attribute name
+    :return: the type of the attribute
+    """
+    return da.document_type.__fields__[attr].type_
+
+
+def _get_topk_fn(embedding_type: Type) -> Callable:
+    """Dynamically import the distance function from the framework-specific module.
+    This will go away once we have a computational backend.
+
+    :param embedding_type: the type of the embedding
+    :param distance_name: the name of the distance function
+    :return: the framework-specific distance function
+    """
+    framework = type_to_framework[embedding_type]
+    return getattr(
+        importlib.import_module(f'docarray.utility.helper.{framework}'),
+        'top_k',
+    )
+
+
+def _get_metric_fn(embedding_type: Type, metric: Union[str, Callable]) -> Callable:
+    """Dynamically import the distance function from the framework-specific module.
+    This will go away once we have a proper computational backend.
+
+    :param embedding_type: the type of the embedding
+    :param metric: the name of the metric, or the metric itself
+    :return: the framework-specific metric
+    """
+    if callable(metric):
+        return metric
+    framework = type_to_framework[embedding_type]
+    return getattr(
+        importlib.import_module(f'docarray.utility.math.metrics.{framework}'),
+        f'{metric}',
+    )
diff --git a/docarray/utility/helper/__init__.py b/docarray/utility/helper/__init__.py
new file mode 100644
index 00000000000..e69de29bb2d
diff --git a/docarray/utility/helper/numpy.py b/docarray/utility/helper/numpy.py
new file mode 100644
index 00000000000..cf2dcadc222
--- /dev/null
+++ b/docarray/utility/helper/numpy.py
@@ -0,0 +1,51 @@
+import warnings
+from typing import Optional, Tuple
+
+import numpy as np
+
+
+def top_k(
+    values: 'np.ndarray',
+    k: int,
+    descending: bool = False,
+    device: Optional[str] = None,
+) -> Tuple['np.ndarray', 'np.ndarray']:
+    """
+    Retrieves the top k smallest values in `values`,
+    and returns them alongside their indices in the input `values`.
+    Can also be used to retrieve the top k largest values,
+    by setting the `descending` flag.
+
+    :param values: Torch tensor of values to rank.
+        Should be of shape (n_queries, n_values_per_query).
+        Inputs of shape (n_values_per_query,) will be expanded
+        to (1, n_values_per_query).
+    :param k: number of values to retrieve
+    :param descending: retrieve largest values instead of smallest values
+    :param device: Not supported for this backend
+    :return: Tuple containing the retrieved values, and their indices.
+        Both ar of shape (n_queries, k)
+    """
+    if device is not None:
+        warnings.warn('`device` is not supported for numpy operations')
+
+    if len(values.shape) == 1:
+        values = np.expand_dims(values, axis=0)
+
+    if descending:
+        values = -values
+
+    if k >= values.shape[1]:
+        idx = values.argsort(axis=1)[:, :k]
+        values = np.take_along_axis(values, idx, axis=1)
+    else:
+        idx_ps = values.argpartition(kth=k, axis=1)[:, :k]
+        values = np.take_along_axis(values, idx_ps, axis=1)
+        idx_fs = values.argsort(axis=1)
+        idx = np.take_along_axis(idx_ps, idx_fs, axis=1)
+        values = np.take_along_axis(values, idx_fs, axis=1)
+
+    if descending:
+        values = -values
+
+    return values, idx
diff --git a/docarray/utility/helper/torch.py b/docarray/utility/helper/torch.py
new file mode 100644
index 00000000000..1e5bb92722a
--- /dev/null
+++ b/docarray/utility/helper/torch.py
@@ -0,0 +1,35 @@
+from typing import Optional, Tuple
+
+import torch
+
+
+def top_k(
+    values: torch.Tensor,
+    k: int,
+    descending: bool = False,
+    device: Optional[str] = None,
+) -> Tuple[torch.Tensor, torch.Tensor]:
+    """
+    Retrieves the top k smallest values in `values`,
+    and returns them alongside their indices in the input `values`.
+    Can also be used to retrieve the top k largest values,
+    by setting the `descending` flag.
+
+    :param values: Torch tensor of values to rank.
+        Should be of shape (n_queries, n_values_per_query).
+        Inputs of shape (n_values_per_query,) will be expanded
+        to (1, n_values_per_query).
+    :param k: number of values to retrieve
+    :param descending: retrieve largest values instead of smallest values
+    :param device: the computational device to use,
+        can be either `cpu` or a `cuda` device.
+    :return: Tuple containing the retrieved values, and their indices.
+        Both ar of shape (n_queries, k)
+    """
+    if device is not None:
+        values = values.to(device)
+    if len(values.shape) <= 1:
+        values = values.view(1, -1)
+    len_values = values.shape[-1] if len(values.shape) > 1 else len(values)
+    k = min(k, len_values)
+    return torch.topk(input=values, k=k, largest=descending, sorted=True, dim=-1)
diff --git a/docarray/utility/math/__init__.py b/docarray/utility/math/__init__.py
new file mode 100644
index 00000000000..e69de29bb2d
diff --git a/docarray/utility/math/metrics/__init__.py b/docarray/utility/math/metrics/__init__.py
new file mode 100644
index 00000000000..e69de29bb2d
diff --git a/docarray/utility/math/metrics/numpy.py b/docarray/utility/math/metrics/numpy.py
new file mode 100644
index 00000000000..0efe5ca40da
--- /dev/null
+++ b/docarray/utility/math/metrics/numpy.py
@@ -0,0 +1,123 @@
+import warnings
+from typing import List, Optional
+
+import numpy as np
+
+
+def _expand_if_single_axis(*matrices: np.ndarray) -> List[np.ndarray]:
+    """Expands arrays that only have one axis, at dim 0.
+    This ensures that all outputs can be treated as matrices, not vectors.
+
+    :param matrices: Matrices to be expanded
+    :return: List of the input matrices,
+        where single axis matrices are expanded at dim 0.
+    """
+    expanded = []
+    for m in matrices:
+        if len(m.shape) == 1:
+            expanded.append(np.expand_dims(m, axis=0))
+        else:
+            expanded.append(m)
+    return expanded
+
+
+def _expand_if_scalar(arr: np.ndarray) -> np.ndarray:
+    if len(arr.shape) == 0:  # avoid scalar output
+        arr = np.expand_dims(arr, axis=0)
+    return arr
+
+
+def cosine_sim(
+    x_mat: np.ndarray,
+    y_mat: np.ndarray,
+    eps: float = 1e-7,
+    device: Optional[str] = None,
+) -> np.ndarray:
+    """Pairwise cosine similarities between all vectors in x_mat and y_mat.
+
+    :param x_mat: np.ndarray of shape (n_vectors, n_dim), where n_vectors is the
+        number of vectors and n_dim is the number of dimensions of each example.
+    :param y_mat: np.ndarray of shape (n_vectors, n_dim), where n_vectors is the
+        number of vectors and n_dim is the number of dimensions of each example.
+    :param eps: a small jitter to avoid divde by zero
+    :param device: Not supported for this backend
+    :return: np.ndarray  of shape (n_vectors, n_vectors) containing all pairwise
+        cosine distances.
+        The index [i_x, i_y] contains the cosine distance between
+        x_mat[i_x] and y_mat[i_y].
+    """
+    if device is not None:
+        warnings.warn('`device` is not supported for numpy operations')
+
+    x_mat, y_mat = _expand_if_single_axis(x_mat, y_mat)
+
+    sims = np.clip(
+        (np.dot(x_mat, y_mat.T) + eps)
+        / (
+            np.outer(np.linalg.norm(x_mat, axis=1), np.linalg.norm(y_mat, axis=1)) + eps
+        ),
+        -1,
+        1,
+    ).squeeze()
+    return _expand_if_scalar(sims)
+
+
+def sqeuclidean_dist(
+    x_mat: np.ndarray,
+    y_mat: np.ndarray,
+    eps: float = 1e-7,
+    device: Optional[str] = None,
+) -> np.ndarray:
+    """Pairwise Squared Euclidian distances between all vectors in x_mat and y_mat.
+
+    :param x_mat: np.ndarray of shape (n_vectors, n_dim), where n_vectors is the
+        number of vectors and n_dim is the number of dimensions of each example.
+    :param y_mat: np.ndarray of shape (n_vectors, n_dim), where n_vectors is the
+        number of vectors and n_dim is the number of dimensions of each example.
+    :param eps: a small jitter to avoid divde by zero
+    :param eps: Negative results (due to numerical inaccuracies) no smaller than
+        -eps will be returned as 0.
+    :param device: Not supported for this backend
+    :return: np.ndarray  of shape (n_vectors, n_vectors) containing all pairwise
+        Squared Euclidian distances.
+        The index [i_x, i_y] contains the cosine Squared Euclidian between
+        x_mat[i_x] and y_mat[i_y].
+    """
+    if device is not None:
+        warnings.warn('`device` is not supported for numpy operations')
+
+    x_mat, y_mat = _expand_if_single_axis(x_mat, y_mat)
+
+    dists = (
+        np.sum(y_mat**2, axis=1)
+        + np.sum(x_mat**2, axis=1)[:, np.newaxis]
+        - 2 * np.dot(x_mat, y_mat.T)
+    ).squeeze()
+
+    # remove numerical artifacts
+    dists = np.where(np.logical_and(dists < 0, dists > -eps), 0, dists)
+    return _expand_if_scalar(dists)
+
+
+def euclidean_dist(
+    x_mat: np.ndarray, y_mat: np.ndarray, device: Optional[str] = None
+) -> np.ndarray:
+    """Pairwise Euclidian distances between all vectors in x_mat and y_mat.
+
+    :param x_mat: np.ndarray of shape (n_vectors, n_dim), where n_vectors is the
+        number of vectors and n_dim is the number of dimensions of each example.
+    :param y_mat: np.ndarray of shape (n_vectors, n_dim), where n_vectors is the
+        number of vectors and n_dim is the number of dimensions of each example.
+    :param eps: a small jitter to avoid divde by zero
+    :param device: Not supported for this backend
+    :return: np.ndarray  of shape (n_vectors, n_vectors) containing all pairwise
+        euclidian distances.
+        The index [i_x, i_y] contains the euclidian distance between
+        x_mat[i_x] and y_mat[i_y].
+    """
+    if device is not None:
+        warnings.warn('`device` is not supported for numpy operations')
+
+    x_mat, y_mat = _expand_if_single_axis(x_mat, y_mat)
+
+    return _expand_if_scalar(np.sqrt(sqeuclidean_dist(x_mat, y_mat)).squeeze())
diff --git a/docarray/utility/math/metrics/torch.py b/docarray/utility/math/metrics/torch.py
new file mode 100644
index 00000000000..2e1b12a50f8
--- /dev/null
+++ b/docarray/utility/math/metrics/torch.py
@@ -0,0 +1,114 @@
+from typing import List, Optional
+
+import torch
+
+
+def _unsqueeze_if_single_axis(*matrices: torch.Tensor) -> List[torch.Tensor]:
+    """Unsqueezes tensors that only have one axis, at dim 0.
+    This ensures that all outputs can be treated as matrices, not vectors.
+
+    :param matrices: Matrices to be unsqueezed
+    :return: List of the input matrices,
+        where single axis matrices are unsqueezed at dim 0.
+    """
+    unsqueezed = []
+    for m in matrices:
+        if len(m.shape) == 1:
+            unsqueezed.append(m.unsqueeze(0))
+        else:
+            unsqueezed.append(m)
+    return unsqueezed
+
+
+def _usqueeze_if_scalar(t: torch.Tensor):
+    if len(t.shape) == 0:  # avoid scalar output
+        t = t.unsqueeze(0)
+    return t
+
+
+def cosine_sim(
+    x_mat: torch.Tensor,
+    y_mat: torch.Tensor,
+    eps: float = 1e-7,
+    device: Optional[str] = None,
+) -> torch.Tensor:
+    """Pairwise cosine similarities between all vectors in x_mat and y_mat.
+
+    :param x_mat: tensor of shape (n_vectors, n_dim), where n_vectors is the
+        number of vectors and n_dim is the number of dimensions of each example.
+    :param y_mat: tensor of shape (n_vectors, n_dim), where n_vectors is the
+        number of vectors and n_dim is the number of dimensions of each example.
+    :param eps: a small jitter to avoid divde by zero
+    :param device: the device to use for pytorch computations.
+        Either 'cpu' or a 'cuda' device.
+        If not provided, the devices of x_mat and y_mat are used.
+    :return: np.ndarray  of shape (n_vectors, n_vectors) containing all pairwise
+        cosine distances.
+        The index [i_x, i_y] contains the cosine distance between
+        x_mat[i_x] and y_mat[i_y].
+    """
+    if device is not None:
+        x_mat = x_mat.to(device)
+        y_mat = y_mat.to(device)
+
+    x_mat, y_mat = _unsqueeze_if_single_axis(x_mat, y_mat)
+
+    a_n, b_n = x_mat.norm(dim=1)[:, None], y_mat.norm(dim=1)[:, None]
+    a_norm = x_mat / torch.clamp(a_n, min=eps)
+    b_norm = y_mat / torch.clamp(b_n, min=eps)
+    sims = torch.mm(a_norm, b_norm.transpose(0, 1)).squeeze()
+    return _usqueeze_if_scalar(sims)
+
+
+def euclidean_dist(
+    x_mat: torch.Tensor, y_mat: torch.Tensor, device: Optional[str] = None
+) -> torch.Tensor:
+    """Pairwise Euclidian distances between all vectors in x_mat and y_mat.
+
+    :param x_mat: tensor of shape (n_vectors, n_dim), where n_vectors is the
+        number of vectors and n_dim is the number of dimensions of each example.
+    :param y_mat: tensor of shape (n_vectors, n_dim), where n_vectors is the
+        number of vectors and n_dim is the number of dimensions of each example.
+    :param device: the device to use for pytorch computations.
+        Either 'cpu' or a 'cuda' device.
+        If not provided, the devices of x_mat and y_mat are used.
+    :return: np.ndarray  of shape (n_vectors, n_vectors) containing all pairwise
+        euclidian distances.
+        The index [i_x, i_y] contains the euclidian distance between
+        x_mat[i_x] and y_mat[i_y].
+    """
+    if device is not None:
+        x_mat = x_mat.to(device)
+        y_mat = y_mat.to(device)
+
+    x_mat, y_mat = _unsqueeze_if_single_axis(x_mat, y_mat)
+
+    dists = torch.cdist(x_mat, y_mat).squeeze()
+    return _usqueeze_if_scalar(dists)
+
+
+def sqeuclidean_dist(
+    x_mat: torch.Tensor, y_mat: torch.Tensor, device: Optional[str] = None
+) -> torch.Tensor:
+    """Pairwise Squared Euclidian distances between all vectors in x_mat and y_mat.
+
+    :param x_mat: tensor of shape (n_vectors, n_dim), where n_vectors is the
+        number of vectors and n_dim is the number of dimensions of each example.
+    :param y_mat: tensor of shape (n_vectors, n_dim), where n_vectors is the
+        number of vectors and n_dim is the number of dimensions of each example.
+    :param eps: a small jitter to avoid divde by zero
+    :param device: the device to use for pytorch computations.
+        Either 'cpu' or a 'cuda' device.
+        If not provided, the devices of x_mat and y_mat are used.
+    :return: np.ndarray  of shape (n_vectors, n_vectors) containing all pairwise
+        Squared Euclidian distances.
+        The index [i_x, i_y] contains the cosine Squared Euclidian between
+        x_mat[i_x] and y_mat[i_y].
+    """
+    if device is not None:
+        x_mat = x_mat.to(device)
+        y_mat = y_mat.to(device)
+
+    x_mat, y_mat = _unsqueeze_if_single_axis(x_mat, y_mat)
+
+    return _usqueeze_if_scalar((torch.cdist(x_mat, y_mat) ** 2).squeeze())
diff --git a/tests/units/util/__init__.py b/tests/units/util/__init__.py
new file mode 100644
index 00000000000..e69de29bb2d
diff --git a/tests/units/util/test_find.py b/tests/units/util/test_find.py
new file mode 100644
index 00000000000..665b880912c
--- /dev/null
+++ b/tests/units/util/test_find.py
@@ -0,0 +1,254 @@
+import numpy as np
+import pytest
+import torch
+
+from docarray import Document, DocumentArray
+from docarray.typing import NdArray, TorchTensor
+from docarray.utility import find, find_batched
+
+
+class TorchDoc(Document):
+    tensor: TorchTensor
+
+
+class NdDoc(Document):
+    tensor: NdArray
+
+
+@pytest.fixture()
+def random_torch_query():
+    return TorchDoc(tensor=torch.rand(128))
+
+
+@pytest.fixture()
+def random_torch_batch_query():
+    return DocumentArray[TorchDoc]([TorchDoc(tensor=torch.rand(128)) for _ in range(5)])
+
+
+@pytest.fixture()
+def random_nd_query():
+    return NdDoc(tensor=np.random.rand(128))
+
+
+@pytest.fixture()
+def random_nd_batch_query():
+    return DocumentArray[NdDoc]([NdDoc(tensor=np.random.rand(128)) for _ in range(5)])
+
+
+@pytest.fixture()
+def random_torch_index():
+    return DocumentArray[TorchDoc](TorchDoc(tensor=torch.rand(128)) for _ in range(10))
+
+
+@pytest.fixture()
+def random_nd_index():
+    return DocumentArray[NdDoc](NdDoc(tensor=np.random.rand(128)) for _ in range(10))
+
+
+@pytest.mark.parametrize('metric', ['cosine_sim', 'euclidean_dist', 'sqeuclidean_dist'])
+def test_find_torch(random_torch_query, random_torch_index, metric):
+    top_k, scores = find(
+        random_torch_index,
+        random_torch_query,
+        embedding_field='tensor',
+        limit=7,
+        metric=metric,
+    )
+    assert len(top_k) == 7
+    assert len(scores) == 7
+    if metric.endswith('_dist'):
+        assert (torch.stack(sorted(scores)) == scores).all()
+    else:
+        assert (torch.stack(sorted(scores, reverse=True)) == scores).all()
+
+
+def test_find_torch_tensor_query(random_torch_query, random_torch_index):
+    query = random_torch_query.tensor
+    top_k, scores = find(
+        random_torch_index,
+        query,
+        embedding_field='tensor',
+        limit=7,
+        metric='cosine_sim',
+    )
+    assert len(top_k) == 7
+    assert len(scores) == 7
+    assert (torch.stack(sorted(scores, reverse=True)) == scores).all()
+
+
+def test_find_torch_stacked(random_torch_query, random_torch_index):
+    random_torch_index.stack()
+    top_k, scores = find(
+        random_torch_index,
+        random_torch_query,
+        embedding_field='tensor',
+        limit=7,
+        metric='cosine_sim',
+    )
+    assert len(top_k) == 7
+    assert len(scores) == 7
+    assert (torch.stack(sorted(scores, reverse=True)) == scores).all()
+
+
+@pytest.mark.parametrize('metric', ['cosine_sim', 'euclidean_dist', 'sqeuclidean_dist'])
+def test_find_np(random_nd_query, random_nd_index, metric):
+    top_k, scores = find(
+        random_nd_index,
+        random_nd_query,
+        embedding_field='tensor',
+        limit=7,
+        metric=metric,
+    )
+    assert len(top_k) == 7
+    assert len(scores) == 7
+    if metric.endswith('_dist'):
+        assert (sorted(scores) == scores).all()
+    else:
+        assert (sorted(scores, reverse=True) == scores).all()
+
+
+def test_find_np_tensor_query(random_nd_query, random_nd_index):
+    query = random_nd_query.tensor
+    top_k, scores = find(
+        random_nd_index,
+        query,
+        embedding_field='tensor',
+        limit=7,
+        metric='cosine_sim',
+    )
+    assert len(top_k) == 7
+    assert len(scores) == 7
+    assert (sorted(scores, reverse=True) == scores).all()
+
+
+def test_find_np_stacked(random_nd_query, random_nd_index):
+    random_nd_index.stack()
+    top_k, scores = find(
+        random_nd_index,
+        random_nd_query,
+        embedding_field='tensor',
+        limit=7,
+        metric='cosine_sim',
+    )
+    assert len(top_k) == 7
+    assert len(scores) == 7
+    assert (sorted(scores, reverse=True) == scores).all()
+
+
+@pytest.mark.parametrize('metric', ['cosine_sim', 'euclidean_dist', 'sqeuclidean_dist'])
+def test_find_batched_torch(random_torch_batch_query, random_torch_index, metric):
+    results = find_batched(
+        random_torch_index,
+        random_torch_batch_query,
+        embedding_field='tensor',
+        limit=7,
+        metric=metric,
+    )
+    assert len(results) == len(random_torch_batch_query)
+    for top_k, scores in results:
+        assert len(top_k) == 7
+        assert len(scores) == 7
+    for sc in [scores for _, scores in results]:
+        if metric.endswith('_dist'):
+            assert (torch.stack(sorted(sc)) == sc).all()
+        else:
+            assert (torch.stack(sorted(sc, reverse=True)) == sc).all()
+
+
+def test_find_batched_torch_tensor_query(random_torch_batch_query, random_torch_index):
+    query = torch.stack(random_torch_batch_query.tensor)
+    results = find_batched(
+        random_torch_index,
+        query,
+        embedding_field='tensor',
+        limit=7,
+        metric='cosine_sim',
+    )
+    assert len(results) == len(random_torch_batch_query)
+    for top_k, scores in results:
+        assert len(top_k) == 7
+        assert len(scores) == 7
+    for sc in [scores for _, scores in results]:
+        assert (torch.stack(sorted(sc, reverse=True)) == sc).all()
+
+
+@pytest.mark.parametrize('stack_what', ['index', 'query', 'both'])
+def test_find_batched_torch_stacked(
+    random_torch_batch_query, random_torch_index, stack_what
+):
+    if stack_what in ('index', 'both'):
+        random_torch_index.stack()
+    if stack_what in ('query', 'both'):
+        random_torch_batch_query.stack()
+
+    results = find_batched(
+        random_torch_index,
+        random_torch_batch_query,
+        embedding_field='tensor',
+        limit=7,
+        metric='cosine_sim',
+    )
+    assert len(results) == len(random_torch_batch_query)
+    for top_k, scores in results:
+        assert len(top_k) == 7
+        assert len(scores) == 7
+    for sc in [scores for _, scores in results]:
+        assert (torch.stack(sorted(sc, reverse=True)) == sc).all()
+
+
+@pytest.mark.parametrize('metric', ['cosine_sim', 'euclidean_dist', 'sqeuclidean_dist'])
+def test_find_batched_np(random_nd_batch_query, random_nd_index, metric):
+    results = find_batched(
+        random_nd_index,
+        random_nd_batch_query,
+        embedding_field='tensor',
+        limit=7,
+        metric=metric,
+    )
+    assert len(results) == len(random_nd_batch_query)
+    for top_k, scores in results:
+        assert len(top_k) == 7
+        assert len(scores) == 7
+    for sc in [scores for _, scores in results]:
+        if metric.endswith('_dist'):
+            assert (sorted(sc) == sc).all()
+        else:
+            assert (sorted(sc, reverse=True) == sc).all()
+
+
+def test_find_batched_np_tensor_query(random_nd_batch_query, random_nd_index):
+    query = np.stack(random_nd_batch_query.tensor)
+    results = find_batched(
+        random_nd_index,
+        query,
+        embedding_field='tensor',
+        limit=7,
+        metric='cosine_sim',
+    )
+    assert len(results) == len(random_nd_batch_query)
+    for top_k, scores in results:
+        assert len(top_k) == 7
+        assert len(scores) == 7
+    for sc in [scores for _, scores in results]:
+        assert (sorted(sc, reverse=True) == sc).all()
+
+
+@pytest.mark.parametrize('stack_what', ['index', 'query', 'both'])
+def test_find_batched_np_stacked(random_nd_batch_query, random_nd_index, stack_what):
+    if stack_what in ('index', 'both'):
+        random_nd_index.stack()
+    if stack_what in ('query', 'both'):
+        random_nd_batch_query.stack()
+    results = find_batched(
+        random_nd_index,
+        random_nd_batch_query,
+        embedding_field='tensor',
+        limit=7,
+        metric='cosine_sim',
+    )
+    assert len(results) == len(random_nd_batch_query)
+    for top_k, scores in results:
+        assert len(top_k) == 7
+        assert len(scores) == 7
+    for sc in [scores for _, scores in results]:
+        assert (sorted(sc, reverse=True) == sc).all()
diff --git a/tests/units/util/test_metrics.py b/tests/units/util/test_metrics.py
new file mode 100644
index 00000000000..c50af26fc76
--- /dev/null
+++ b/tests/units/util/test_metrics.py
@@ -0,0 +1,169 @@
+import numpy as np
+import torch
+
+import docarray.utility.math.metrics.numpy as np_distances
+import docarray.utility.math.metrics.torch as torch_distances
+
+
+def test_cosine_sim_torch():
+    a = torch.rand(128)
+    b = torch.rand(128)
+    assert torch_distances.cosine_sim(a, b).shape == (1,)
+    assert torch_distances.cosine_sim(a, b) == torch_distances.cosine_sim(b, a)
+    torch.testing.assert_close(torch_distances.cosine_sim(a, a), torch.ones(1))
+
+    a = torch.rand(10, 3)
+    b = torch.rand(5, 3)
+    assert torch_distances.cosine_sim(a, b).shape == (10, 5)
+    assert torch_distances.cosine_sim(b, a).shape == (5, 10)
+    diag_dists = torch.diagonal(torch_distances.cosine_sim(b, b))  # self-comparisons
+    torch.testing.assert_allclose(diag_dists, torch.ones(5))
+
+
+def test_cosine_sim_np():
+    a = np.random.rand(128)
+    b = np.random.rand(128)
+    assert np_distances.cosine_sim(a, b).shape == (1,)
+    assert np_distances.cosine_sim(a, b) == np_distances.cosine_sim(b, a)
+    np.testing.assert_array_almost_equal(np_distances.cosine_sim(a, a), np.ones((1,)))
+
+    a = np.random.rand(10, 3)
+    b = np.random.rand(5, 3)
+    assert np_distances.cosine_sim(a, b).shape == (10, 5)
+    assert np_distances.cosine_sim(b, a).shape == (5, 10)
+    diag_dists = np.diagonal(np_distances.cosine_sim(b, b))  # self-comparisons
+    np.testing.assert_array_almost_equal(diag_dists, np.ones((5,)))
+
+
+def test_cosine_sim_compare():
+    a = torch.rand(128)
+    b = torch.rand(128)
+    torch.testing.assert_close(
+        torch_distances.cosine_sim(a, b),
+        torch.from_numpy(np_distances.cosine_sim(a.numpy(), b.numpy())),
+    )
+
+    a = torch.rand(10, 3)
+    b = torch.rand(5, 3)
+    torch.testing.assert_close(
+        torch_distances.cosine_sim(a, b),
+        torch.from_numpy(np_distances.cosine_sim(a.numpy(), b.numpy())),
+    )
+
+
+def test_euclidean_dist_torch():
+    a = torch.rand(128)
+    b = torch.rand(128)
+    assert torch_distances.euclidean_dist(a, b).shape == (1,)
+    assert torch_distances.euclidean_dist(a, b) == torch_distances.euclidean_dist(b, a)
+    torch.testing.assert_close(torch_distances.euclidean_dist(a, a), torch.zeros(1))
+
+    a = torch.rand(10, 3)
+    b = torch.rand(5, 3)
+    assert torch_distances.euclidean_dist(a, b).shape == (10, 5)
+    assert torch_distances.euclidean_dist(b, a).shape == (5, 10)
+    diag_dists = torch.diagonal(
+        torch_distances.euclidean_dist(b, b)
+    )  # self-comparisons
+    torch.testing.assert_allclose(diag_dists, torch.zeros(5))
+
+    a = torch.tensor([0.0, 2.0, 0.0])
+    b = torch.tensor([0.0, 0.0, 2.0])
+    desired_output_singleton = torch.sqrt(torch.tensor([2.0**2.0 + 2.0**2.0]))
+    torch.testing.assert_close(
+        torch_distances.euclidean_dist(a, b), desired_output_singleton
+    )
+
+    a = torch.tensor([[0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
+    b = torch.tensor([[0.0, 0.0, 2.0], [0.0, 2.0, 0.0]])
+    desired_output_singleton = torch.tensor(
+        [[desired_output_singleton.item(), 0.0], [0.0, desired_output_singleton.item()]]
+    )
+    torch.testing.assert_close(
+        torch_distances.euclidean_dist(a, b), desired_output_singleton
+    )
+
+
+def test_euclidean_dist_np():
+    a = np.random.rand(128)
+    b = np.random.rand(128)
+    assert np_distances.euclidean_dist(a, b).shape == (1,)
+    assert np_distances.euclidean_dist(a, b) == np_distances.euclidean_dist(b, a)
+    np.testing.assert_array_almost_equal(
+        np_distances.euclidean_dist(a, a), np.zeros((1,))
+    )
+
+    a = np.random.rand(10, 3)
+    b = np.random.rand(5, 3)
+    assert np_distances.euclidean_dist(a, b).shape == (10, 5)
+    assert np_distances.euclidean_dist(b, a).shape == (5, 10)
+    diag_dists = np.diagonal(np_distances.euclidean_dist(b, b))  # self-comparisons
+    np.testing.assert_array_almost_equal(diag_dists, np.zeros((5,)))
+
+    a = np.array([0.0, 2.0, 0.0])
+    b = np.array([0.0, 0.0, 2.0])
+    desired_output_singleton = np.sqrt(np.array([2.0**2.0 + 2.0**2.0]))
+    np.testing.assert_array_almost_equal(
+        np_distances.euclidean_dist(a, b), desired_output_singleton
+    )
+
+    a = np.array([[0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
+    b = np.array([[0.0, 0.0, 2.0], [0.0, 2.0, 0.0]])
+    desired_output_singleton = np.array(
+        [[desired_output_singleton.item(), 0.0], [0.0, desired_output_singleton.item()]]
+    )
+    np.testing.assert_array_almost_equal(
+        np_distances.euclidean_dist(a, b), desired_output_singleton
+    )
+
+
+def test_euclidean_dist_compare():
+    a = torch.rand(128)
+    b = torch.rand(128)
+    torch.testing.assert_close(
+        torch_distances.euclidean_dist(a, b),
+        torch.from_numpy(np_distances.euclidean_dist(a.numpy(), b.numpy())).to(
+            torch.float32
+        ),
+    )
+
+    a = torch.rand(10, 3)
+    b = torch.rand(5, 3)
+    torch.testing.assert_close(
+        torch_distances.euclidean_dist(a, b),
+        torch.from_numpy(np_distances.euclidean_dist(a.numpy(), b.numpy())),
+    )
+
+
+def test_sqeuclidean_dist_torch():
+    a = torch.rand(128)
+    b = torch.rand(128)
+    assert torch_distances.sqeuclidean_dist(a, b).shape == (1,)
+    torch.testing.assert_close(
+        torch_distances.sqeuclidean_dist(a, b),
+        torch_distances.euclidean_dist(a, b) ** 2,
+    )
+
+    a = torch.rand(10, 3)
+    b = torch.rand(5, 3)
+    assert torch_distances.sqeuclidean_dist(a, b).shape == (10, 5)
+    torch.testing.assert_close(
+        torch_distances.sqeuclidean_dist(a, b),
+        torch_distances.euclidean_dist(a, b) ** 2,
+    )
+
+
+def test_sqeuclidea_dist_np():
+    a = np.random.rand(128)
+    b = np.random.rand(128)
+    assert np_distances.sqeuclidean_dist(a, b).shape == (1,)
+    np.testing.assert_array_almost_equal(
+        np_distances.sqeuclidean_dist(a, b), np_distances.euclidean_dist(a, b) ** 2
+    )
+
+    a = np.random.rand(10, 3)
+    b = np.random.rand(5, 3)
+    assert np_distances.sqeuclidean_dist(a, b).shape == (10, 5)
+    np.testing.assert_array_almost_equal(
+        np_distances.sqeuclidean_dist(a, b), np_distances.euclidean_dist(a, b) ** 2
+    )
diff --git a/tests/units/util/test_topk.py b/tests/units/util/test_topk.py
new file mode 100644
index 00000000000..2bd8d3b520a
--- /dev/null
+++ b/tests/units/util/test_topk.py
@@ -0,0 +1,49 @@
+import numpy as np
+import torch
+
+from docarray.utility.helper.numpy import top_k as np_top_k
+from docarray.utility.helper.torch import top_k as torch_top_k
+
+
+def test_topk_torch():
+    a = torch.tensor([1, 4, 2, 7, 4, 9, 2])
+    vals, indices = torch_top_k(a, 3)
+    assert vals.shape == (1, 3)
+    assert indices.shape == (1, 3)
+    assert (vals.squeeze() == torch.tensor([1, 2, 2])).all()
+    assert (indices.squeeze() == torch.tensor([0, 2, 6])).all() or (
+        indices.squeeze() == torch.tensor([0, 6, 2])
+    ).all()
+
+    a = torch.tensor([[1, 4, 2, 7, 4, 9, 2], [11, 6, 2, 7, 3, 10, 4]])
+    vals, indices = torch_top_k(a, 3)
+    assert vals.shape == (2, 3)
+    assert indices.shape == (2, 3)
+    assert (vals[0] == torch.tensor([1, 2, 2])).all()
+    assert (indices[0] == torch.tensor([0, 2, 6])).all() or (
+        indices[0] == torch.tensor([0, 6, 2])
+    ).all()
+    assert (vals[1] == torch.tensor([2, 3, 4])).all()
+    assert (indices[1] == torch.tensor([2, 4, 6])).all()
+
+
+def test_topk_numpy():
+    a = np.array([1, 4, 2, 7, 4, 9, 2])
+    vals, indices = np_top_k(a, 3)
+    assert vals.shape == (1, 3)
+    assert indices.shape == (1, 3)
+    assert (vals.squeeze() == np.array([1, 2, 2])).all()
+    assert (indices.squeeze() == np.array([0, 2, 6])).all() or (
+        indices.squeeze() == np.array([0, 6, 2])
+    ).all()
+
+    a = np.array([[1, 4, 2, 7, 4, 9, 2], [11, 6, 2, 7, 3, 10, 4]])
+    vals, indices = np_top_k(a, 3)
+    assert vals.shape == (2, 3)
+    assert indices.shape == (2, 3)
+    assert (vals[0] == np.array([1, 2, 2])).all()
+    assert (indices[0] == np.array([0, 2, 6])).all() or (
+        indices[0] == np.array([0, 6, 2])
+    ).all()
+    assert (vals[1] == np.array([2, 3, 4])).all()
+    assert (indices[1] == np.array([2, 4, 6])).all()