octave_packages/linear-algebra-2.2.0/smwsolve.m

   1 ## Copyright (C) 2009 VZLU Prague, a.s., Czech Republic
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn{Function File} {@var{x} =} smwsolve (@var{a}, @var{u}, @var{v}, @var{b})
  18 ## @deftypefnx{Function File} {} smwsolve (@var{solver}, @var{u}, @var{v}, @var{b})
  19 ## Solves the square system @code{(A + U*V')*X == B}, where @var{u} and @var{v} are
  20 ## matrices with several columns, using the Sherman-Morrison-Woodbury formula,
  21 ## so that a system with @var{a} as left-hand side is actually solved. This is
  22 ## especially advantageous if @var{a} is diagonal, sparse, triangular or
  23 ## positive definite.
  24 ## @var{a} can be sparse or full, the other matrices are expected to be full.
  25 ## Instead of a matrix @var{a}, a user may alternatively provide a function
  26 ## @var{solver} that performs the left division operation.
  27 ## @end deftypefn
  28
  29 ## Author: Jaroslav Hajek <highegg@gmail.com>
  30
  31 function x = smwsolve (a, u, v, b)
  32
  33   if (nargin != 4)
  34     print_usage ();
  35   endif
  36
  37   n = columns (u);
  38
  39   if (n != columns (v) || rows (a) != rows (u) || columns (a) != rows (v))
  40     error ("smwsolve: dimension mismatch");
  41   elseif (! issquare (a))
  42     error ("smwsolve: need a square matrix");
  43   endif
  44
  45
  46   nc = columns (b);
  47   n = columns (u);
  48
  49   if (ismatrix (a))
  50     xx = a \ [b, u];
  51   elseif (isa (a, "function_handle"))
  52     xx = a ([b, u]);
  53     if (rows (xx) != rows (a) || columns (xx) != (nc + n))
  54       error ("smwsolve: invalid result from a solver function");
  55     endif
  56   else
  57     error ("smwsolve: a must be a matrix or function handle");
  58   endif
  59
  60   x = xx(:,1:nc);
  61   y = xx(:,nc+1:nc+n);
  62
  63   vxx = v' * xx;
  64   vx = vxx(:,1:nc);
  65   vy = vxx(:,nc+1:nc+n);
  66
  67   x = x - y * ((eye (n) + vy) \ vx);
  68
  69 endfunction
  70
  71 %!test
  72 %! A = 2.1*eye (10);
  73 %! u = rand (10, 2); u /= diag (norm (u, "cols"));
  74 %! v = rand (10, 2); v /= diag (norm (v, "cols"));
  75 %! b = rand (10, 2);
  76 %! x1 = (A + u*v') \ b;
  77 %! x2 = smwsolve (A, u, v, b);
  78 %! assert (x1, x2, 1e-13);