## 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:

• Build the covariance matrix in block form given the value of the integral. Otherwise also calculate the value of the integral
• Find the Cholesky in block form. Take care of negative EV.
• Rebuild the matrix

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
```