Mercurial > hg > octave-lyh

comparison scripts/special-matrix/invhilb.m @ 3889:ac24529a78a0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2002-04-04 23:03:15 by jwe]
author jwe
date Thu, 04 Apr 2002 23:03:15 +0000
parents ae7adbb591e8
children 2168f4a0e88d
comparison
equal deleted inserted replaced
3888:70ebd3d672a1 3889:ac24529a78a0
1 ## Copyright (C) 1996, 1997 John W. Eaton 1 ## Copyright (C) 2002 Dirk Laurie
2 ## 2 ##
3 ## This file is part of Octave. 3 ## This file is part of Octave.
4 ## 4 ##
5 ## Octave is free software; you can redistribute it and/or modify it 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 6 ## under the terms of the GNU General Public License as published by
17 ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA 17 ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
18 ## 02111-1307, USA. 18 ## 02111-1307, USA.
19 19
20 ## -*- texinfo -*- 20 ## -*- texinfo -*-
21 ## @deftypefn {Function File} {} invhilb (@var{n}) 21 ## @deftypefn {Function File} {} invhilb (@var{n})
22 ## Return the inverse of a Hilbert matrix of order @var{n}. This is exact. 22 ## Return the inverse of a Hilbert matrix of order @var{n}. This can be
23 ## Compare with the numerical calculation of @code{inverse (hilb (n))}, 23 ## computed computed exactly using
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))},
24 ## which suffers from the ill-conditioning of the Hilbert matrix, and the 62 ## which suffers from the ill-conditioning of the Hilbert matrix, and the
25 ## finite precision of your computer's floating point arithmetic. 63 ## finite precision of your computer's floating point arithmetic.
64 ##
26 ## @end deftypefn 65 ## @end deftypefn
27 ## @seealso{hankel, vander, sylvester_matrix, hilb, and toeplitz} 66 ## @seealso{hankel, vander, sylvester_matrix, hilb, and toeplitz}
28 67
29 ## Author: jwe 68 ## Author: Dirk Laurie <dirk@siegfried.wisk.sun.ac.za>
30 69
31 function retval = invhilb (n) 70 function retval = invhilb (n)
32 71
33 if (nargin != 1) 72 if (nargin != 1)
34 usage ("invhilb (n)"); 73 usage ("invhilb (n)");
35 endif 74 endif
36 75
37 nmax = length (n); 76 nmax = length (n);
38 if (nmax == 1) 77 if (nmax == 1)
39 retval = zeros (n); 78
40 for l = 1:n 79 ## The point about the second formula above is that when vectorized,
41 for k = l:n 80 ## p(k) is evaluated for k=1:n which involves O(n) calls to bincoeff
42 tmp = 1; 81 ## instead of O(n^2).
43 for i = 1:n 82 ##
44 tmp = tmp * (i + k - 1); 83 ## We evaluate the expression as (-1)^(i+j)*(p(i)*p(j))/(i+j-1) except
45 endfor 84 ## when p(i)*p(j) would overflow. In cases where p(i)*p(j) is an exact
46 for i = 1:n 85 ## machine number, the result is also exact. Otherwise we calculate
47 if (i != k) 86 ## (-1)^(i+j)*p(i)*(p(j)/(i+j-1)).
48 tmp = tmp * (l + i - 1); 87 ##
49 endif 88 ## The Octave bincoeff routine uses transcendental functions (lgamma
50 endfor 89 ## and exp) rather than multiplications, for the sake of speed.
51 for i = 1:n 90 ## However, it rounds the answer to the nearest integer, which
52 if (i != l) 91 ## justifies the claim about exactness made above.
53 tmp = tmp / (i - l); 92
54 endif 93 retval = zeros (n);
55 endfor 94 k = [1:n];
56 for i = 1:n 95 p = k .* bincoeff (k+n-1, k-1) .* bincoeff (n, k);
57 if (i != k) 96 p(2:2:n) = -p(2:2:n);
58 tmp = tmp / (i - k); 97 if (n < 203)
59 endif 98 for l = 1:n
60 endfor 99 retval(l,:) = (p(l) * p) ./ [l:l+n-1];
61 retval (k, l) = tmp;
62 retval (l, k) = tmp;
63 endfor 100 endfor
64 endfor 101 else
102 for l = 1:n
103 retval(l,:) = p(l) * (p ./ [l:l+n-1]);
104 endfor
105 endif
65 else 106 else
66 error ("hilb: expecting scalar argument, found something else"); 107 error ("invhilb: expecting scalar argument, found something else");
67 endif 108 endif
68 109
69 endfunction 110 endfunction