model_orthonormalize_qr.m

function onvec = model_orthonormalize_qr(model, model_data, vec,epsilon)
% function onvec = model_orthonormalize_qr(model, model_data, vec[,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
% additionally, a reordering is performed, such that the onb is at
% most matched with original vectors

if nargin < 4
  epsilon = 1e-7; % => for nonlinear evolution this value is required.
end;
%epsilon = 1e-10; => for linear evolution this value was OK.
%epsilon = 1e-12; -> This causes error by returning non-orthogonal 
%                    vectors (e.g. two identical vectors)
%                    So take epsilon larger than 1e-12!

if isempty(vec)
  onvec = zeros(size(vec));
  return;
end;

A = model.get_inner_product_matrix(model_data);

% 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(B,0) produces an "economy size" decomposition in which E
%    is a permutation vector, so that B(:,E) = Q*R.

[Q, R, E]= qr(R_M * vec,0 );
%[Q, R]= qr(R_M * vec,0 );

%keyboard;

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

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

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

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

% resort vectors such that original order is recovered most:
K = abs(vec'* A * onvec);
% search in each row maximum entry
ind = [];
for i = 1:size(vec,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(onvec,2));
all_ind(ind) = 0;
rest_ind = find(all_ind);
% resort:
if ~isempty(rest_ind)
  onvec = [onvec(:,ind),onvec(:,rest_ind)];
else
  onvec = onvec(:,ind);  
end;

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