Default rank revealing tolerance

In few MATLAB functions such as RANK, PINV, ORTH there is a parameter TOL.

MATLAB uses the default value as

% tol = max(size(A)) * eps(norm(A))

where A the matrix under study.

I have two quetions:

First question: what is the theoretical basis of such formula? Especially why the max(size(A)) factor?

In my simple test here, it seems the error or the singular value grows like a square root of the dimension m. Here is the script to create a matrix A of size m x 4 of rank 3. I use SVD to compute the singulars values, s(4) ideally should be 0, s1(1:3) about unity, but I see experimentally s(4) grow like sqrt(m) as showed in this script:

m0 = 1e5;
s4 = [];
m = [];
while m0 < 1e8
    [Q, ~] =qr(rand(m0,3),0);
    X = rand(3,4);
    A = Q * X;
    [~,s,~] = svd(A,0,'vector');
    s4(end+1) = s(4);
    m(end+1) = m0;
    m0 = ceil(m0*1.1);
end

P = polyfit(log(m),log(s4),1);
p = P(1);
c=exp(P(2));

figure
loglog(m, s4);
hold on
h = loglog(m, c*m.^p);
xlabel('m')
ylabel('s(4)')
legend(h, sprintf('%g*m^{%g}', c, p))

PS: There is a floor noise about 1e-16 for m < 1e5, that could affect the fit, so I take the m value above 1e5 up to 1e8 .

Intuitively the square-root relationship seems reasonable to me since we could consider truncation errors in each component are somewhat random and independent.

Second question. The function lsqminnorm also have a default tol but is it NOT specified beside telling it computed from QR decomposition. It describes how the rank k is used to approximate the solution, but I don't see how the tol value is defined to estimate in turn the rank k. I would expect something like

% tol = eps(abs(R(1,1))) * max(size(A));

where R is the triangular matrix returned by qr, meaning abs(R(1,1)) = max(vecnorm(A)).

Can someone (TMW staff) shed a light and give the method to determine the default tol value implemented within lsqminnorm?

[Paul]

TMW Support team addresses your first question in this question

[Bruno Luong]

@Paul Thanks. The reference is LINPACK userg guide. I can't access to the latest version, but some kind of error bound is derived in the older documents. I see that the error bound of singular value is scaled by

"p(m,n) is a modestly growing function of m and n."

where m and n are matrix dimension. Someone then decide p(m,n) = max(m,n)/ Experimentally I would think the sqrt of that is perhaps a better choice.

[Paul]

I have no actual opininon on this subject.

I did some (i.e., very little) searching around and haven't yet seen a sqrt in setting the default tolerance in other packages. What I have seen is this tolerance, which is not the same as Matlab's documented tolerance

tol = max(size(A)) * norm(A) * eps;

(where norm(A) is expressed as max(svd(A))). source 1 , source 2, source 3 (note that source 3 uses an absolute tolerance as well). I'm assuming that eps in those other packages has the same meaning as in Matlab.

[Bruno Luong]

Why you say it's not the same?

1. norm(A) is max(svd(A))
2. and eps(norm(A)) is norm(A)*eps(1) (or if there is a difference (?), it is negligible)

[Paul]

I was just pointing out that eps(norm(A)) ~= norm(A)*eps

rng(100)
A = 1e10*rand(20);
eps(norm(A))

ans = 1.5259e-05

norm(A)*eps

ans = 2.2118e-05

[Bruno Luong]

Fair enoigh, they differ by factor of 2 at most. As long as it is proportional to norm(A) I'm fine, so both are OK for me.

[Bruno Luong]

https://www.netlib.org/lapack/lug/node97.html eqt UERRBD(i)

https://core.ac.uk/download/pdf/216041984.pdf Lemme 4.2, for definition of modestly growing

[Paul]

Thanks for the links. I checked the LINPACK User's Guide from siam.org and I didn't see the TMW formula in Chapter 11, Section 1 (at least not explicitly).