ITE / code / shared / embedded / ITL / incompleteCholeskyMulti.m

% The following function computes the Incomplete Cholesky decomposition of
% the Gram matrix generated by the random VECTOR samples
%
% Input:   "samples" should be COLUMN vectors of length (nxm)
%				'kernelSize' is a scalar for the kernel size.
%               "precision" is a scalar for the precision
%
% Output: 'choleskyMarix' contains nxd Cholesky matrix
%
% Default:  precision = 10^-6.
%
% The code uses incompleteCholesky when m is 1.
%
% Author: Sohan Seth (sohan@cnel.ufl.edu)	Date: 01.06.2009

function choleskyMatrix  = incompleteCholeskyMulti(samples,kernelSize,precision)

if nargin == 2
    precision = 10^-6;
end

[n m] = size(samples);
if m == 1
    choleskyMatrix  = incompleteCholesky(samples,kernelSize,precision);
    return;
end

twoSigmaSquare = 2*kernelSize^2;
permuteVector = 1:n;
partialTraceCholeskyMatrix = zeros(1,n);
% choleskyMatrix = zeros(n, n/100);
choleskyMatrix(:,1) = exp(-sum((samples - repmat(samples(1,:),n,1)).^2,2)/twoSigmaSquare);
diagCholeskyMatrix = zeros(1, n);
indexOfDiagonalEntries = sub2ind([n,n], 1:n, 1:n);

for countCol=1:n

    if countCol == 1
        diagCholeskyMatrix(countCol:n) = (exp(-sum((samples(permuteVector(countCol:n),:) - samples(permuteVector(countCol:n),:)).^2,2)/twoSigmaSquare))';
    else
        diagCholeskyMatrix(countCol:n) = (exp(-sum((samples(permuteVector(countCol:n),:) - samples(permuteVector(countCol:n),:)).^2,2)/twoSigmaSquare))' ...
            - sum(choleskyMatrix(countCol:n,1:countCol-1).^2,2)';
    end

    partialTraceCholeskyMatrix(countCol) = sum(diagCholeskyMatrix(countCol:n));

    if  partialTraceCholeskyMatrix(countCol) <= 0
        warning(['Negative diagonal entry', num2str(partialTraceCholeskyMatrix(countCol))])
    end

    if  partialTraceCholeskyMatrix(countCol) <= precision
        break
    else

        [bestElementValue bestElementIndex] = max(diagCholeskyMatrix(countCol:n));

        bestElementIndex = bestElementIndex+countCol-1;

        temp1 = permuteVector(countCol);
        permuteVector(countCol) = permuteVector(bestElementIndex);
        permuteVector(bestElementIndex) = temp1;

        temp2 = choleskyMatrix(countCol, 1:countCol-1);
        choleskyMatrix(countCol, 1:countCol-1) = choleskyMatrix(bestElementIndex, 1:countCol-1);
        choleskyMatrix(bestElementIndex, 1:countCol-1) = temp2;

        bestElementValue = sqrt(bestElementValue);

        choleskyMatrix(countCol,countCol) = bestElementValue;

        choleskyMatrix(countCol+1:n,countCol) = (exp(-sum((repmat(samples(permuteVector(countCol),:),n-countCol,1) - samples(permuteVector(countCol+1:n),:)).^2,2)/twoSigmaSquare) ...
            - choleskyMatrix(countCol+1:n,1:countCol-1)*choleskyMatrix(countCol,1:countCol-1)')/bestElementValue;

    end
end

if countCol ~= n
    choleskyMatrix = choleskyMatrix(:,1:countCol-1);
end

temp = [permuteVector',choleskyMatrix];
temp = sortrows(temp,1);
choleskyMatrix = temp(:,2:end);
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.