def mat_mult_along(x, y, axis_x, axis_y=None):
    """Multiply matrices spanning arbitrary axes.

    Returns (sum over n)(x_n*y_n) where n runs of axis_x in x and axis_y in y.

    Args:
        x (array): left matrix
        y (array): right matrix
        axis_x (int): axis of x along which to sum
        axis_y (int): axis of y along which to sum. Defaults to axis_x+1

    Returns:
        matrix product
    """
    if axis_y is None:
        axis_y = axis_x + 1
    assert x.shape[axis_y] == y.shape[axis_x]
    return sum(x.take([n], axis_y)*y.take([n], axis_x) for n in range(x.shape[axis_y]))