get_orthonormal(m,n)

function answer = get_orthonormal(m,n)
% Produces an m x n set of orthonormal vectors, 
% (thats n vectors, each of length m)
% 
% Inputs should be two scalars, m and n, where n is smaller than 
% or equal to m.
%
% Example: >> get_orthonormal(5,4)
%
% ans =
%     0.1503   -0.0884   -0.0530    0.8839
%    -0.4370   -0.7322   -0.1961   -0.2207
%    -0.3539    0.3098    0.7467   -0.0890
%     0.7890   -0.1023    0.0798   -0.3701
%    -0.1968    0.5913   -0.6283   -0.1585



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% CHECK USER INPUT
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


if ( (nargin==2) && (m>n) && (isnumeric(m)*isnumeric(n)) )
    
elseif ( nargin==1 && isnumeric(m) && length(m)==1 )
    
    n=m;
    
else
   error('Incorrect Inputs. Please read help text in m-file.')
end




% to get n orthogonal vectors (each of size m), we will first get a larger mxm 
% set of orthogonal vectors, and then just trim the set so it is 
% of size mxn,
%
% to get an mxm set of orthogonal vectors, 
% we can exploit the fact that the eigenvectors
% from distinct eigenspaces (corresponding to different eigenvalues) of a
% symmetric matrix are orthogonal



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Create the orthonormal vectors
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

count=0;
while (count==0)

    % generate an mxm matrix A, then make a symmetric mxm matrix B
    A=rand(m);
    B=A'*A ;

    % now take the eigenvectors of B, 
    % eigenvectors from different eigenspaces (corresponding to different
    % eigenvalues) will be orthogonal

    % there is a chance that there will be repeat eigenvalues, which would give
    % rise to non-orthogonal eigenvectors (though they will still be independent)
    % we will check for this below
    % if this does happen, we will just start the loop over again 
    % and hope that the next randomly created symmetric matrix will not
    % have repeat eigenvalues,
    % (another approach would be to put the non-orthogonal vectors
    % through the gram-schmidt process and get orthogonal vectors that way)

    % since matlab returns unit length eigenvectors, they will also be
    % orthonormal

    [P,D] = eig(B) ;

    % can double check the orthonormality, by taking the difference between 
    % P'P and I, if it is non-zero, then there is an error (caused by repeat
    % eigenvalues), repeat the loop over again

    if ((P'*P - eye(m))>eps) 
        % error, vectors not orthonormal, repeat the random matrix draw again
        count=0
    else
        % we want the first n of these orthonormal columns
        answer=P(:,1:n) ;
        count=1;
    end


end