子空間夾角計(jì)算-MATLAB源代碼修正
2017-01-16 by:CAE仿真在線 來源:互聯(lián)網(wǎng)
function theta = subspace(A,B)
%SUBSPACE Angle between subspaces.
% SUBSPACE(A,B) finds the angle
between two subspaces specified by the
% columns of A and B.
%
% If the angle is small, the
two spaces are nearly linearly dependent.
% In a physical experiment
described by some observations A, and a second
% realization of the experiment
described by B, SUBSPACE(A,B) gives a
% measure of the amount of new
information afforded by the second
% experiment not associated
with statistical errors of fluctuations.
%
% Class support for inputs A,
B:
%
float: double, single
% The algorithm used here
ensures that small angles are computed
% accurately, and it allows
subspaces of different dimensions following
% the definition in [2]. The
first issue is crucial. The second issue is
% not so important; but since
the definition from [2] coinsides with the
% standard definition when the
dimensions are equal, there should be no
% confusion - and subspaces
with different dimensions may arise in
% problems where the dimension
is computed as the numerical rank of some
% inaccurate matrix.
% References:
% [1] A. Bjorck
& G. Golub, Numerical methods for computing
%
angles between linear subspaces, Math. Comp. 27 (1973),
%
pp. 579-594.
% [2] P.-A. Wedin, On angles
between subspaces of a finite
%
dimensional inner product space, in B. Kagstrom &
A. Ruhe (Eds.),
%
Matrix Pencils, Lecture Notes in Mathematics 973, Springer,
1983,
%
pp. 263-285.
% Thanks to Per Christian
Hansen
% Copyright 1984-2007 The
MathWorks, Inc.
% $Revision: 5.10.4.3
$ $Date: 2007/09/18 02:15:38 $
% Compute orthonormal bases, using SVD in "orth" to avoid
problems
% when A and/or B is nearly rank deficient.
A = orth(A);
B = orth(B);
%Check rank and swap
if size(A,2) < size(B,2)
tmp = A; A = B; B = tmp;
end
% Compute the projection the most accurate way, according to
[1].
for k=1:size(A,2)
B = B -
A(:,k)*(A(:,k)'*B);
end
% Make sure it's magnitude is less than 1.
theta = asin(min(ones(superiorfloat(A,B)),(norm(B))));
%說明
以下代碼
for k=1:size(A,2)
B = B -
A(:,k)*(A(:,k)'*B);
end
相當(dāng)于如下公式計(jì)算:
B = (eye(size(A,1))-A*A')*B;
這樣可以減少循環(huán)次數(shù),提高計(jì)算速度;
矩陣B的特征值即為夾角正弦值
theta = asin(min(ones(superiorfloat(A,B)),(norm(B))));
superiorfloat() 返回 'single','double'...即A,B的數(shù)據(jù)類型
ones('double') 返回 1
以上代碼說明,夾角正弦值應(yīng)該小于等于1,當(dāng)B的特征值大于1時(shí),取1。
相關(guān)標(biāo)簽搜索:子空間夾角計(jì)算-MATLAB源代碼修正 MatLab培訓(xùn) MatLab培訓(xùn)課程 MatLab在線視頻教程 MatLab技術(shù)學(xué)習(xí)教程 MatLab軟件教程 MatLab資料下載 MatLab代做 MatLab基礎(chǔ)知識(shí) Fluent、CFX流體分析 HFSS電磁分析 Ansys培訓(xùn) Abaqus培訓(xùn)