dualsparsegeneralfeatureslm2.m

a function inside machine learning

function  [newTrainX, newTestX, subspaceInfo] = dualSparseGeneralFeaturesLM2(trainX, trainY, testX, params)
%A function to compute general features in the dual space, which will not
%require the kernel matrix to be entirely in memory. Designed for datasets
%with a large number of examples. 
%TODO: Rescale alignment based directions 

if (nargin ~= 4)
    disp('Incorrect number of inputs');
    help dualSparseGeneralFeaturesLM2;
    return;
end

useSparse = issparse(trainX); 

if useSparse 
    disp('Using sparse representation'); 
    dataRank = sprank(trainX); 
else 
    dataRank = rank(trainX); 
end 

%Store all the parameters 
chunkSize = params.chunkSize; 
dualSparseMeasureFunction = char(params.dualSparseMeasureFunction); 
T = min(params.iterations, dataRank); 
kernelFunction = char(params.kernelFunctionName); 
kernelParams = params; 
cacheSize = params.cacheSize; 

if isfield(params, 'normalise')
    normaliseFeatures = params.normalise;
else 
    normaliseFeatures = 1; 
end

fprintf('Using chunks of size %d\n', chunkSize); 



numTrainExamples = size(trainX, 1); 
numTestExamples = size(testX, 1); 
numFeatures = size(trainX, 2); 
numTrainChunks = ceil(numTrainExamples/chunkSize);  %Maximum number of chunks
numTestChunks = ceil(numTestExamples/chunkSize); 

alpha = 0.0001; %This number is added to the diagonal of matrices to make them non singular
tol = 10^-3; 

Yj = trainY; 

%Some temp variables
numCols = min(numTrainExamples, cacheSize); 
fullMeasures = zeros(numTrainExamples, 1);
partialMeasures = zeros(numCols, 1);
tempKjCol = zeros(numTrainExamples, 1); 

kernelStatus = 1; 
lastFullKernel = 0; %The last time the full kernel matrix was evaluated 

normSqTau = zeros(T, 1);

%Note that Kj = XjX' 
if ~useSparse 
    Kj = zeros(chunkSize, numTrainExamples);
    b = zeros(numTrainExamples, T); 
    tau = zeros(numTrainExamples, T);
    KTau = zeros(numTrainExamples, T); %tau'K     
    tempX = zeros(chunkSize, numFeatures);
    tempK = zeros(numTrainExamples, chunkSize);
    tempB = zeros(numTrainChunks, 1);
    reducedKTau = zeros(numCols, T); 
else 
    Kj = sparse(chunkSize, numTrainExamples);
    b = sparse(numTrainExamples, T);
    tau = sparse(numTrainExamples, T);
    KTau = sparse(numTrainExamples, T); %tau'K
    tempX = sparse(chunkSize, numFeatures);
    tempK = sparse(numTrainExamples, chunkSize);
    tempB = sparse(numTrainChunks, 1);
    reducedKTau = sparse(numCols, T);
end

%Now deflate and find new dual directions
for i=1:T
    fprintf('Iteration %d\n', i);
    
    %Compute the optimum tau, b for the the kernel
    %matrix
    if kernelStatus == 1
        previousMeasure = -1; 
        
        for j=1:numTrainChunks
            fprintf('Processing kernel chunk %d\n', j);
            startExample = ((j-1)*chunkSize)+1;
            endExample = min(j*chunkSize, numTrainExamples);
            tempChunkSise = endExample-startExample+1;

            if ~useSparse
                sumTauTauK = zeros(numTrainExamples, tempChunkSise);
            else
                sumTauTauK = sparse(numTrainExamples, tempChunkSise);
            end

            rowIndices = startExample:endExample;
            tempX = trainX(rowIndices, :); 
            tempK = feval(kernelFunction, trainX, tempX, kernelParams);
            
            %Need to recompute Ktau for all the missed iterations 
            for k=lastFullKernel:i-2  %i-2 since we have the if statement below 
                KTau(rowIndices, k) = tempK'*tau(:, k);
            end
                        
            %Need to compute KTau
            if i > 1
                KTau(rowIndices, i-1) = tempK'*tau(:, i-1);
            end

            %This part of the kernel matrix with all rows, partial cols
            sumTauTauK = tau(:, 1:i-1)*(diag(1./normSqTau(1:i-1))*KTau(rowIndices, 1:i-1)');
            Kj = tempK - sumTauTauK;

            %We want to keep the columns that have a high measure (e.g.
            %alignment, covariance) 
            [fullMeasures(rowIndices), bs] = feval(dualSparseMeasureFunction, tempK, Kj, trainY, Yj, rowIndices'); 
            [maxMeasure, k] = max(abs(fullMeasures(rowIndices))); 
            
            if maxMeasure > previousMeasure 
                previousMeasure = maxMeasure; 
                tempKjCol = Kj(:, k); 
                tempbi = bs(k);
            end
            
            %At the last chunk figure out which columns are the best
            if j == numTrainChunks
                fullMeasuresAndIndices = [fullMeasures , (1:numTrainExamples)'];
                fullMeasuresAndIndices = sortrows(fullMeasuresAndIndices);
                [maxMeasure, k] = max(fullMeasuresAndIndices(:, 1));
                maxMeasureIndex = fullMeasuresAndIndices(k, 2);
                minMeasure = fullMeasuresAndIndices(1, 1);
                %minMeasure = 0.02;  %Force recomputation 
                
                b(maxMeasureIndex, i) = tempbi;
                tau(:, i) = tempKjCol*tempbi;
                
                relevantColumns = unique((fullMeasuresAndIndices(end-numCols+1:end, 2) > 0).*fullMeasuresAndIndices(end-numCols+1:end, 2));
                
                if relevantColumns(1) == 0 %The first elment might be zero 
                    relevantColumns = relevantColumns(2:end);  
                end
                
                lastFullKernel = i; 
                kernelStatus = -1; %means we need to compute a new kernel using only n best columns 
            end
            clear sumTauTauK Kj; 
        end
    else
        if kernelStatus == -1
            disp('Using reduced kernel matrix'); 
            tempX = trainX(relevantColumns, :); 
            tempK = feval(kernelFunction, trainX, tempX, kernelParams);
            
            %Need to compute the reduced KTau
            for j=1:i-2  %i-2 since we have the if statement below 
                reducedKTau(:, j) = tempK'*tau(:, j);
            end
            
            sumTauTauK = tau(:, 1:i-2)*(diag(1./normSqTau(1:i-2))*reducedKTau(:, 1:i-2)');
            Kj = tempK - sumTauTauK;
            
            kernelStatus = 0;  %Means we have computed the new reduced kernel 
        end
        
        %Need to compute KTau
        if i > 1
            reducedKTau(:, i-1) = tempK'*tau(:, i-1);
        end
        
        %Deflate Kj 
        Kj = Kj - tau(:, i-1)*(tau(:, i-1)'*Kj)/normSqTau(i-1); 
        
        %Compute the measure of the last tau here
        partialMeasures = feval(dualSparseMeasureFunction, tempK, Kj, trainY, Yj, relevantColumns); 
        [maxMeasure, maxMeasureIndex] = max(partialMeasures); 
        
        exampleIndex = relevantColumns(maxMeasureIndex); 
        
        %Need to compute the scaling factor a 
        tempExample = trainX(exampleIndex, :); 
        a = sqrt(abs(feval(kernelFunction, tempExample, tempExample, kernelParams)));         
        b(exampleIndex, i) = 1/a;
        tau(:, i) = Kj(:, maxMeasureIndex)/a;

        if maxMeasure < minMeasure
            fprintf('Measure has dropped to %f, which is lower than minimum of %f\n', full(maxMeasure), full(minMeasure)); 
            kernelStatus = 1;
        end
               
    end
    
    %Matlab memory 
    clear sumTauTauK tempX;
    pack; 

    normSqTau(i) = tau(:, i)'*tau(:, i); 
    %Deflate Yj
    Yj = Yj -  tau(:, i)*((tau(:, i)'*Yj)/normSqTau(i));
end

clear Kj tempK; 

%We need to compute trainK * b. Since b is sparse, our life is slightly
%easier. Also need testTrainK * b; 
[exampleIndices, nzElements] = findNonZeroElements(b);
diagBElements = diag(nzElements);

tempTrainK = feval(kernelFunction, trainX, trainX(exampleIndices, :), kernelParams);
trainKb = tempTrainK*diagBElements;

tempTestTrainK = feval(kernelFunction, testX, trainX(exampleIndices, :), kernelParams);
testTrainKb = tempTestTrainK*diagBElements;

%Z should be sparse hopefully - otherwise we're stuffed 
Q = (tau'*trainKb + eye(T)*alpha);
Z = (diag(diag(tau'*tau)) + eye(T)*alpha)\Q; 

newTrainX = trainKb/Z;
newTestX = testTrainKb/Z;

if normaliseFeatures == 1
    disp('Features are normalised');
    [newTrainX, newTestX] = normalise(newTrainX, newTestX);
end

subspaceInfo = struct; 
subspaceInfo.b = sparse(b); 
%subspaceInfo.tau = tau;