Cleanup in linalg: subspace angle, distance and Frobenius inner

linalg/frobenius_inner.m

@@ -1,10 +1,21 @@
 function d=frobenius_inner(A,B)
-% FROBENIUS_INNER Short description of frobenius_inner.
-%   FROBENIUS_INNER Long description of frobenius_inner.
+% FROBENIUS_INNER Computes the Frobenius inner product.
+%   D=FROBENIUS_INNER(A,B) computes the Frobenius inner product between
+%   matrices A and B. The Frobenius inner product is given by 
+%     A:B := tr(A'*B) = sum_ij A_ij B_ij
+%
+% Note:
+%   The computation is of course not carried out by the trace formula,
+%   which is very inefficient, but by a more efficient method. 
 %
 % Example (<a href="matlab:run_example frobenius_inner">run</a>)
+%   A = rand(4,5);
+%   B = rand(4,5);
+%   fprintf('A:B=%g\n', frobenius_inner(A, B));
+%   % computing the Frobenius norm
+%   fprintf('||A||_F=%g=%g\n', sqrt(frobenius_inner(A, A)), norm(A,'fro'));
 %
-% See also
+% See also NORM, TRACE
 
 %   Elmar Zander
 %   Copyright 2010, Inst. of Scientific Computing, TU Braunschweig

linalg/subspace_angles.m

@@ -0,0 +1,47 @@
+function [theta,U,V]=subspace_angles(A,B)
+% SUBSPACE_ANGLES Compute the principal angles between subspaces.
+%   THETA=SUBSPACE_ANGLES(A, B) computes the principal angles THETA between
+%   subspaces spanned by the matrices A and B. For a definition see [1]
+%   section 12.4.3.
+%
+%   [THETA,U,V]=SUBSPACE_ANGLES(A, B) also computes the principal vectors
+%   and returns them in U and V.
+%
+% References
+%   [1] G. Golub and C. F. Van Loan, Matrix Computations, 3rd ed. John
+%       Hopkins University Press
+%
+% Example (<a href="matlab:run_example subspace_angles">run</a>)
+%   A = rand(7, 3);
+%   B = rand(7, 5);
+%   [theta, U, V] = subspace_angles(A, B)
+%
+% See also SUBSPACE, SUBSPACE_DISTANCE
+
+%   Elmar Zander
+%   Copyright 2010, Inst. of Scientific Computing, TU Braunschweig
+%
+%   This program 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. 
+%   See the GNU General Public License for more details. You should have
+%   received a copy of the GNU General Public License along with this
+%   program.  If not, see <http://www.gnu.org/licenses/>.
+
+
+[QA,RA]=qr(A,0); %#ok<NASGU>
+[QB,RB]=qr(B,0); %#ok<NASGU>
+C=QA'*QB;
+if nargout>1
+    [Y,S,Z]=svd(C);
+    q=min(size(A,2),size(B,2));
+    U=QA*Y(:,1:q);
+    V=QB*Z(:,1:q);
+    s=diag(S);
+else
+    s=svd(C);
+end
+% clamp s between 0 and 1
+s=max(min(s,1),0);
+theta=acos(s);

linalg/subspace_distance.m

@@ -1,10 +1,35 @@
-function dist=subspace_distance(A,B)
-% SUBSPACE_DISTANCE Short description of subspace_distance.
-%   SUBSPACE_DISTANCE Long description of subspace_distance.
+function dist=subspace_distance(A, B, varargin)
+% SUBSPACE_DISTANCE Compute distance between subspaces.
+%   DIST=SUBSPACE_DISTANCE(A, B) computed the distance between the
+%   subspaces A and B.
 %
+% Options
+%   type: {'mlab'}, 'gv', 'wwf'
+%       Compute the distance compatible to the matlab subspace functions
+%       ('mlab', default), according to [1] section 2.6.3 ('gv'), or
+%       according to 2 ('wwf'). Note, that the last method only works, if
+%       the subspaces have the same size. Note further, that in that case
+%       method 'mlab' and 'gv' coincide.
+%  
+% References
+%   [1] G. Golub and C. F. Van Loan, Matrix Computations, 3rd ed. John
+%       Hopkins University Press, 
+%   [2] Xi-Chen Sun, Xiao Wang, Qiansheng Cheng, Jufu Feng: On Subspace
+%       Distances, LMAM, School of Mathematical Sciences,Peking
+%       www.math.pku.edu.cn:8000/var/preprint/620.pdf‎ 
+% 
 % Example (<a href="matlab:run_example subspace_distance">run</a>)
+%   A=rand(10,3);
+%   B=rand(10,5);
+%   subspace_distance(A, B)
+%   subspace_distance(A, B, 'type', 'gv')
 %
-% See also
+%   B=rand(10,3);
+%   subspace_distance(A, B)
+%   subspace_distance(A, B, 'type', 'gv')
+%   subspace_distance(A, B, 'type', 'wwf')
+%
+% See also SUBSPACE, SUBSPACE_ANGLES
 
 %   Elmar Zander
 %   Copyright 2010, Inst. of Scientific Computing, TU Braunschweig
@@ -18,17 +43,28 @@
 %   received a copy of the GNU General Public License along with this
 %   program.  If not, see <http://www.gnu.org/licenses/>.
 
-if ~iscell(A)
-    A=orth(A);
-    B=orth(B);
-    if size(A,2)<size(B,2)
-        dist=norm(A-B*B'*A);
-    else
-        dist=norm(B-A*A'*B);
-    end
-    % should give the same as the sin(subspace(A1,A2)) where subspace is a
-    % matlab buildin for computing angles between subspaces
-else
-    % A={A1,A2}
-
+options=varargin2options(varargin);
+[type,options]=get_option(options, 'type', 'mlab');
+check_unsupported_options(options, mfilename);
+
+U=orth(A);
+V=orth(B);
+switch type
+    case 'mlab'
+        if size(U,2)<size(V,2)
+            dist=norm(U-V*V'*U, 2);
+        else
+            dist=norm(V-U*U'*V, 2);
+        end
+    case 'gv'
+        dist=norm(V*V'-U*U', 2);
+    case 'wwf'
+        [m, n]=size(U);
+        assert(all(size(V)==[m,n]));
+        dist=sqrt(max(min(m,n)-trace(V'*U*U'*V),0));
+        % Note: probably there is an error in [2] and it should read
+        % min(m,n) instead of max(m,n). Further, we need to clamp the value
+        % to be always positive.
+    otherwise
+        error('sglib:subspace_distance', 'Unknown distance type: %s', type);
 end