cond function fails for non-square matrices

[stevengj]

Currently, linalg.jl defines cond(a::AbstractMatrix) = cond(a, 2). This only works for square matrices. The correct definition (as well as being more efficient) is:

cond(a::AbstractMatrix) = let s = svdvals(a); max(s) / min(s); end

And cond(a,2) should also do this, of course.

Note also that cond(a,p) fails with an exception (Singular system) if a is singular, whereas it should really give Inf in that case (or NaN if a is all zeros).

[ViralBShah]

To fix:
cond(A,p) returns Inf for all zeros input except when p = 2.

[toivoh]

What does it mean to define the condition number for a non-square matrix?
The way that I understand it, the condition number for matrix A is the least upper bound on relative error amplification when trying to solve one of the linear equation systems

A*x = y  # recovery subject to left condition number?
u*A = v  # recovery subject to right condition number?

If A is nonsquare, it will be possible to recover at most one of x and u above, given y and v. So it seems that we should really talk about left and right condition numbers in this case?

svdvals currently seems to leave out the trailing zeros beyond min(m,n) for an m x n matrix, so the definition above would give the optimistic answer, which seems a bit dangerous to me. Can we have lcond and rcond for left and right condition numbers (whichever would be which :), and let cond return Inf for any nonsquare matrix?

[ViralBShah]

cond returning Inf for non-square does not seem like a good idea, as that would silently suggest that the matrix is ill-conditioned, when it may actually not be ill-conditioned.

[stevengj]

@ViralBShah, your proposed solution seems wrong to me, in two respects. First, you also have to handle non-zero singular matrices, e.g. cond(diagm([1,0])) should give Inf. Second, the definition of the condition number gives 0*Inf for all-zero matrices, which in IEEE arithmetic should mean that NaN is returned in that case (Matlab gives inf, but that seems wrong to me...it may be a historical artifact of pre-IEEE arithmetic).

@toivoh, the definition of condition number for non-square matrices is totally standard, see e.g. Trefethen's linear-algebra textbook. It is the upper bound on the relative error amplification to when computing the function f(x) = A*x. Given a choice of vector and matrix norms, this is cond(A) = |A| * (sup_x |x|/|Ax|), which is finite as long as A has full column rank (i.e. |Ax| > 0 for |x| > 0).

For the case of the L2 vector norm and the corresponding induced matrix norm, this gives cond(A) = (sup_x |Ax|/|x|) * (sup_x |x|/|Ax|), which equals the ratio of the maximum and minimum singular values.

For any arbitrary vector norm and the corresponding induced matrix norm, if inv(A) exists, one obtains cond(A) = |A| * |inv(A)|, which is the formula you are using. One could also in principle compute cond(A) for non-square A in any arbitrary norm, but in practice it is too complicated to be worth it (Matlab doesn't bother).

[Note that, if inv(A) exists, then cond(A) is also the upper bound on the relative error amplification of f(x) = A\x.]

[ViralBShah]

I currently have this:

julia> cond(diagm([1,0]))
Inf

julia> cond(zeros(5,5))
Inf

I just committed a fix that returns Inf when the result is NaN, which is consistent with Matlab. Wouldn't it be better if the result is Inf for a matrix of all zeros as a usability thing, since you can then check for all ill-conditioned matrices with isinf?

[stevengj]

@ViralBShah, you are begging the question. Being "ill-conditioned" means that the condition number is large. Since the condition number of a matrix of all zeros does not seem to be well-defined, it is not at all clear that it is "ill-conditioned".

For any nonzero matrix A, consider the limit of cond(s_A) as the real number s goes to zero. Depending on A, you can get anything from 1 to infinity from this limit, which means that cond(0_A) is a classic indeterminate form in analysis.

[Except for the case of a 1x1 matrix, in which case the limit is 1. So, arguably, cond(zeros) should return 1 for a 1x1 matrix and NaN in all other cases.]

[Note that f(x)=zeros_x is computed exactly for all x, so the *absolute_ forward error is zero. The problem is that the relative error is 0/0 which is not well defined, even as a limit.]

[ViralBShah]

Ok. I am going to leave this issue open, until we get this right.

[JeffBezanson]

You could also throw an error instead of giving NaN.

[stevengj]

@JeffBezanson, since NaN is expressly provided by IEEE for representing indeterminate forms like 0/0, I'm not sure an Julia exception is called for.

[ViralBShah]

Thinking out aloud, while IEEE provides NaN for this reason and the usage in this case would be correct, is there a reason to use julia exceptions to provide a friendlier interface for julia users?

Certainly worth also checking with @alanedelman on his views.

[toivoh]

@stevengj: Fair enough. So the corollary is that cond(A) = Inf (or possibly NaN) for fat matrices, right?

About taking limits: Since the condition number is an upper bound, I think that the only prudent thing to do would be to take the limit superior, which will indeed be Inf taken over matrices that converge to zero.

I also do not see how anyone would be helped by letting cond(0) = 1 instead of Inf. Is zero not an ill-conditioned operator? If you have an example where this would help things, please share.

[JeffBezanson]

OK, this might be a case where NaN is the best approach. Is it still justifiable to give a domain error for, e.g. sqrt(-1), or does the exact same argument work in that case? I feel like in those cases an error is more useful, since otherwise you can be left wondering where the NaN arose. Is this defensibly a "different case"?