<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def apply_M(self, ax, ay):
 """Linear operator that converts ax, ay to abcd. """

        jac = numpy.array(
            [[self.dx.dot(ax), self.dy.dot(ax)], [self.dx.dot(ay), self.dy.dot(ay)]]
        )

        # jacs and J are of shape (2, 2, k). M must be of the same shape and
        # contain the result of the k 2x2 dot products. Perhaps there's a
        # dot() for this.
        M = numpy.einsum("ijl,jkl->ikl", jac, self.J)
        # M = numpy.array([
        #     [
        #         jac[0][0]*self.J[0][0] + jac[0][1]*self.J[1][0],
        #         jac[0][0]*self.J[0][1] + jac[0][1]*self.J[1][1],
        #     ],
        #     [
        #         jac[1][0]*self.J[0][0] + jac[1][1]*self.J[1][0],
        #         jac[1][0]*self.J[0][1] + jac[1][1]*self.J[1][1],
        #     ],
        #     ])

        # One could use
        #
        #     M = numpy.moveaxis(M, -1, 0)
        #     _, sigma, _ = numpy.linalg.svd(M)
        #
        # but computing the singular values explicitly via
        # <https://scicomp.stackexchange.com/a/14103/3980> is faster and more
        # explicit.
        a = (M[0, 0] + M[1, 1]) / 2
        b = (M[0, 0] - M[1, 1]) / 2
        c = (M[1, 0] + M[0, 1]) / 2
        d = (M[1, 0] - M[0, 1]) / 2

        return a, b, c, d