[Math] Conditional Number growth of Hilbert Matrix: Theoretical vs MatLab.

hilbert-matricesMATLABmatrices

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)"

  1. Compute the condition number of the Hilbert matrices Hn ∈ R, N×N, for all N = 1, . . . , 50.
  2. Using Matlab to calculate the log of condition number of Hn versus N. (The blue 'x')
  3. 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:

enter image description here

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')

Best Answer

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}. $$

Related Question