I'm going to suggest a new function here. I don't feel strongly about broadcasting either way, but there have been requests for the covariance matrix of the fitted parameters and I often end up computing that myself as it is more informative and allows for generating error bars on curve fits and such. Unfortunately, I can't see if DGELSD provides sufficient information to do that, the documentation simply says that the input design matrix is destroyed. DGELSS looks like an option but I can't find it in our fallback routines. I suppose we could implement it for ourselves. Multiple return options would be required as computing the covariance could take significant time.

That sounds like an orthogonal suggestion here - the only point of this change is to enable broadcasting, since lstsq is one of the few remaining unbroadcast functions in np.linalg.

Unfortunately, the meaning of resids has to be tweaked every so slightly to make those semantics possible.

It seemed to me the main change is that resids previously had (in some sense) a keepdims in there, while now it does not? I.e. you do not add a dimension of size 1 at the end anymore.
For all practical purposes that probably does not matter much hopefully and I guess without it maybe it generalizes slight less well.
Although, I am not immediately sure, is that change for the simplest call actually necessary (I bet you guys discussed it before)?

It seemed to me the main change is that resids previously had (in some sense) a keepdims in there, while now it does not? I.e. you do not add a dimension of size 1 at the end anymore.

It's a little worse than that. The old behavior was to return any of:

• A result of shape (0,) for an underconstrained solution (for any input shape)
• A result of shape (1,) for a well-constrained solution with input shape (M,)
• A result of shape (K,) for a well-constrained solution with input shape (M, K)

The new behavior is

• A result of shape (..., N,) for input shape (..., N, M)
• A result of shape () for input shape (M,)

So there are two incompatible changes here, not just one.

I suppose I was guessing that the majority of use cases is likely currently either in the (K,) category (unchanged) or in the (1,) category which changed to (). The old underconstrained one seems strange, but of course it is not impossible someone actually checks len(constr) or so to see whether it was underconstrained... So in that sense, breaking the (1,) case may be better anyway, because at least it should break loudly when it breaks.

So the only option would be keeping the (0,) and (1,) shapes even though they are ugly... It is super hard to say how disruptive that change could be. It seems a question of boldness, this is the right thing if we are bold enough to risk breaking someone.
With release notes, I think I may be happy if we expect that almost all code that actually could break by this, should break loudly.

That sounds like an orthogonal suggestion here

Yes, I thought it might be more useful in the long run :) Computing the sum of squared residuals also has a noticeable impact on performance, which has also been the subject of complaints.

But back to the case at hand with some explanation for other reviewers.

The problem with the stacked case it that the dimensions of the returned residual array must be determined by M, N, and K and cannot depend on runtime results like matrix rank. DEGLS will always compute residuals if M > N, regardless of if the matrix is rank N or not, there just might be some missing residuals if it isn't rank N, so I would just return the results from those residuals in that case, although that becomes problematical when there are few residuals. When the matrix has rank >= M, zeros seems the appropriate return if we are not going to return an empty array. In all other cases NaN seems appropriate, meaning "unknown".

Unfortunately, we fixed the function rather than the documentation in #9986, so there is some question as to whether we should retain the NaN array return when the matrix rank is < N regardless.

The upshot is that we need to make the change to deal with stacked arrays. It is enough of a corner case that I think we can do that, having zero/NaN arrays should deal with the case of downstream users being caught without being aware that something has changed.

In short pseudo code

if M > N:
    # There are residuals returned, use them
    # Note that this result should be divided by (M - N)
    # to get unbiased estimate of variance if the model is any good.
    return sum_of_squares
elseif rank(a) >= M:
   # No residuals returned, but they are zero
    return zeros
else:
   #  No residuals returned, but they exist
   return NaNs

And an array always returned for the residuals.

@IvanYashchuk since you've just spent a bunch of time dealing with this residuals shape issue for PyTorch, would you mind having a look at this PR?