How does the default tolerance of "rank" work in MATLAB R2018a and newer?

I am checking the function "rank" in MATLAB, and notice it has a tolerance. I notice that the value of this tolerance is obtained as:

"The rank of a matrix A is computed as the number of singular values that are larger than a tolerance. By default, the tolerance is max(size(A))*eps(norm(A)). However, you can specify a different tolerance with the command rank(A,tol)."

Are there any references available explaining this value?

Answer 1

[MathWorks Support Team]

In general, the rank decision is based on how many singular values are not 0. In numerics, however, one rarely ever sees real zeros. Hence, we need to truncate by some 'tol'.

This 'tol' now depends on norm(A) which is, in fact, just the largest singular value.

So eps(norm(A)), as described in help of EPS, "is the positive distance from ABS(norm(A)) to the next larger in magnitude floating point number of the same precision as X."

In other words, based on the largest singular value, which is norm(A), we cannot be more accurate as values smaller than this cannot be captured as they are smaller than the distance to the next floating-point number.

Additionally, we need to scale this value in some way. And there are various ways. A very simple way is just by the number of cols/rows in the input, which is represented by max(size(A)).

This can be interpreted as:

Each row/col we process in the algorithm is allowed to introduce a round-off error of eps(norm(A)) without being numerically incorrect, hence, for very large matrices, we still allow some significant tolerance.

 

Here are some examples:

We create a block matrix with Identity (eye) and a small value p on the diagonal. The singular values are all 1 and p. Depending on p, the rank changes. We print p, rank r, and the tol used in rank:

>> p=eps, A=blkdiag(eye(10),p);
r=rank(A);
tol=max(size(A))*eps(norm(A))
p =
   2.2204e-16
r =
    10
tol =
   2.4425e-15


You see the tol is much larger then p, hence, we have rank 10 and not 11. If we now make p larger than tol, we get rank 11:
>> p=12*eps, A=blkdiag(eye(10),p); r=rank(A), tol=max(size(A))*eps(norm(A))
p =
   2.6645e-15
r =
11
tol =
   2.4425e-15


However, if we now make the matrix 101x101, we must allow more round-off errors and the same p is not large enough:
>> p=12*eps, A=blkdiag(eye(100),p); r=rank(A), tol=max(size(A))*eps(norm(A))
p =
   2.6645e-15
r =
   100
tol =
   2.2427e-14


As we have 10x as many cols/rows, let's make p 10x larger and we get rank 101:

>> p=120*eps, A=blkdiag(eye(100),p); rank(A), max(size(A))*eps(norm(A))
p =
   2.6645e-14
ans =
101
ans =
   2.2427e-14

The reference of the algorithm that is used to compute the tolerance for "rank" is in Linpack user guide section Chapter 11, section 1.

Because the default tolerance is in a way 'randomly' chosen, and some users might need their own 'tol', "rank" gives the option to do so. This is just one way of scaling. Some applications might need better tuning, and this can be done using the 'tol' input.