function final = testfcn
% data1 = load('E'); % Not used
data2 = load('R');
data3 = load('P');
% tem1 = data1.tem1;
tem2 = data2.tem2;
tem4 = data3.tem4;
sparse = zeros(14596, 2209);
final = zeros(144, 2209, 'uint8');
tic
% Normalize the dictionary - outside the ii loop:
D = zeros(size(tem2));
for index = 1:size(D, 2)
D(:,index) = tem2(:,index) ./ norm(tem2(:,index));
% v = dictionary(:, index);
% dictionary(:,index) = v ./ sqrt(v' * v); % Not exactly the same!!!
end
warning('OFF', 'MATLAB:nearlySingularMatrix');
parfor ii = 1:15 % size(tem4,2)
t = 5;
s = tem4(:, ii);
r = s;
atoms = zeros(size(D));
coefs = zeros(size(D, 2),1);
dictionary = D;
index = [];
while t > 1e-15 && any(dictionary(:))
% inner_product = dictionary' * r; % Dot Product
% [m, ind] = max(abs(inner_product));
ind = GetMaxIndex(r, dictionary);
index = [index, ind]; % Ugly: growing vector!
atoms(:, ind) = dictionary(:, ind); % Select maximum atom
dictionary(:, ind) = 0;
x = (atoms(:, index).' * atoms(:, index)) \ (atoms(:, index).' * r);
r = r - atoms(:, index) * x;
% v = atoms(:, index); % Not exactly the same!!!
% x = (v.' * v) \ (v.' * r);
% r = r - v * x;
t = norm(r);
end
coefs(index) = x; % Outside the loop
sparse(:,ii) = coefs;
sig = D * coefs;
final(:, ii) = (((max(s)) - min(s)) / (max(sig) - min(sig))) * ...
(sig - min(sig)) + min(s);
end
warning('ON', 'MATLAB:nearlySingularMatrix');
toc
And GetMaxIndex is a C-mex function:
// Author: Jan Simon 2017, License: CC BY-SA 3.0
#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
// GetMaxIndex(r, dict)
// [~, ind] = max(abs(dict' * r));
double *r, *dict; // Inputs
size_t nr, ndict;
double *rEnd, *s, dot, maxValue; // Working
size_t i, maxIndex;
r = mxGetPr(prhs[0]); // Obtain inputs
nr = mxGetNumberOfElements(prhs[0]);
rEnd = r + nr;
dict = mxGetPr(prhs[1]);
ndict = mxGetN(prhs[1]);
maxValue = 0; // Calculation
maxIndex = 0;
for (i = 0; i < ndict; i++) {
dot = 0.0; // Dot product
for (s = r; s < rEnd; dot += *s++ * *dict++) ; // empty loop
dot = fabs(dot);
if (dot > maxValue) { // Remember maximum index
maxValue = dot;
maxIndex = i;
}
}
plhs[0] = mxCreateDoubleScalar((double) (maxIndex + 1));
return;
}
Please note the 2 comments "Not exactly the same!!!". I've commented the sections out, although they are faster and differ in the magnitude of eps() only in the resutls. But this changes the output final already, which means, that your algorithm is instable: The result depends on tiny rounding effects. Does this satisfy your needs?
For the first 15 iterations for ii this needs 10.8 sec instead of 31.18 sec - R2016b/64/Win7 without Parallel Computing Toolbox. The full problem will need 150 times longer divided by the number of cores.
