orthonormalize_qr.m

function Xon = orthonormalize_qr(X,A,epsilon)
% function Xon = orthonormalize_qr(X[,A,epsilon])
%
% orthonormalization of vectors in columns of 
% matrix X to Xon by qr-decomposition of vectors and 
% the inner product matrix A can be specified additionally, i.e.
% <x,x> = x' * A * x. Then the Cholesky factoriziation of A is 
% performed A = R_M' * R_M
%
% Then a QR-decomposition of R_M * X is performed which yields the
% desired new vectors

% Bernard Haasdonk 22.8.2008

if nargin < 3
  epsilon = 1e-7; % => for nonlinear evolution this value is required.
end;

if isempty(X)
  Xon = zeros(size(X));
  return;
end;

if nargin<2
  A = 1;
end;

% check on identity of vectors (can happen, that numerics does not detect
% this afterwards due to rounding errors!!)

for i=1:(size(X,2)-1)
  for j=(i+1):size(X,2)
    if isequal(X(:,i),X(:,j))
      X(:,j) = 0;
    end;
  end;
end;

Xon = X;

% incomplete cholesky of inner-product matrix:
if matlab_version> 7.9
  opts.shape = 'upper';
  opts.type = 'ict';
  R_M = ichol(sparse(A), opts);
else
     R_M = cholinc(sparse(A),'inf');
end;

% qr decomposition of R_M * X with permutation indices E
[Q, R, E]= qr(R_M * Xon,0 );

% sort such that linear independent first and Q * R = R_M * Xon
Xon = Xon(:,E);

% search nonvanishing diagonal entries of R
ind = find(abs(diag(R))>epsilon);
Xon = Xon(:,ind);
Rind = R(ind,ind); 
Xon = Xon(:,ind) / Rind;

% resort vectors such that original order is recovered most:
K = abs(X'* A * Xon);
% search in each row maximum entry
ind = [];
for i = 1:size(Xon,2)
  [max_val, max_i] = max(K(i,:)); 
  if length(max_i)>1
    error('found multiple matching vectors!');
  end;
  if max_val > 0
    ind = [ind, max_i];
  end;
  K(:,max_i) = 0; % "delete" from further search
end;
all_ind = ones(1,size(Xon,2));
all_ind(ind) = 0;
rest_ind = find(all_ind);
% resort:
if ~isempty(rest_ind)
  Xon = [Xon(:,ind),Xon(:,rest_ind)];
else
  Xon = Xon(:,ind);  
end;

%disp('check orthonormality!!');
%keyboard;
%K = Xon'*A*Xon;
%if max(max(abs(eye(size(Xon,2))-K)))>1e-3
%  K
%  error(['check orthonormalization!! Non orthogonal vectors' ...%
%        ' generated']
%end;

% eliminate zero-columns
%nsqr = sum(Xon.^2);
%i = find(nsqr > 0.1);
%Xon = Xon(:,i);

%resort zero-columns to back
%nsqr = sum(Xon.^2);
%[s,i] = sort(-nsqr);
%Xon = Xon(:,i);
%| \docupdate