COND for homogeneous equations?

Normally, when we want to compute error tendency in the solution of a system of heterogeneous linear equations,

A*x=b

we compute the condition number cond(A) or rcond(A). I am curious to know if there is an analogous way to assess homogeneous equations,

A*x=0

In computer vision, it is common for a desired solution to be the unique (ideally) null vector of such a system. But how do we measure how robustly one-dimensional (or more) the null space of a matrix is? It's occurred to me that you could do it by looking at the ratio between the smallest two singular values,

S=svd(A);
metric = S(end-1)/S(end);

but I was curious if there was a more established metric, and a corresponding MATLAB builtin that computes it.

EDIT: And additionally, how would we relate this metric to accuracy of the solution? If A is rank n-1, can we put bounds on the error in computing the null space basis vector, x? And what about rank n-m, m>1?

[John D'Errico]

I think you should write a tool for this purpose. Call it something like singularityReport perhaps. It could generate various measures about the matrix, terribly useful for someone who does not really understand the SVD.

How close is the matrix to singularity?

If singular, what is the dimension of the null space, and what further perturbation of the matrix would be required to increase the dimension of the null space?

If not singular, how will solving a linear system using this matrix amplify any noise in the problem?

You could make it so the tool returns a set of measures in a struct, but provide a verbosity parameter to allow the tool to explain those measures in more detail. This would make it a useful tool to someone fairly new to numerical linear algebra who wants to solve a linear system.

Answer 1

[John D'Errico]

The question is, I know you are singular, but just HOW singular are you? ;-) It seems to me there are several questions of interest.

- Is A numerically singular at all? Cond gives us a measure of that.

- What is the rank of the null space? Null gives us an opinion about that, but not a metric.

A measure of how close to entrance into the null space is to consider the tolerance that null would use. That depends on the largest singular value. Any singular value that lies below some fraction of the largest singular value, results in the corresponding singular vector getting dumped into the null space.

tol = max(m,n) * eps(max(S));

So I would look at the smallest singular value that lies above that tolerance. Then compute something like this:

min(S(S > tol))/tol

So how close is that value to the tolerance for entrance into the null space?

Note that this applies to a matrix that is numerically full rank too. For example...

A = vander(1:10);
[m,n] = size(A);

tol = max(m,n) * eps(max(S))
tol =
    2.38418579101562e-06

min(S(S > tol))/tol
ans =
        216.509958943425

So A is close to singular, but not there yet.

But if I push things a bit harder...

A = vander(1:15);
[m,n] = size(A);
S = svd(A)
S =
      3.15839869049552e+16
          66664300967337.7
          322942655858.374
          2723747542.51185
          36045801.6655343
           713911.28106567
          20777.0142721889
          891.531934617208
          58.0952625517871
           6.2867119286002
          1.42076230762098
          0.36061079623119
        0.0322030740265822
       0.00106225764978576
      1.22304251430976e-05

tol = max(m,n) * eps(max(S))
tol =
    60

So it looks like 7 singular values fall below tol. Easy to compute the rank of the null space.

sum(S<tol)
ans =
     7

And the next one is not that far behind.

min(S(S > tol))/tol
ans =
        14.8588655769535

Or, perhaps report it in terms of the log10 of that, which gives the answer in terms of how many powers of 10 does it lie above the tolerance.

log10(min(S(S > tol))/tol)
ans =
          1.17198565380708

Another one. I recall that magic squares (as produced by magic, for order > 2) are singular for even order.

A = magic(4);
[m,n] = size(A);
S = svd(A)
S =
                        34
          17.8885438199983
          4.47213595499958
      4.17280694870505e-16

tol = max(m,n) * eps(max(S))
tol =
       2.8421709430404e-14

min(S(S > tol))/tol
ans =
   1.5735e+14

So the rank of the null space was 1, but the next singular value is nowhere close to that tolerance.

So I'd think you need to report two measures.

1. How many singular values fall below tol?

2. How close is the next one to tolerance failure?

[John D'Errico]

I'd be interested if you manage to come up with some tool that would investigate this problem.

Answer 2

[Christine Tobler]

That's an interesting question, I'm not sure what the right way to do this is. In a first step, the matrix condition number describes how much an error in the right-hand side b can be amplified in the solution vector x. So if there is mathematically exact x_exact and b_exact with A*x_exact = b_exact, then the best numerical solution x will satisfy

 || A*x - b_exact|| / ||b_exact|| <= cond(A) * ||x - x_exact|| / ||x_exact||.

When b is always zero, this relative error bound does not make sense. Also, in this bound we are assuming that A can be represented exactly, which may not make sense for the homogeneous equation. So strictly speaking, the first step would be to figure out what kind of error in the solution x we want to estimate.

That said, I think the ratio of the last and second-to-last singular value you propose is an intuitive measure that should work great. I'd just propose to use them the other way round, S(end) / S(end-1), so that a matrix with exactly rank n-1 will show a value of zero.

Answer 3

[Matt J]

Thank you John and Christine for your responses. It's good input. Now that I think of it, though, I think Christine was right about my question not being fully formulated. We really want a metric that, analogous to cond(A) will give accuracy bounds on the null vector of a rank n-1 matrix (and maybe something more general for rank n-m).

[Christine Tobler]

I think John is right that there need to be two measures to check here, and that we should look at more of the singular values. I would just describe the two conditions a bit differently:

1. Does a solution to A*x = 0 exist? (that is, is the matrix singular?)-> we can measure this using cond(A), i.e., s(end)/s(1). This answers "How small a (relative) perturbation of A could make it exactly singular?"
2. Is the solution to A*x = 0 unique? Here's where the second smallest singular value comes in, which we could measure in two ways:

a) s(end-1)/s(1): "How small a (relative) perturbation of A would cause it to be rank n-2?"

b) s(end-1)/s(end): "If the solution x = V(:,end) gives us norm(A*x), how much smaller will norm(A*x2) be using x2 = V(:, end-1)?"

For square matrices, existence and uniqueness of the solution are coupled, so there it's enough to look at the condition number alone.