1 function onvec = model_orthonormalize_qr(model, model_data, vec,epsilon)
2 %
function onvec = model_orthonormalize_qr(model, model_data, vec[,epsilon])
4 % orthonormalization of vectors in columns of
5 % matrix X to Xon by qr-decomposition of vectors and
6 % the inner product matrix A can be specified additionally, i.e.
7 % <x,x> = x' * A * x. Then the Cholesky factoriziation of A is
8 % performed A = R_M' * R_M
10 % Then a QR-decomposition of R_M * X is performed which yields the
12 % additionally, a reordering is performed, such that the onb is at
13 % most matched with original vectors
16 epsilon = 1e-7; % =>
for nonlinear evolution
this value is required.
18 %epsilon = 1e-10; =>
for linear evolution
this value was OK.
19 %epsilon = 1e-12; -> This causes error by returning non-orthogonal
20 % vectors (e.g. two identical vectors)
21 % So take epsilon larger than 1e-12!
24 onvec = zeros(size(vec));
28 A = model.get_inner_product_matrix(model_data);
30 % incomplete cholesky of inner-product matrix:
33 R_M = ichol(sparse(A), opts);
35 % qr decomposition of R_M * X with permutation indices E
36 % [Q,R,E] = QR(B,0) produces an "economy size" decomposition in which E
37 % is a permutation vector, so that B(:,E) = Q*R.
39 [Q, R, E]= qr(R_M * vec,0 );
40 %[Q, R]= qr(R_M * vec,0 );
44 % sort such that linear independent first and Q * R = R_M * Xon
47 % search nonvanishing diagonal entries of R
48 ind = find(abs(diag(R))>epsilon);
51 onvec = vec2(:,ind) / Rind;
53 % eliminate zero-columns
54 %nsqr = sum(onvec.^2);
55 %i = find(nsqr > 0.1);
58 %resort zero-columns to back
59 %nsqr = sum(onvec.^2);
63 % resort vectors such that original order is recovered most:
64 K = abs(vec'* A * onvec);
65 % search in each row maximum entry
68 [max_val, max_i] = max(K(i,:));
70 error('found multiple matching vectors!');
75 K(:,max_i) = 0; % "delete" from further search
77 all_ind = ones(1,size(onvec,2));
79 rest_ind = find(all_ind);
82 onvec = [onvec(:,ind),onvec(:,rest_ind)];
88 K = model.get_inner_product_matrix(model_data);
89 if max(max(abs(eye(size(onvec,2))-onvec'*K*onvec)))>1e-3
91 max(max(abs(eye(size(onvec,2))-onvec'*K*onvec)))
92 %error(['check orthonormalization!! Non orthogonal vectors' ...