It seems I am there and have my non-positive covariance matrix repaired.
However, I just wanted to ask the people here whether you think I did
everything correct?
I am doing it the following way: Qf is my non-positive covariance matrix,
via Matlab's chol() I test whether it really is non-positive. I perfom
'modified choleski decomposition', see code at bottom of email, to calculate
D with all off-diagonal elements equals zero. However, I replace their
main-diagonal negative values by zero and REPEAT the whole thing until Qf is
positive and p=0.
However, the QUESTION is whether I can do it repeatedly like that, Grewal
and Andrews dont write anything about repeated execution. I tested the
initial and final covariance matrix and the maximum error,
max(max(Qf-save_Qf)), I can observe is within 10^-18.
The second QUESTION I have is why C=chol(Qf) return something
different from C=utchol(Qf), see code at bottom? They are both upper
triangular,
i think, but their values are clearly different. Is that just a numerical
issue?
% -----------------------
% repair covariance matrix
[C,p] = chol(Qf);
while p>0 % test if non-positive
[U,D]=udcomp(Qf); % from Grewal and Andrews, p. 222
%[U,D]=schur(Qf); % doesn't work because D has off-diagonal
non-zero terms
for i=1:length(D)
if D(i,i)<0 D(i,i)=0; end % set negative entries to zero, see
Grewal and Andrews, p. 293
end
Qf = U*D*U';
[C,p] = chol(Qf);
end
C=utchol(Qf); % from Grewal and Andrews
Thanks for your help in advance,
Tobias
-------------------
Cod for udcomp:
function [U,D]=udcomp(M)
m=length(M);
for j=m:-1:1
for i=j:-1:1
sigma=M(i,j);
for k=j+1:m
sigma=sigma-U(i,k)*D(k,k)*U(j,k);
end
if i==j
D(j,j)=sigma;
U(j,j)=1;
else
U(i,j)=sigma/D(j,j);
end
end
end
-------------------
Cod for utchol:
function C = utchol(P)
%
% for P symmetric and positive definite,
% computes upper triangular C such that
% C*C' = P
%
[n,m] = size(P);
if (n-m) error('non-square argument'); end;
for j=m:-1:1,
for i=j:-1:1,
sigma = P(i,j);
for k=j+1:m,
sigma = sigma - C(i,k)*C(j,k);
end;
C(j,i) = 0;
if (i==j)
C(i,j) = sqrt(max([0,sigma]));
elseif (C(j,j) == 0)
C(i,j) = 0;
else
C(i,j) = sigma/C(j,j);
end;
end;
end;
------------------
PS: I need C for the generation of discrete-time process noise, that is,
w=C*randn(n,1). Which C should I use, I would think the one from utchol(),
but
why?
.
.
=================================================================
Instructions for joining and leaving this list, remarks about the
problem of INAPPROPRIATE MESSAGES, and archives are available at:
. http://jse.stat.ncsu.edu/ .
=================================================================
