[SciPy-dev] Hilbert-matrix

[Robert Kern]

2009/5/20 Stéfan van der Walt <stefan at sun.ac.za>:
> Hi all,
>
> The Ill-conditioned Hilbert-matrix is often studied by students in
> numerical linear algebra, and as such I thought it might deserve a
> place in scipy.linalg.
>
> def hilbert(N):
>    x, y = np.ogrid[1.:N+1, 1.:N+1]
>    return 1 / (x + y - 1)
>
> The trickier issue, however, would be calculating the inverse -- does
> anybody have suggestions on how to implement it accurately?

There is an analytic formula involving combinations:

  http://mathworld.wolfram.com/HilbertMatrix.html

You could probably implement that straightforwardly in
double-precision with gammaln followed by exponentiation. That would
probably be about as good as any double-precision inversion routine
should be expected to get, best-case, since for largish matrices
double-precision won't even accurately represent the correct answer. A
fully-correct bigint expansion is neat, but not particularly necessary
to demonstrate instability.

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
  -- Umberto Eco