Contents

Find the Cholesky decomposition of the covariance matrix

Commented out: old cold to build the covariance if the j values are not passed to the function. Uses residue calculus and helper function.

Steps:

Christel Hohenegger

Last updated: 05/01/2016

function matC=CovMat(parameters,jInt)

%function matC=CovMat(parameters,kk,diagindex)

numTsteps = parameters.numTsteps;
numBlocks = min(2^10,numTsteps);
count = numTsteps/numBlocks;

Get the results from residue calculus (commented out)

Needed poles, residues, t-vector

%{
residueCalc = ResCalc(parameters,kk);
tVec = residueCalc.tVec;
resVec = residueCalc.resVec;
omegaVec = residueCalc.omegaVec;
%}

Build the covariance matrix in block form

Commented out: old code if the values of the j integral is not passed to the function. Caclulated instead using residue calculus.

matA = zeros(numBlocks,numBlocks,count,count);
jVal = zeros(numBlocks+1,count);
%{
numKernels = parameters.numKernels;
dtND = parameters.dtND;
tND=dtND*(0:numTsteps);
if diagindex == 1
    cGLE = 1/(2*pi);
else
    cGLE = 1/(2*pi*sqrt(2));
end
aGLE = kk^2;
bGLE = parameters.beta*kk^2;
%}
for iindex=1:count
    jVal(:,iindex) = jInt(1+(iindex-1)*numBlocks:1+iindex*numBlocks);
%{
    tVal = tND(1+(iindex-1)*numBlocks:1+iindex*numBlocks);
    for jindex=1:length(tVal)
        jVal(jindex,iindex) = cGLE^2*(tVal(jindex)*tVec(numKernels+1)...
            /(aGLE*tVec(numKernels+1)+1i*bGLE/numKernels*tVec(numKernels))...
            +sum(imag(resVec./omegaVec.^2 ...
            .*(exp(1i*tVal(jindex)*omegaVec)-1))));
    end
%}
end
%jVal(1,1)=0;
for iindex=1:count
    for jindex=1:iindex
        matR3 = repmat(jVal(2:end,jindex)',numBlocks,1);
        matR1 = repmat(jVal(2:end,iindex),1,numBlocks);
        if iindex==jindex
            matR2 = toeplitz(jVal(1:end-1,1));
        else
            matR2 = toeplitz(jVal(1:end-1,iindex-jindex+1),...
                flipud(jVal(2:end,iindex-jindex))');
        end
        matA(:,:,iindex,jindex) = matR1-matR2+matR3;
    end
end

Find the Cholesky

Find negative eigenvalues and replace them by a positive tolerance. This is done by decomposing the matrix, finding the appropriate value and rebuilding the matrix (diagaonal decomposition).

matL = zeros(numBlocks,numBlocks,count,count);
tol = 1e-8;
for iindex=1:count
    for jindex=1:iindex
        if iindex==jindex
            temp = 0;
            for kindex=1:iindex-1
                temp = temp+squeeze(matL(:,:,iindex,kindex))...
                    *squeeze(matL(:,:,iindex,kindex))';
            end
            try
                matL(:,:,iindex,jindex) = chol(squeeze(...
                    matA(:,:,iindex,jindex))-temp,'lower');
            catch
                [matV,matD]=eig(squeeze(matA(:,:,iindex,jindex))...
                    -temp,'nobalance');
                diag(matD);
                vecD = diag(matD);
                fprintf('ERROR: Matrix is not positive definite.\n')
                fprintf('Largest negative EV %6.4e .\n',min(vecD))
                vecD(vecD<0)=tol;
                matC = matV*diag(vecD)/matV;
                matL(:,:,iindex,jindex) = chol(matC,'lower');
            end
        else
            temp = 0;
            for kindex=1:jindex-1
                temp = temp+squeeze(matL(:,:,iindex,kindex))...
                    *squeeze(matL(:,:,jindex,kindex))';
            end
            matL(:,:,iindex,jindex) = (squeeze(matA(:,:,iindex,jindex))...
                -temp)/(squeeze(matL(:,:,jindex,jindex))');
        end
    end
end

Rebuild the matrix

matC = zeros(numTsteps,numTsteps);
for iindex=1:count
    for jindex=1:iindex
        matC(1+(iindex-1)*numBlocks:iindex*numBlocks,...
            1+(jindex-1)*numBlocks:jindex*numBlocks) = ...
            matL(:,:,iindex,jindex);
    end
end

Published with MATLAB® R2015b