dualsparsegeneralfeatureslm.m

a function inside machine learning

function  [newTrainX, newTestX, subspaceInfo] = dualSparseGeneralFeaturesLM(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. 

if (nargin ~= 4)
    disp('Incorrect number of inputs');
    help dualSparseGeneralFeaturesLM;
    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; 
sparseDualFeatureDirection = char(params.dualFeatureDirection); 
T = min(params.iterations, dataRank); 
kernelFunction = char(params.kernelFunctionName); 
kernelParams = params; 

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

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

%Our dual projection directions and scalings
b = zeros(numTrainExamples, T); 
tau = zeros(numTrainExamples, T); 
KTau = zeros(numTrainExamples, T); %tau'K
normSqTau = zeros(T, 1); 

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

Yj = trainY; 

%Note that Kj = XjX' 
if ~useSparse 
    Kj = zeros(chunkSize, numTrainExamples);
    tempX = zeros(chunkSize, numFeatures);
    tempK = zeros(numTrainExamples, chunkSize);
    tempB = zeros(numTrainChunks, 1);
else 
    Kj = sparse(chunkSize, numTrainExamples);
    tempX = sparse(chunkSize, numFeatures);
    tempK = sparse(numTrainExamples, chunkSize);
    tempB = sparse(numTrainChunks, 1);
end

%Now deflate and find new dual directions
for i=1:T
    fprintf('Iteration %d\n', i);

    %Compute the optimum tau, b, and a for the the kernel
    %matrix
    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 compute KTau
        if i > 1
            KTau(rowIndices, i-1) = tempK'*tau(:, i-1);
        end 
        
        sumTauTauK = tau(:, 1:i-1)*(diag(1./normSqTau(1:i-1))*KTau(rowIndices, 1:i-1)');

        %This part of the kernel matrix with all rows, partial cols
        Kj = tempK - sumTauTauK; 
        
        [b(:, i), tau(:, i)] = feval(sparseDualFeatureDirection, tempK, Kj, trainY, Yj, tau(:, i), b(:, i), startExample);
        
        clear sumTauTauK;
    end

    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)))\Q;  %inv(tau'tau) * 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; 
%subspaceInfo.Z = b/Z;  %b * inv(Z)