A.M. Yaglom, in his "Correlation Theory of Stationary and Related Random Functions. Volume I: Basic Results", ISBN: 0-387-96268-9, presents the damped cosine function (equation (2.116) of page 127) $$ B(\tau)=C\,e^{-\alpha\tau}\cos(\omega_0\tau-\psi). $$ Consider only $\tau\geq0$ or, as it is stated originally in the reference, $|\tau|$. It is said that such function is a correlation function (it can be understood also as covariance function) whenever $C>0$, $\alpha>0$, $\omega_0>0$, and (equation (2.117) of page 128, and statement of page 129) $$ |\psi|\leq\arctan\left(\frac{\alpha}{\omega_0}\right). $$
Well, if $B(\tau)$ is a covariance function, a covariance matrix derived from it must be positive definite. In page 129 it is presented the function for the particular values $$ R(\tau)=e^{-0.06\tau}\cos(0.21\tau-0.08), $$ where $C=1$, $\alpha=0.06$, $\omega_0=0.21$, and $\psi=0.08$. As $C>0$, $\alpha>0$, $\omega_0>0$, and $|\psi|\leq\arctan(\alpha/\omega_0)\rightarrow0.08\leq0.2783$, the function is a covariance function.
Consider now the below Octave/Matlab code, which produces the images at the bottom. If we try to do the Cholesky decomposition of the matrix we can see that it is not positive definite.
I can't understand why, but also I'm pretty sure I'm misunderstanding something. I'm working on least squares collocation, where covariance matrices derived from covariance functions are needed, and must be positive definite. Can someone explain me what am I misunderstanding?
%Function parameters in A.M. Yaglom, p. 129 C = 1.0; alpha = 0.06; omega = 0.21; p = 0.08; %Function according to Yaglom: R(t) = C*exp(-alpha*t).*cos(omega*t-p) val1 = tan(p); val2 = alpha/omega; fprintf(1,'For positive-definiteness must be tan(p)<=alpha/omega:\n'); fprintf(1,'\ttan(%f)<=%f/%f -> %f<=%f\n',p,alpha,omega,val1,val2); if val1<=val2 fprintf(1,'\tFunction IS positive-definite\n'); else fprintf(1,'\tFunction NOT positive-definite\n'); end %Function plot t = 0.0:0.1:40.0; R = C*exp(-alpha*t).*cos(omega*t-p); figure(1); plot(t,R,'*'); title('Yaglom'); xlabel('t'); ylabel('R'); %Evaluation points forming the covariance matrix x = 1.0:25.0; y = 1.0:25.0; [xx,yy] = meshgrid(x,y); xx = xx(:)'; yy = yy(:)'; n = length(xx); %Yaglom covariance matrix covY = zeros(n); covS = zeros(n); for i=1:n %Distances t = sqrt((xx(i)-xx(i:end)).^2+(yy(i)-yy(i:end)).^2); %Covariances cY = C*exp(-alpha*t).*cos(omega*t-p); %Assign to the matrix covY(i,i:end) = cY; covY(i:end,i) = cY'; end %Plot figure(2); imagesc(covY); colorbar; axis('equal'); title('Yaglom'); %Try Cholesky factorization for checking positive-definiteness try chol(covY); fprintf(1,'Yaglom covariance matrix IS positive-definite\n'); catch fprintf(1,'Yaglom covariance matrix NOT positive-definite\n'); end The results of the script are:
and
