def bellman_operator(self, v, Tv=None, sigma=None):
        """
        The Bellman operator, which computes and returns the updated
        value function `Tv` for a value function `v`.

        Parameters
        ----------
        v : array_like(float, ndim=1)
            Value function vector, of length n.

        Tv : ndarray(float, ndim=1), optional(default=None)
            Optional output array for Tv.

        sigma : ndarray(int, ndim=1), optional(default=None)
            If not None, the v-greedy policy vector is stored in this
            array. Must be of length n.

        Returns
        -------
        Tv : ndarray(float, ndim=1)
            Updated value function vector, of length n.

        """
        vals = self.R + self.beta * self.Q.dot(v)  # Shape: (L,) or (n, m)

        if Tv is None:
            Tv = np.empty(self.num_states)
        self.s_wise_max(vals, out=Tv, out_argmax=sigma)
        return Tv