Mercurial > hg > octave-nkf
annotate scripts/linear-algebra/subspace.m @ 9565:fe57b638e48c
adapt octave-bug.cc.in and mkoctfile.cc.in to recent configure changes
| author | Benjamin Lindner <lindnerb@users.sourceforge.net> |
|---|---|
| date | Wed, 26 Aug 2009 11:31:26 -0400 |
| parents | eb63fbe60fab |
| children | c776f063fefe |
| rev | line source |
|---|---|
| 8920 | 1 ## Copyright (C) 2008, 2009 VZLU Prague, a.s., Czech Republic |
| 7611 | 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} {@var{angle} =} subspace (@var{a}, @var{B}) | |
| 21 ## Determine the largest principal angle between two subspaces | |
| 22 ## spanned by columns of matrices @var{a} and @var{b}. | |
| 23 ## @end deftypefn | |
| 24 | |
| 25 ## Author: Jaroslav Hajek <highegg@gmail.com> | |
| 26 | |
| 27 ## reference: | |
| 28 ## [1] Andrew V. Knyazev, Merico E. Argentati: | |
| 29 ## Principal Angles between Subspaces in an A-Based Scalar Product: | |
| 30 ## Algorithms and Perturbation Estimates. | |
| 31 ## SIAM Journal on Scientific Computing, Vol. 23 no. 6, pp. 2008-2040 | |
| 32 ## | |
| 33 ## other texts are also around... | |
| 34 | |
| 35 function ang = subspace (a, b) | |
| 36 | |
|
7612
c1702f963a5e
error check for subspace.m
Jaroslav Hajek <highegg@gmail.com>
parents:
7611
diff
changeset
|
37 if (nargin != 2) |
|
c1702f963a5e
error check for subspace.m
Jaroslav Hajek <highegg@gmail.com>
parents:
7611
diff
changeset
|
38 print_usage (); |
|
c1702f963a5e
error check for subspace.m
Jaroslav Hajek <highegg@gmail.com>
parents:
7611
diff
changeset
|
39 elseif (ndims (a) != 2 || ndims (b) != 2) |
|
c1702f963a5e
error check for subspace.m
Jaroslav Hajek <highegg@gmail.com>
parents:
7611
diff
changeset
|
40 error ("subspace: expecting A and B to be 2-dimensional arrays"); |
|
7616
fb4fa07bc364
more checks for subspace
Jaroslav Hajek <highegg@gmail.com>
parents:
7612
diff
changeset
|
41 elseif (rows (a) != rows (b)) |
| 8664 | 42 error ("subspace: column dimensions of a and b must match"); |
|
7612
c1702f963a5e
error check for subspace.m
Jaroslav Hajek <highegg@gmail.com>
parents:
7611
diff
changeset
|
43 endif |
|
c1702f963a5e
error check for subspace.m
Jaroslav Hajek <highegg@gmail.com>
parents:
7611
diff
changeset
|
44 |
| 7611 | 45 a = orth (a); |
| 46 b = orth (b); | |
| 47 c = a'*b; | |
| 48 scos = min (svd (c)); | |
| 49 if (scos^2 > 1/2) | |
| 50 if (columns (a) >= columns (b)) | |
| 51 c = b - a*c; | |
| 52 else | |
| 53 c = a - b*c'; | |
| 54 endif | |
| 55 ssin = max (svd (c)); | |
| 56 ang = asin (min (ssin, 1)); | |
| 57 else | |
| 58 ang = acos (scos); | |
| 59 endif | |
| 60 | |
| 61 endfunction |