octave_packages/m/special-matrix/invhilb.m

   1 ## Copyright (C) 1993-2012 Dirk Laurie
   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 3 of the License, or (at
   8 ## your option) 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, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {} invhilb (@var{n})
  21 ## Return the inverse of the Hilbert matrix of order @var{n}.  This can be
  22 ## computed exactly using
  23 ## @tex
  24 ## $$\eqalign{
  25 ##   A_{ij} &= -1^{i+j} (i+j-1)
  26 ##              \left( \matrix{n+i-1 \cr n-j } \right)
  27 ##              \left( \matrix{n+j-1 \cr n-i } \right)
  28 ##              \left( \matrix{i+j-2 \cr i-2 } \right)^2 \cr
  29 ##          &= { p(i)p(j) \over (i+j-1) }
  30 ## }$$
  31 ## where
  32 ## $$
  33 ##   p(k) = -1^k \left( \matrix{ k+n-1 \cr k-1 } \right)
  34 ##               \left( \matrix{ n \cr k } \right)
  35 ## $$
  36 ## @end tex
  37 ## @ifnottex
  38 ##
  39 ## @example
  40 ## @group
  41 ##
  42 ##             (i+j)         /n+i-1\  /n+j-1\   /i+j-2\ 2
  43 ##  A(i,j) = -1      (i+j-1)(       )(       ) (       )
  44 ##                           \ n-j /  \ n-i /   \ i-2 /
  45 ##
  46 ##         = p(i) p(j) / (i+j-1)
  47 ##
  48 ## @end group
  49 ## @end example
  50 ##
  51 ## @noindent
  52 ## where
  53 ##
  54 ## @example
  55 ## @group
  56 ##              k  /k+n-1\   /n\
  57 ##     p(k) = -1  (       ) (   )
  58 ##                 \ k-1 /   \k/
  59 ## @end group
  60 ## @end example
  61 ##
  62 ## @end ifnottex
  63 ## The validity of this formula can easily be checked by expanding
  64 ## the binomial coefficients in both formulas as factorials.  It can
  65 ## be derived more directly via the theory of Cauchy matrices.
  66 ## See J. W. Demmel, @cite{Applied Numerical Linear Algebra}, p. 92.
  67 ##
  68 ## Compare this with the numerical calculation of @code{inverse (hilb (n))},
  69 ## which suffers from the ill-conditioning of the Hilbert matrix, and the
  70 ## finite precision of your computer's floating point arithmetic.
  71 ## @seealso{hilb}
  72 ## @end deftypefn
  73
  74 ## Author: Dirk Laurie <dlaurie@na-net.ornl.gov>
  75
  76 function retval = invhilb (n)
  77
  78   if (nargin != 1)
  79     print_usage ();
  80   elseif (! isscalar (n))
  81     error ("invhilb: N must be a scalar integer");
  82   endif
  83
  84   ## The point about the second formula above is that when vectorized,
  85   ## p(k) is evaluated for k=1:n which involves O(n) calls to bincoeff
  86   ## instead of O(n^2).
  87   ##
  88   ## We evaluate the expression as (-1)^(i+j)*(p(i)*p(j))/(i+j-1) except
  89   ## when p(i)*p(j) would overflow.  In cases where p(i)*p(j) is an exact
  90   ## machine number, the result is also exact.  Otherwise we calculate
  91   ## (-1)^(i+j)*p(i)*(p(j)/(i+j-1)).
  92   ##
  93   ## The Octave bincoeff routine uses transcendental functions (gammaln
  94   ## and exp) rather than multiplications, for the sake of speed.
  95   ## However, it rounds the answer to the nearest integer, which
  96   ## justifies the claim about exactness made above.
  97
  98   retval = zeros (n);
  99   k = [1:n];
 100   p = k .* bincoeff (k+n-1, k-1) .* bincoeff (n, k);
 101   p(2:2:n) = -p(2:2:n);
 102   if (n < 203)
 103     for l = 1:n
 104       retval(l,:) = (p(l) * p) ./ [l:l+n-1];
 105     endfor
 106   else
 107     for l = 1:n
 108       retval(l,:) = p(l) * (p ./ [l:l+n-1]);
 109     endfor
 110   endif
 111
 112 endfunction
 113
 114
 115 %!assert (invhilb (1), 1)
 116 %!assert (invhilb (2), [4, -6; -6, 12])
 117 %!test
 118 %! result4 = [16  , -120 , 240  , -140;
 119 %!            -120, 1200 , -2700, 1680;
 120 %!            240 , -2700, 6480 , -4200;
 121 %!            -140, 1680 , -4200, 2800];
 122 %! assert (invhilb (4), result4);
 123 %!assert (abs (invhilb (7) * hilb (7) - eye (7)) < sqrt (eps))
 124
 125 %!error invhilb ()
 126 %!error invhilb (1, 2)
 127 %!error <N must be a scalar integer> invhilb ([1, 2])
 128