Return the inverse of the Hilbert matrix of order n.
This can be computed exactly using
(i+j) /n+i-1\ /n+j-1\ /i+j-2\ 2 A(i,j) = -1 (i+j-1)( )( ) ( ) \ n-j / \ n-i / \ i-2 / = p(i) p(j) / (i+j-1)
where
k /k+n-1\ /n\ p(k) = -1 ( ) ( ) \ k-1 / \k/
The validity of this formula can easily be checked by expanding the binomial coefficients in both formulas as factorials. It can be derived more directly via the theory of Cauchy matrices. See J. W. Demmel, Applied Numerical Linear Algebra, p. 92.
Compare this with the numerical calculation of inverse (hilb (n))
,
which suffers from the ill-conditioning of the Hilbert matrix, and the
finite precision of your computer’s floating point arithmetic.
See also: hilb.
Package: octave