compute projection matrix from grassmann coordinates THIS IS NO USER FUNCTION
0001 function projectionMatrix = computeProjectionFromGrass(x, grassStruct) 0002 % compute projection matrix from grassmann coordinates 0003 % 0004 % THIS IS NO USER FUNCTION 0005 0006 % The elk-library: convex geometry applied to crystallization modeling. 0007 % Copyright (C) 2013 Alexander Reinhold 0008 % 0009 % This program is free software: you can redistribute it and/or modify it 0010 % under the terms of the GNU General Public License as published by the 0011 % Free Software Foundation, either version 3 of the License, or (at your 0012 % option) any later version. 0013 % 0014 % This program is distributed in the hope that it will be useful, but 0015 % WITHOUT ANY WARRANTY; without even the implied warranty of 0016 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0017 % General Public License for more details. 0018 % 0019 % You should have received a copy of the GNU General Public License along 0020 % with this program. If not, see <http://www.gnu.org/licenses/>. 0021 0022 %% projection matrix in scaled-PCA space 0023 % set tangent space in grassmann manifold (see documentation for details) 0024 tanMatrix = grassStruct.baseNull * reshape(x, ... 0025 grassStruct.nDim - grassStruct.targetDim, ... 0026 grassStruct.targetDim); 0027 [U S V] = svd(tanMatrix,'econ'); 0028 baseMatrix = grassStruct.baseRange * V * diag(cos(diag(S))) * V'... 0029 + U * diag(sin(diag(S))) * V'; 0030 0031 % we check orthonormality by comparing (baseMatrix'*baseMatrix) which should 0032 % result in a unity matrix by summing up all squarred matrix elements of 0033 % the difference: 0034 if abs(norm(baseMatrix'*baseMatrix - eye(size(baseMatrix,2)),'fro')) > grassStruct.zerotol 0035 error('elk:decomposition:numericalError', ['The basis obtained from ' ... 0036 'the Grassmann manifold is not orthogonal. This should not happen.']); 0037 % This would be the code to obtain an orthonormal baseMatrix again: 0038 baseMatrix = orth(baseMatrix); 0039 end 0040 0041 % transform basis to real-data space 0042 baseMatrix = grassStruct.mapToRealSpace * baseMatrix; 0043 % after this coordinate transformation, the base Matrix is NOT orthonormal 0044 % anymore and the projection must be computed as: 0045 projectionMatrix = baseMatrix*inv(baseMatrix'*baseMatrix)*baseMatrix';