|
3889
|
1 ## Copyright (C) 2002 Dirk Laurie |
|
2313
|
2 ## |
|
|
3 ## This file is part of Octave. |
|
|
4 ## |
|
|
5 ## Octave is free software; you can redistribute it and/or modify it |
|
|
6 ## under the terms of the GNU General Public License as published by |
|
|
7 ## the Free Software Foundation; either version 2, or (at your option) |
|
|
8 ## any later version. |
|
|
9 ## |
|
|
10 ## Octave is distributed in the hope that it will be useful, but |
|
|
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
13 ## General Public License for more details. |
|
|
14 ## |
|
|
15 ## You should have received a copy of the GNU General Public License |
|
|
16 ## along with Octave; see the file COPYING. If not, write to the Free |
|
5307
|
17 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
|
|
18 ## 02110-1301, USA. |
|
245
|
19 |
|
3369
|
20 ## -*- texinfo -*- |
|
|
21 ## @deftypefn {Function File} {} invhilb (@var{n}) |
|
3889
|
22 ## Return the inverse of a Hilbert matrix of order @var{n}. This can be |
|
5435
|
23 ## computed exactly using |
|
3889
|
24 ## @tex |
|
|
25 ## $$\eqalign{ |
|
|
26 ## A_{ij} &= -1^{i+j} (i+j-1) |
|
|
27 ## \left( \matrix{n+i-1 \cr n-j } \right) |
|
|
28 ## \left( \matrix{n+j-1 \cr n-i } \right) |
|
|
29 ## \left( \matrix{i+j-2 \cr i-2 } \right)^2 \cr |
|
|
30 ## &= { p(i)p(j) \over (i+j-1) } |
|
|
31 ## }$$ |
|
|
32 ## where |
|
|
33 ## $$ |
|
|
34 ## p(k) = -1^k \left( \matrix{ k+n-1 \cr k-1 } \right) |
|
|
35 ## \left( \matrix{ n \cr k } \right) |
|
|
36 ##$$ |
|
|
37 ## @end tex |
|
|
38 ## @ifinfo |
|
|
39 ## @example |
|
|
40 ## |
|
|
41 ## (i+j) /n+i-1\ /n+j-1\ /i+j-2\ 2 |
|
|
42 ## A(i,j) = -1 (i+j-1)( )( ) ( ) |
|
|
43 ## \ n-j / \ n-i / \ i-2 / |
|
|
44 ## |
|
|
45 ## = p(i) p(j) / (i+j-1) |
|
|
46 ## |
|
|
47 ## @end example |
|
|
48 ## where |
|
|
49 ## @example |
|
|
50 ## k /k+n-1\ /n\ |
|
|
51 ## p(k) = -1 ( ) ( ) |
|
|
52 ## \ k-1 / \k/ |
|
|
53 ## @end example |
|
|
54 ## @end ifinfo |
|
|
55 ## |
|
|
56 ## The validity of this formula can easily be checked by expanding |
|
|
57 ## the binomial coefficients in both formulas as factorials. It can |
|
|
58 ## be derived more directly via the theory of Cauchy matrices: |
|
|
59 ## see J. W. Demmel, Applied Numerical Linear Algebra, page 92. |
|
|
60 ## |
|
|
61 ## Compare this with the numerical calculation of @code{inverse (hilb (n))}, |
|
3369
|
62 ## which suffers from the ill-conditioning of the Hilbert matrix, and the |
|
|
63 ## finite precision of your computer's floating point arithmetic. |
|
5642
|
64 ## @seealso{hankel, vander, sylvester_matrix, hilb, toeplitz} |
|
3369
|
65 ## @end deftypefn |
|
4
|
66 |
|
5132
|
67 ## Author: Dirk Laurie <dlaurie@na-net.ornl.gov> |
|
2314
|
68 |
|
2311
|
69 function retval = invhilb (n) |
|
4
|
70 |
|
|
71 if (nargin != 1) |
|
904
|
72 usage ("invhilb (n)"); |
|
4
|
73 endif |
|
|
74 |
|
|
75 nmax = length (n); |
|
|
76 if (nmax == 1) |
|
3889
|
77 |
|
|
78 ## The point about the second formula above is that when vectorized, |
|
|
79 ## p(k) is evaluated for k=1:n which involves O(n) calls to bincoeff |
|
|
80 ## instead of O(n^2). |
|
|
81 ## |
|
|
82 ## We evaluate the expression as (-1)^(i+j)*(p(i)*p(j))/(i+j-1) except |
|
|
83 ## when p(i)*p(j) would overflow. In cases where p(i)*p(j) is an exact |
|
|
84 ## machine number, the result is also exact. Otherwise we calculate |
|
|
85 ## (-1)^(i+j)*p(i)*(p(j)/(i+j-1)). |
|
|
86 ## |
|
4031
|
87 ## The Octave bincoeff routine uses transcendental functions (gammaln |
|
3889
|
88 ## and exp) rather than multiplications, for the sake of speed. |
|
|
89 ## However, it rounds the answer to the nearest integer, which |
|
|
90 ## justifies the claim about exactness made above. |
|
|
91 |
|
|
92 retval = zeros (n); |
|
|
93 k = [1:n]; |
|
|
94 p = k .* bincoeff (k+n-1, k-1) .* bincoeff (n, k); |
|
|
95 p(2:2:n) = -p(2:2:n); |
|
|
96 if (n < 203) |
|
|
97 for l = 1:n |
|
|
98 retval(l,:) = (p(l) * p) ./ [l:l+n-1]; |
|
4
|
99 endfor |
|
3889
|
100 else |
|
|
101 for l = 1:n |
|
|
102 retval(l,:) = p(l) * (p ./ [l:l+n-1]); |
|
|
103 endfor |
|
|
104 endif |
|
4
|
105 else |
|
3889
|
106 error ("invhilb: expecting scalar argument, found something else"); |
|
4
|
107 endif |
|
|
108 |
|
|
109 endfunction |
