Conditional Number growth of Hilbert Matrix: Theoretical vs MatLab. I need to investigate how the condition number of the Hilbert matrix grows with the size N. The Matlab command is: "cond(hilb(N),2)"


*

*Compute the condition number of the Hilbert matrices Hn ∈ R, N×N, for all N = 1, . . . , 50.

*Using Matlab to calculate the log of condition number of Hn versus N. (The blue 'x')

*Compare this with the anticipated theoretical growth (The red line) of $$O\left(\frac{(1+\sqrt{2})^{4N}}{\sqrt{N}}\right) $$


I got a plot like this:

When N = 13, the Condition Number reaches the maximum. The Condition Number does not continue to grow when N>13. Why does the Condition Number stop growing after N=13? 
% generate Hilbert matrices and compute cond number with 2-norm

N=50; % maximum size of a matrix
condofH = []; % conditional number of Hilbert Matrix
N_it= zeros(1,N); 


% compute the cond number of Hn
for n = 1:N
    Hn = hilb(n);
    N_it(n)=n;
    condofH = [condofH cond(Hn,2)];
end

% at this point we have a vector condofH that contains the condition
% number of the Hilber matrices from 1x1 to 50x50.
% plot on the same graph the theoretical growth line. 


% Theoretical growth of condofH
x = 1:50;
y = (1+sqrt(2)).^(4*x)./(sqrt(x));

% plot
plot(N_it, log(y));
plot(N_it, log(condofH),'x', N_it,log(y));


% plot labels
plot(N_it, log(condofH),'x', N_it,log(y))
title('Conditional Number growth of Hilbert Matrix: Theoretical vs Matlab')
xlabel('N', 'fontsize', 16)
ylabel('log(cond(Hn))','fontsize', 16)
lgd = legend ('Location', 'northwest')
legend('MatLab', 'Theoretical')
legend('show')

 A: In order to understand what exactly happens here, one has to take a close look at how the command cond is implemented in details, which I have not done (see the update at the end). However, one possible explanation for the phenomenon is the following: suppose you have to calculate the function $g(x)=x$, and the calculation is implemented as $g(x)=1/f(x)$ where $f(x)=1/x$. The function $g(x)$ has the linear growth (theoretically), but with floating point calculations for large $x$ (when $1/x$ becomes smaller than the MATLAB machine epsilon eps=$2^{-52}$) you will get round-off errors dominating in the denominator
$$
g(x)\approx\frac{1}{\frac{1}{x}+{\rm eps}}\approx \frac{1}{{\rm eps}}\approx 2^{52}.
$$
I guess that something similar happens in the command cond.
Update: After having a look at cond.m I confirm the phenomenon. MATLAB calculates the condition number for $p=2$ via SVD, that is, for Hilbert matrices it becomes
$$
k(H_n)=\frac{\lambda_\max(H_n)}{\lambda_\min(H_n)}.
$$
When $n\to\infty$ we have $\lambda_\max\to\pi$ (e.g. see this question) and $\lambda_\min\to 0$. Calculation with floating points for large $n$ looks as
$$
k(H_n)\approx\frac{\pi}{\lambda_\min+{\rm eps}}\approx\frac{\pi}{\rm eps}\approx 2^{52}.
$$
A: MatLab calculates Condition Number in this way:
$$Cond(H_n) = ||H_n^{−1}||⋅||H_n|| $$
One explanation for this is that MatLab uses: invhilb(n) to generate Hn-1. According to the MatLab documentation, invhilb(n) gives exact inverse of the Hilbert matrix without roundoff error when N<15. For N>=15, invhilb(n) involves round-off errors and Cond(Hn) Algorithm starts to loose robustness. This explains why Log(cond(Hn)) does not continue to grow after N=14.
For N>=15, I predict MatLab generates Hn-1 using inv(Hn) The documentation indicates that inv() performs an LU decomposition of Hn. Therefore this algorithm does round-off in $O(2/3N^3)$ times. It suffers from accumulation of roundoff errors.
