A.15  Gram-Schmidt orthonormalization procedure

The Gram-Schmidt is an automatic procedure to obtain a set of $N$ orthonormal vectors $\mathbf{\text{y}}$ from an input set $\{{\mathbf{\text{x}}}_1,{\mathbf{\text{x}}}_2,\dots ,{\mathbf{\text{x}}}_M\}$ composed by $M$ vectors. If the $M$ vectors are linearly independent, then $N=M$. If there is linear dependency among the vectors, then $N<M$. Listing A.4 illustrates the procedure.

In summary, the first basis function is ${\mathbf{\text{y}}}_1={\mathbf{\text{x}}}_1∕\parallel {\mathbf{\text{x}}}_1\parallel$. Because ${\mathbf{\text{y}}}_1$ is the normalized first input vector, it can properly represent ${\mathbf{\text{x}}}_1$. The next step is to iteratively (in a loop) add new basis functions to represent the remaining vectors. For example, if ${\mathbf{\text{x}}}_2$ were colinear to ${\mathbf{\text{y}}}_1$, it would not be necessary to enlarge the basis set (${\mathbf{\text{y}}}_1$ could represent both ${\mathbf{\text{x}}}_1$ and ${\mathbf{\text{x}}}_2$). If that is not the case, ${\mathbf{\text{y}}}_2$ cannot be simply ${\mathbf{\text{y}}}_2={\mathbf{\text{x}}}_2∕\parallel {\mathbf{\text{x}}}_2\parallel$ because, eventually, ${\mathbf{\text{y}}}_1$ and ${\mathbf{\text{y}}}_2$ can be linearly dependent. Therefore, first ${\mathbf{\text{x}}}_2$ is projected in ${\mathbf{\text{y}}}_1$ and ${\mathbf{\text{y}}}_2$ is the (normalized) error vector of this projection. This process is repeated. For example, when obtaining ${\mathbf{\text{y}}}_i$ to represent a vector ${\mathbf{\text{x}}}_j$, one first takes in account the projection of ${\mathbf{\text{x}}}_j$ into all previously selected basis $1,2,\dots ,i-1$. If the error is not numerically negligible (given a tolerance), this suggests that ${\mathbf{\text{x}}}_j$ does not reside in the space spanned by the current set of basis functions and the normalized projection error is added to this set.

1function [Ah,A]=ak_gram_schmidt(x,tol) 
2% function [Ah,A]=ak_gram_schmidt(x,tol) 
3%Gram-Schmidt orthonormalization procedure (x is real) 
4%Inputs:  x - real matrix: each row is an input vector 
5%         tol - tolerance for stopping procedure (and SVD) 
6%Outputs: Ah  - direct matrix, which is the Hermitian (transpose in 
7%               this case) of A 
8%         A   - inverse matrix with basis functions as columns 
9%Example of usage: 
10% x=[1 1; 0 2; 1 0; 1 4] 
11% [Ah,A]=ak_gram_schmidt(x,1e-12) %basis functions are columns of A 
12% testVector=[10; 10]; %test: vector in same direction as first basis 
13% testCoefficients=Ah*testVector %result is [14.1421; 0] 
14 
15%% Notes about the code: 
16%1) The inner product sum(y(m,:).*x(k,:)) is alternatively 
17%   calculated as y(m,:)*x(k,:)' (row vectors) 
18%2) In general, a projection is projectionOverBasis=p_xy= 
19%             ( <y,x> / ||y||^2 ) * y = 
20%             ((y(m,:)*x(k,:)')/(y(m,:)*y(m,:)'))*y(m,:); 
21%but in the code y(m,:)*y(m,:)'=1 and the denominator is 1 
22if nargin<2 
23    tol=min(max(size(x)')*eps(max(x))); %find tolerance 
24end 
25[m,dimension]=size(x); 
26if dimension>m 
27    warning(['Input dimension larger than number of vectors, ' ... 
28        'which is ok if you are aware about']); 
29end 
30N=rank(x,tol); %note: rank is slow because it uses SVD 
31y=zeros(N,dimension); %pre-allocate space 
32%pick first vector and normalize it 
33y(1,:)=x(1,:)/sqrt(sum(x(1,:).^2)); 
34numBasis = 1; 
35for k=2:m 
36    errorVector=x(k,:); %error (or target) vector 
37    for m=1:numBasis%go over all previously selected basis 
38        %p_xy = <x,y> * y; %recall, in this case: ||y||=1 
39        projectionOverBasis = ((y(m,:)*x(k,:)'))*y(m,:); 
40        %update target: 
41        errorVector = errorVector - projectionOverBasis; 
42    end 
43    magErrorVector = sqrt(sum(errorVector.^2)); 
44    if ( magErrorVector > tol ) 
45        %keep the new vector in basis set 
46        numBasis = numBasis + 1; 
47        y(numBasis,:)= errorVector / magErrorVector; 
48        if (numBasis >= N) 
49            break; %Abort. Reached final number of vectors 
50        end 
51    end 
52end 
53A=transpose(y);  %basis functions are the columns of inverse matrix A 
54Ah=A'; %the direct matrix transform is the Hermitian of A
  

The result of the Gram-Schmidt procedure depends on the order of the input vectors. Rearranging these vectors may produce a different (still valid) solution. An example of Gram-Schmidt execution can be found in Application 2.1.

The following section presents another procedure, which has similarities to the Gram-Schmidt and very interesting properties.