view scripts/linear-algebra/subspace.m @ 27978:a4268efb7334

maint: Eliminate single "Author:" lines from code base. * color-picker.cc, color-picker.h, octave-qscintilla.cc, octave-qscintilla.h, octave-txt-lexer.cc, octave-txt-lexer.h, balance.cc, kron.cc, lookup.cc, ls-hdf5.cc, ls-mat5.cc, ls-oct-text.cc, sylvester.cc, tsearch.cc, jit-ir.cc, jit-ir.h, jit-typeinfo.cc, jit-typeinfo.h, jit-util.cc, jit-util.h, pt-jit.cc, pt-jit.h, Array-jit.cc, bsxfun.h, url-transfer.cc, url-transfer.h, acosd.m, acotd.m, acscd.m, asecd.m, asind.m, atand.m, cosd.m, cotd.m, cscd.m, secd.m, sind.m, tand.m, cart2pol.m, cart2sph.m, curl.m, del2.m, divergence.m, flip.m, pol2cart.m, randi.m, sortrows.m, sph2cart.m, griddata3.m, griddatan.m, rectint.m, voronoi.m, dialog.m, guidata.m, guihandles.m, __file_filter__.m, __fltk_file_filter__.m, __get_funcname__.m, __is_function__.m, __uigetdir_fltk__.m, __uigetfile_fltk__.m, __uiobject_split_args__.m, __uiputfile_fltk__.m, uibuttongroup.m, uicontextmenu.m, uicontrol.m, uigetdir.m, uigetfile.m, uimenu.m, uipanel.m, uipushtool.m, uiputfile.m, uiresume.m, uitoggletool.m, uitoolbar.m, uiwait.m, doc.m, autumn.m, bone.m, cmpermute.m, cmunique.m, cool.m, copper.m, cubehelix.m, flag.m, frame2im.m, hot.m, hsv.m, hsv2rgb.m, im2frame.m, imfinfo.m, imformats.m, iscolormap.m, jet.m, pink.m, prism.m, __imfinfo__.m, rainbow.m, rgb2hsv.m, spinmap.m, spring.m, viridis.m, white.m, winter.m, importdata.m, strmatch.m, bandwidth.m, isbanded.m, isdiag.m, istril.m, istriu.m, krylov.m, linsolve.m, logm.m, lscov.m, subspace.m, bunzip2.m, compare_versions.m, edit.m, fullfile.m, getfield.m, gunzip.m, inputParser.m, license.m, nargchk.m, narginchk.m, nthargout.m, python.m, setfield.m, tar.m, unpack.m, ver.m, zip.m, odeplot.m, fminbnd.m, fminunc.m, fsolve.m, fzero.m, glpk.m, humps.m, pathdef.m, savepath.m, diffuse.m, lighting.m, material.m, shading.m, specular.m, contour.m, ellipsoid.m, isocaps.m, isocolors.m, isonormals.m, isosurface.m, light.m, pcolor.m, __interp_cube__.m, __marching_cube__.m, __patch__.m, __stem__.m, reducepatch.m, reducevolume.m, ribbon.m, shrinkfaces.m, slice.m, stem.m, surf.m, surfl.m, tetramesh.m, waterfall.m, __actual_axis_position__.m, __next_line_color__.m, allchild.m, figure.m, findall.m, findobj.m, gcf.m, hggroup.m, isprop.m, ndgrid.m, __add_default_menu__.m, saveas.m, __splinefit__.m, prefdir.m, preferences.m, profexplore.m, profexport.m, profile.m, profshow.m, setdiff.m, freqz.m, freqz_plot.m, sinc.m, unwrap.m, pcr.m, qmr.m, sprandn.m, betaincinv.m, cosint.m, ellipke.m, factor.m, gammaincinv.m, legendre.m, nchoosek.m, primes.m, invhilb.m, toeplitz.m, erase.m, mat2str.m, strtrim.m, validatestring.m, fail.m, compare_plot_demos.m, dump_demos.m, html_compare_plot_demos.m, addtodate.m, datenum.m, jit.tst: Eliminate single "Author:" lines from code base.
author Rik <rik@octave.org>
date Tue, 21 Jan 2020 12:53:13 -0800
parents bd51beb6205e
children b09432b20a84
line wrap: on
line source

########################################################################
##
## Copyright (C) 2008-2020 The Octave Project Developers
##
## See the file COPYRIGHT.md in the top-level directory of this
## distribution or <https://octave.org/copyright/>.
##
## This file is part of Octave.
##
## Octave is free software: you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <https://www.gnu.org/licenses/>.
##
########################################################################

## -*- texinfo -*-
## @deftypefn {} {@var{angle} =} subspace (@var{A}, @var{B})
## Determine the largest principal angle between two subspaces
## spanned by the columns of matrices @var{A} and @var{B}.
## @end deftypefn

## reference:
## [1]  Andrew V. Knyazev, Merico E. Argentati:
## Principal Angles between Subspaces in an A-Based Scalar Product:
## Algorithms and Perturbation Estimates.
## SIAM Journal on Scientific Computing, Vol. 23 no. 6, pp. 2008-2040
##
## other texts are also around...

function ang = subspace (A, B)

  if (nargin != 2)
    print_usage ();
  elseif (ndims (A) != 2 || ndims (B) != 2)
    error ("subspace: A and B must be 2-dimensional arrays");
  elseif (rows (A) != rows (B))
    error ("subspace: column dimensions of A and B must match");
  endif

  A = orth (A);
  B = orth (B);
  c = A'*B;
  scos = min (svd (c));
  if (scos^2 > 1/2)
    if (columns (A) >= columns (B))
      c = B - A*c;
    else
      c = A - B*c';
    endif
    ssin = max (svd (c));
    ang = asin (min (ssin, 1));
  else
    ang = acos (scos);
  endif

endfunction


%!assert (subspace (1, 1), 0)
%!assert (subspace ([1, 0]', [1, 1; 0, 1]'), 0, 3*eps)
%!assert (subspace ([1, 0, 1]', [1, 1, 0; 1, -1, 0]'), pi/4, 3*eps)
%!assert (subspace ([1 5 0 0; -3 2 0 0]', [0 0 4 2; 0 0 4 3]'), pi/2)
%!assert (subspace ([1 1 1 1; 1 2 3 4]', [1 -1 -1 1]'), pi/2)

%!test
%! ## For small angle between subspaces
%! theta = pi/200;
%! Ry = [cos(theta), 0, sin(theta);0, 1, 0;-sin(theta), 0, cos(theta)];
%! a = Ry*[3*e, 0, 0]';
%! b = [1, 1, 0; 1, -1, 0]';
%! assert (theta, subspace (a, b), eps);