> > > But, backward compatibility aside, could we have ONLY Scalars?
> > Well, it is hard to write functions that work on N-Dimensions
> > (where N
> > can be 0), if the 0-D array does not exist. 

So as a (silly) example, the following does not generalize to 0d, even
though it should:

def weird_normalize_by_trace_inplace(stacked_matrices)
    """Devides matrices by their trace but retains sign
    (works in-place, and thus e.g. not for integer arrays)

    Parameters
    ----------
    stacked_matrices : (..., N, M) ndarray
    """
    assert stacked_matrices.shape[-1] == stacked_matrices.shape[-2]

    trace = np.trace(stacked_matrices, axis1=-2, axis2=-1)
    trace[trace < 0] *= -1
    stacked_matrices /= trace

Sure that function does not make sense and you could rewrite it, but
the fact is that in that function you want to conditionally modify
trace in-place, but trace can be 0d and the "conditional" modification
breaks down.