def MatrixSolveLs(a, rhs, l2_reg):
    """
    Matrix least-squares solve op.
    """
    r = np.empty(rhs.shape).astype(a.dtype)
    for coord in np.ndindex(a.shape[:-2]):
        pos = coord + (Ellipsis,)
        r[pos] = np.linalg.lstsq(a[pos], rhs[pos])[0]

    return r,