Let cond of an empty matrix return zero

[martinholters]

julia> cond(ones(0,0), 1)
0.0 # no change

julia> cond(ones(0,0), 2)
ERROR: ArgumentError: reducing over an empty collection is not allowed # master
0.0 # this PR

julia> cond(ones(0,0), Inf)
0.0 # no change

This is also consistent with opnorm:

julia> opnorm(ones(0,0), 2) * opnorm(inv(ones(0,0)), 2)
0.0

Fixes JuliaLang/LinearAlgebra.jl#778.

[KlausC]

I propose to do isempty(A) && return zero(eltype(A)) just before the call to svdvals.
That would cover also the cases, when svdvals does not succeed (for example if the element type is BigFloat).
Also add test cases for cond(ones(BigFloat, 0, 0), p) for p in (1, 2, Inf)
and for cond(ones(10, 0), 2) and cond(ones(0,3), 2)

[antoine-levitt]

What's the rationale for this? Condition numbers are always greater than 1 so returning 0 feels very weird here.

[KlausC]

What's the rationale for this?

I think there are at least 3 good reasons.

1. the common definition for all square matrices is cond(A) := opnorm(A) * opnorm(A^-1). (assuming A^-1 == A for empty A)
2. cond(A, 1) == 0 == cond(A, Inf) already for empty A
3. it is done the same way in MATLAB

[martinholters]

The main motivation for me here was consistency, see also the original description. In fact, I had been skeptical about this giving zero myself.

But there is a point to be made why it indeed should be zero: The docstring of cond says: "computed using the operator p-norm". For square matrices, I take that to mean cond(A, p) = opnorm(A, p) * opnorm(inv(A), p). And opnorm for an empty matrix gives zero, naturally leading to a condition number of zero. I think the operator norm for an empty matrix should indeed be zero: The smallest non-negative c such that norm(A*x) <= c*norm(x) for any x is certainly zero if both A and x are empty. And the empty matrix should arguable also be its own inverse, leading to a condition number of zero.

But there is also an argument when using the definition cond(A) = maximum(svdvals(A)) / minimum(svdvals(A)) (for p==2): The singular values are non-negative, so I could safely add init arguments as follows: cond(A) = maximum(svdvals(A); init=0.0) / minimum(svdvals(A); init=Inf). This, too, would yield zero for empty A.

[antoine-levitt]

OK, the part that I find questionable is that the inverse of an empty matrix should be well-defined. But if that is the way it is then indeed this definition is more consistent.

[oscardssmith]

It seems well defined to me. The empty matrix is a diagonal matrix with 1s on it's diagonals and 0s elsewhere. Thus, the inverse of an empty matrix is itself.

[antoine-levitt]

The empty matrix is a diagonal matrix with 1s on it's diagonals and 0s elsewhere

That makes no sense to me: where are you taking this from? As far as I can see, the semantics julia defines on empty arrays are those of the "neutral" element for reductions, eg sum, prod, det, tr, etc. I would imagine that eg pinv of the empty array is well-defined, but inv would be an error.

[martinholters]

Certainly empty_square_matrix * empty_square_matrix * x == x for any x (of appropriate size, i.e. empty), so empty_square_matrix appears to be its own inverse.

[antoine-levitt]

Well OK, the algebra of empty matrices is well-defined as the single-element algebra, where every operation maps to the single empty matrix, so from that point of view it makes sense (although I would imagine in practice it would create more bugs than it simplifies code).

Cond is different though, and it doesn't seem to me that there is an unambiguous way of defining it. Eg if you take the definition https://wikimedia.org/api/rest_v1/media/math/render/svg/613bce563088d5e506597b022fb9b09461e39199 seriously, you get -inf. It could also plausibly be 1 (the "best case") or +inf (the "worst case"). So it would perhaps be wiser not to define it.

[martinholters]

Cond is different though, and it doesn't seem to me that there is an unambiguous way of defining it.

Yeah, but the choice made for cond has been documented (if maybe somewhat vague) by referring to the operator norm.

Would you argue that we should also change the behavior for p != 2 here? To also error perhaps? I'm not completely opposed to that, but obviously, that would be a breaking change.

[antoine-levitt]

Well naively I would err on the side of caution and error (even if p != 2). But if there's precedent for this (both p != 2) and the matlab convention), maybe it makes sense to make this change. It's unlikely to matter in practice, hopefully.

[DilumAluthge]

We have moved the LinearAlgebra stdlib to an external repo: https://github.com/JuliaLang/LinearAlgebra.jl

@martinholters If you think that this PR is still relevant, please open a new PR on the LinearAlgebra.jl repo.