Wed May 20 19:35:57 CDT 2009
2009/5/20 Stéfan van der Walt <firstname.lastname@example.org>:
> 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:
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.
"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
More information about the Scipy-dev