diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index f76d6850bf0..1a240756f3b 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -34,6 +34,16 @@
   [1 1 0]
   ```
 
+* `expectation_value` was added to `qml.math` to calculate the expectation value of a matrix for pure states.
+  [(#4484)](https://github.com/PennyLaneAI/pennylane/pull/4484)
+
+  ```pycon
+  >>> state_vector = [1/np.sqrt(2), 0, 1/np.sqrt(2), 0]
+  >>> operator_matrix = qml.matrix(qml.PauliZ(0), wire_order=[0,1])
+  >>> qml.math.expectation_value(operator_matrix, state_vector)
+  tensor(-2.23711432e-17+0.j, requires_grad=True)
+  ```
+
 * The `default.tensor` device is introduced to perform tensor network simulations of quantum circuits using the `mps` (Matrix Product State) method.
   [(#5699)](https://github.com/PennyLaneAI/pennylane/pull/5699)
 
@@ -470,6 +480,7 @@ Tarun Kumar Allamsetty,
 Guillermo Alonso-Linaje,
 Utkarsh Azad,
 Lillian M. A. Frederiksen,
+Ludmila Botelho,
 Gabriel Bottrill,
 Astral Cai,
 Ahmed Darwish,
diff --git a/pennylane/math/__init__.py b/pennylane/math/__init__.py
index feeba2bf5a5..731979daac9 100644
--- a/pennylane/math/__init__.py
+++ b/pennylane/math/__init__.py
@@ -68,6 +68,7 @@
 from .quantum import (
     cov_matrix,
     dm_from_state_vector,
+    expectation_value,
     marginal_prob,
     mutual_info,
     partial_trace,
@@ -151,6 +152,7 @@ def __getattr__(name):
     "dot",
     "einsum",
     "expand_matrix",
+    "expectation_value",
     "eye",
     "fidelity",
     "fidelity_statevector",
diff --git a/pennylane/math/quantum.py b/pennylane/math/quantum.py
index a1fdeabaf2d..f4748ee7743 100644
--- a/pennylane/math/quantum.py
+++ b/pennylane/math/quantum.py
@@ -472,7 +472,10 @@ def reduce_statevector(state, indices, check_state=False, c_dtype="complex128"):
         [ABC[i + 1] for i in sorted(indices)] + [ABC[num_wires + i + 1] for i in sorted(indices)]
     )
     density_matrix = einsum(
-        f"a{indices1},a{indices2}->a{target}", state, np.conj(state), optimize="greedy"
+        f"a{indices1},a{indices2}->a{target}",
+        state,
+        np.conj(state),
+        optimize="greedy",
     )
 
     # Return the reduced density matrix by using numpy tensor product
@@ -514,7 +517,10 @@ def dm_from_state_vector(state, check_state=False, c_dtype="complex128"):
     """
     num_wires = int(np.log2(np.shape(state)[-1]))
     return reduce_statevector(
-        state, indices=list(range(num_wires)), check_state=check_state, c_dtype=c_dtype
+        state,
+        indices=list(range(num_wires)),
+        check_state=check_state,
+        c_dtype=c_dtype,
     )
 
 
@@ -663,7 +669,14 @@ def _compute_vn_entropy(density_matrix, base=None):
 
 
 # pylint: disable=too-many-arguments
-def mutual_info(state, indices0, indices1, base=None, check_state=False, c_dtype="complex128"):
+def mutual_info(
+    state,
+    indices0,
+    indices1,
+    base=None,
+    check_state=False,
+    c_dtype="complex128",
+):
     r"""Compute the mutual information between two subsystems given a state:
 
     .. math::
@@ -719,29 +732,139 @@ def mutual_info(state, indices0, indices1, base=None, check_state=False, c_dtype
         raise ValueError("Subsystems for computing mutual information must not overlap.")
 
     return _compute_mutual_info(
-        state, indices0, indices1, base=base, check_state=check_state, c_dtype=c_dtype
+        state,
+        indices0,
+        indices1,
+        base=base,
+        check_state=check_state,
+        c_dtype=c_dtype,
     )
 
 
 # pylint: disable=too-many-arguments
 def _compute_mutual_info(
-    state, indices0, indices1, base=None, check_state=False, c_dtype="complex128"
+    state,
+    indices0,
+    indices1,
+    base=None,
+    check_state=False,
+    c_dtype="complex128",
 ):
     """Compute the mutual information between the subsystems."""
     all_indices = sorted([*indices0, *indices1])
     vn_entropy_1 = vn_entropy(
-        state, indices=indices0, base=base, check_state=check_state, c_dtype=c_dtype
+        state,
+        indices=indices0,
+        base=base,
+        check_state=check_state,
+        c_dtype=c_dtype,
     )
     vn_entropy_2 = vn_entropy(
-        state, indices=indices1, base=base, check_state=check_state, c_dtype=c_dtype
+        state,
+        indices=indices1,
+        base=base,
+        check_state=check_state,
+        c_dtype=c_dtype,
     )
     vn_entropy_12 = vn_entropy(
-        state, indices=all_indices, base=base, check_state=check_state, c_dtype=c_dtype
+        state,
+        indices=all_indices,
+        base=base,
+        check_state=check_state,
+        c_dtype=c_dtype,
     )
 
     return vn_entropy_1 + vn_entropy_2 - vn_entropy_12
 
 
+def _check_hermitian_operator(operators):
+    """Check the shape, and if the matrix is hermitian."""
+    dim = operators.shape[-1]
+
+    if (
+        len(operators.shape) not in (2, 3)
+        or operators.shape[-2] != dim
+        or not np.log2(dim).is_integer()
+    ):
+        raise ValueError(
+            "Operator matrix must be of shape (2**wires,2**wires) "
+            "or (batch_dim, 2**wires, 2**wires)."
+        )
+
+    if len(operators.shape) == 2:
+        operators = qml.math.stack([operators])
+
+    if not is_abstract(operators):
+        for ops in operators:
+            conj_trans = np.transpose(np.conj(ops))
+            if not allclose(ops, conj_trans):
+                raise ValueError("The matrix is not Hermitian.")
+
+
+def expectation_value(
+    operator_matrix, state_vector, check_state=False, check_operator=False, c_dtype="complex128"
+):
+    r"""Compute the expectation value of an operator with respect to a pure state.
+
+    The expectation value is the probabilistic expected result of an experiment.
+    Given a pure state, i.e., a state which can be represented as a single
+    vector :math:`\ket{\psi}` in the Hilbert space, the expectation value of an
+    operator :math:`A` can computed as
+
+    .. math::
+        \langle A \rangle_\psi = \bra{\psi} A \ket{\psi}
+
+
+    Args:
+        operator_matrix (tensor_like): operator matrix with shape ``(2**N, 2**N)`` or ``(batch_dim, 2**N, 2**N)``.
+        state_vector (tensor_like): state vector with shape ``(2**N)`` or ``(batch_dim, 2**N)``.
+        check_state (bool): if True, the function will check the validity of the state vector
+            via its shape and the norm.
+        check_operator (bool): if True, the function will check the validity of the operator
+            via its shape and whether it is hermitian.
+        c_dtype (str): complex floating point precision type.
+
+    Returns:
+        float: Expectation value of the operator for the state vector.
+
+    **Example**
+
+    The expectation value for any operator can obtained by passing their matrix representation as an argument.
+    For example, for a 2 qubit state, we can compute the expectation value of the operator :math:`Z \otimes I` as
+
+
+    >>> state_vector = [1/np.sqrt(2), 0, 1/np.sqrt(2), 0]
+    >>> operator_matrix = qml.matrix(qml.PauliZ(0), wire_order=[0,1])
+    >>> qml.math.expectation_value(operator_matrix, state_vector)
+    tensor(-2.23711432e-17+0.j, requires_grad=True)
+
+    .. seealso:: :func:`pennylane.math.fidelity`
+
+    """
+    state_vector = cast(state_vector, dtype=c_dtype)
+    operator_matrix = cast(operator_matrix, dtype=c_dtype)
+
+    if check_state:
+        _check_state_vector(state_vector)
+
+    if check_operator:
+        _check_hermitian_operator(operator_matrix)
+
+    if qml.math.shape(operator_matrix)[-1] != qml.math.shape(state_vector)[-1]:
+        raise qml.QuantumFunctionError(
+            "The operator and the state vector must have the same number of wires."
+        )
+
+    # The overlap <psi|A|psi>
+    expval = qml.math.einsum(
+        "...i,...i->...",
+        qml.math.conj(state_vector),
+        qml.math.einsum("...ji,...i->...j", operator_matrix, state_vector, optimize="greedy"),
+        optimize="greedy",
+    )
+    return expval
+
+
 # pylint: disable=too-many-arguments
 def vn_entanglement_entropy(
     state, indices0, indices1, base=None, check_state=False, c_dtype="complex128"
@@ -884,7 +1007,10 @@ def _compute_relative_entropy(rho, sigma, base=None):
     # the matrix of inner products between eigenvectors of rho and eigenvectors
     # of sigma; this is a doubly stochastic matrix
     rel = qml.math.einsum(
-        f"{indices_rho},{indices_sig}->{target}", np.conj(u_rho), u_sig, optimize="greedy"
+        f"{indices_rho},{indices_sig}->{target}",
+        np.conj(u_rho),
+        u_sig,
+        optimize="greedy",
     )
     rel = np.abs(rel) ** 2
 
diff --git a/tests/math/test_expectation_value.py b/tests/math/test_expectation_value.py
new file mode 100644
index 00000000000..ff07406e68a
--- /dev/null
+++ b/tests/math/test_expectation_value.py
@@ -0,0 +1,200 @@
+# Copyright 2024 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Unit tests for differentiable matrix-vector expectation values.
+"""
+
+import numpy as onp
+import pytest
+
+import pennylane as qml
+from pennylane import numpy as np
+
+pytestmark = pytest.mark.all_interfaces
+
+tf = pytest.importorskip("tensorflow", minversion="2.1")
+torch = pytest.importorskip("torch")
+jax = pytest.importorskip("jax")
+jnp = pytest.importorskip("jax.numpy")
+
+
+class TestExpectationValueMath:
+    """Tests for Expectation value of a operator for a state vector."""
+
+    ops_vs_vecstates = [
+        ([[1, 0], [0, 0]], [1, 0], 1),
+        ([[0, 1], [1, 0]], [0, 1], 0),
+        ([[0.5, 0.5], [0.5, 0.5]], [1, 1] / np.sqrt(2), 1),
+        (
+            [[0.40975111, 0.40751457], [0.40751457, 0.59024889]],
+            [0.8660254, 0.5],
+            0.8077935208042251,
+        ),
+        (
+            [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, -1, 0], [0, 0, 0, -1]],
+            [1, 0, 1, 0] / np.sqrt(2),
+            0,
+        ),
+    ]
+
+    array_funcs = [
+        lambda x: x,
+        onp.array,
+        np.array,
+        jnp.array,
+        torch.tensor,
+        tf.Variable,
+        tf.constant,
+    ]
+
+    @pytest.mark.parametrize("operator_and_states", ops_vs_vecstates)
+    @pytest.mark.parametrize("func", array_funcs)
+    def test_mat_expectation_value(self, operator_and_states, func):
+        """Test the expectation value of a operator for a vector state."""
+        ops, state_vectors, expected = operator_and_states
+        ops = func(ops)
+        state_vectors = func(state_vectors)
+        overlap = qml.math.expectation_value(ops, state_vectors)
+        assert qml.math.allclose(expected, overlap)
+
+    state_wrong_amp = [
+        ([[1, 0], [0, 1]], [0.5, 0]),
+        ([[1, 0], [0, 1]], [26, 70]),
+        (
+            [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, -1, 0], [0, 0, 0, -1]],
+            [1, 2, 0, 0],
+        ),
+        (
+            [[0, 1, 0, 0], [1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]],
+            [1, 1, 1, 1],
+        ),
+    ]
+
+    @pytest.mark.parametrize("ops,state_vectors", state_wrong_amp)
+    def test_state_vector_wrong_amplitudes(self, ops, state_vectors):
+        """Test that a message is raised when a state does not have right norm"""
+        with pytest.raises(ValueError, match="Sum of amplitudes-squared does not equal one."):
+            qml.math.expectation_value(ops, state_vectors, check_state=True)
+
+    operator_wrong_shape = [
+        ([[1, 0, 0], [0, 0, 1]], [0, 1]),
+        (
+            [[1, 0, 0, 0, 0], [0, 1, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, 0, 1, 0], [0, 0, 0, 0, 1]],
+            [1, 0, 0, 0],
+        ),
+        ([[0, 1, 0], [1, 0, 0], [0, 0, 0], [0, 0, 1]], [1, 0]),
+    ]
+
+    @pytest.mark.parametrize("ops,state_vectors", operator_wrong_shape)
+    def test_operator_wrong_shape(self, ops, state_vectors):
+        """Test that a message is raised when the state does not have the right shape."""
+        with pytest.raises(ValueError, match="Operator matrix must be of shape"):
+            qml.math.expectation_value(ops, state_vectors, check_operator=True)
+
+    operator_non_hermitian = [
+        ([[1, 1], [0, 0]], [0, 1]),
+        (
+            [[0, 0, 0, 1], [0, 1, 0, 1], [1, 0, 1, 0], [1, 0, 0, 1]],
+            [1, 0, 0, 0],
+        ),
+    ]
+
+    @pytest.mark.parametrize("ops,state_vectors", operator_non_hermitian)
+    def test_operator_non_hermitian(self, ops, state_vectors):
+        """Test that a message is raised when the state does not have the right shape."""
+        with pytest.raises(ValueError, match="The matrix is not Hermitian"):
+            qml.math.expectation_value(ops, state_vectors, check_operator=True)
+
+    state_wrong_shape = [
+        ([[1, 0], [0, 1]], [0, 1, 0]),
+        ([[1, 0], [0, 1]], [0, 0, 0, 0, 1]),
+        (
+            [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, -1, 0], [0, 0, 0, -1]],
+            [1, 0, 0, 0, 0],
+        ),
+        ([[0, 1, 0, 0], [1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]], [1, 0, 0]),
+    ]
+
+    @pytest.mark.parametrize("ops,state_vectors", state_wrong_shape)
+    def test_state_vector_wrong_shape(self, ops, state_vectors):
+        """Test that a message is raised when the state does not have the right shape."""
+        with pytest.raises(ValueError, match="State vector must be of shape"):
+            qml.math.expectation_value(ops, state_vectors, check_state=True)
+
+    def test_same_number_wires_dm(self):
+        """Test that the two states must act on the same number of wires"""
+        ops = np.diag([0, 1, 0, 0])
+        state_vectors = [1, 0]
+        with pytest.raises(
+            qml.QuantumFunctionError,
+            match="The operator and the state vector must have the same number of wires.",
+        ):
+            qml.math.expectation_value(ops, state_vectors, check_state=True, check_operator=True)
+
+    @pytest.mark.parametrize("func", array_funcs)
+    def test_broadcast_op_sv(self, func):
+        """Test simultaneous broadcasting of operators and state vectors works."""
+        ops = qml.math.stack(
+            (
+                func([[1, 0], [0, 0]]),
+                func([[0.5, 0.5], [0.5, 0.5]]),
+                func([[0, 0], [0, 1]]),
+                func([[0.40975111, 0.40751457], [0.40751457, 0.59024889]]),
+            )
+        )
+        state_vectors = qml.math.stack(
+            [
+                func([0, 1]),
+                func([1, 0]),
+                func([1, 1] / np.sqrt(2)),
+                func([0.8660254, 0.5]),
+            ]
+        )
+        expected = [0, 0.5, 0.5, 0.8077935208042251]
+
+        overlap = qml.math.expectation_value(ops, state_vectors)
+        assert qml.math.allclose(overlap, expected)
+
+    @pytest.mark.parametrize("func", array_funcs)
+    def test_broadcast_op_unbatched(self, func):
+        """Test broadcasting works for expectation values when the operators input is unbatched"""
+        ops = func([[1, 0], [0, 0]])
+        state_vectors = qml.math.stack(
+            [
+                func([0, 1]),
+                func([1, 0]),
+                func([1, 1] / np.sqrt(2)),
+                func([0.8660254, 0.5]),
+            ]
+        )
+        expected = [0, 1, 0.5, 0.7499999934451599]
+
+        overlap = qml.math.expectation_value(ops, state_vectors)
+        assert qml.math.allclose(overlap, expected)
+
+    @pytest.mark.parametrize("func", array_funcs)
+    def test_broadcast_sv_unbatched(self, func):
+        """Test broadcasting works for expectation values when the state vector input is unbatched"""
+        ops = qml.math.stack(
+            (
+                func([[1, 0], [0, 0]]),
+                func([[0.5, 0.5], [0.5, 0.5]]),
+                func([[0, 0], [0, 1]]),
+                func([[0.40975111, 0.40751457], [0.40751457, 0.59024889]]),
+            )
+        )
+        state_vectors = func([1, 0])
+        expected = [1, 0.5, 0, 0.40975111]
+
+        overlap = qml.math.expectation_value(ops, state_vectors)
+        assert qml.math.allclose(overlap, expected)