octave_packages/linear-algebra-2.2.0/cod.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{q}, @var{r}, @var{z}] =} cod (@var{a})
  18 ## @deftypefnx{Function File} {[@var{q}, @var{r}, @var{z}, @var{p}] =} cod (@var{a})
  19 ## @deftypefnx{Function File} {[@dots{}] =} cod (@var{a}, '0')
  20 ## Computes the complete orthogonal decomposition (COD) of the matrix @var{a}:
  21 ## @example
  22 ##   @var{a} = @var{q}*@var{r}*@var{z}'
  23 ## @end example
  24 ## Let @var{a} be an M-by-N matrix, and let @code{K = min(M, N)}.
  25 ## Then @var{q} is M-by-M orthogonal, @var{z} is N-by-N orthogonal,
  26 ## and @var{r} is M-by-N such that @code{@var{r}(:,1:K)} is upper
  27 ## trapezoidal and @code{@var{r}(:,K+1:N)} is zero.
  28 ## The additional @var{p} output argument specifies that pivoting should be used in
  29 ## the first step (QR decomposition). In this case,
  30 ## @example
  31 ##   @var{a}*@var{p} = @var{q}*@var{r}*@var{z}'
  32 ## @end example
  33 ## If a second argument of '0' is given, an economy-sized factorization is returned
  34 ## so that @var{r} is K-by-K.
  35 ##
  36 ## @emph{NOTE}: This is currently implemented by double QR factorization plus some
  37 ## tricky manipulations, and is not as efficient as using xRZTZF from LAPACK.
  38 ## @seealso{qr}
  39 ## @end deftypefn
  40
  41 ## Author: Jaroslav Hajek <highegg@gmail.com>
  42
  43 function [q, r, z, p] = cod (a, varargin)
  44
  45   if (nargin < 1 || nargin > 2 || nargout > 4 || ! ismatrix (a))
  46     print_usage ();
  47   endif
  48
  49   [m, n] = size (a);
  50   k = min (m, n);
  51   economy = nargin == 2;
  52   pivoted = nargout == 4;
  53
  54   ## Compute the initial QR decomposition
  55   if (pivoted)
  56     [q, r, p] = qr (a, varargin{:});
  57   else
  58     [q, r] = qr (a, varargin{:});
  59   endif
  60
  61   if (m >= n)
  62     ## In this case, Z is identity, and we're finished.
  63     z = eye (n, class (a));
  64   else
  65     ## Permutation inverts row order.
  66     pr = eye (m) (m:-1:1, :);
  67     ## Permutation inverts first m columns order.
  68     pc = eye (n) ([m:-1:1, m+1:n], :);
  69     ## Make n-by-m matrix, invert first m columns
  70     r = (pr * r * pc')';
  71     ## QR factorize again.
  72     [z, r] = qr (r, varargin{:});
  73     ## Recover final R and Z
  74     if (economy)
  75       r = pr * r' * pr';
  76       z = pc * z * pr';
  77     else
  78       r = pr * r' * pc';
  79       z = pc * z * pc';
  80     endif
  81   endif
  82
  83 endfunction
  84
  85 %!test
  86 %! a = rand (5, 10);
  87 %! [q, r, z] = cod (a);
  88 %! assert (norm (q*r*z' - a) / norm (a) < 1e-10);
  89 %!test
  90 %! a = rand (5, 10) + i * rand (5, 10);
  91 %! [q, r, z] = cod (a);
  92 %! assert (norm (q*r*z' - a) / norm (a) < 1e-10);
  93 %!test
  94 %! a = rand (5, 10);
  95 %! [q, r, z, p] = cod (a);
  96 %! assert (norm (q*r*z' - a*p) / norm (a) < 1e-10);