Performance of vectorizing code to create a sparse matrix with a single 1 per row from a vector of indexes

Question

I have a large column vector y containing integer values from 1 to 10. I wanted to convert it to a matrix where each row is full of 0s except for a 1 at the index given by the value at the respective row of y.

This example should make it clearer:

y = [3; 4; 1; 10; 9; 9; 4; 2; ...]

% gets converted to:

Y = [
    0 0 1 0 0 0 0 0 0 0;
    0 0 0 1 0 0 0 0 0 0;
    1 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 1;
    0 0 0 0 0 0 0 0 1 0;
    0 0 0 0 0 0 0 0 1 0;
    0 0 0 1 0 0 0 0 0 0;
    0 1 0 0 0 0 0 0 0 0;
    ...
    ]

I have written the following code for this (it works):

m = length(y);
Y = zeros(m, 10);
for i = 1:m
    Y(i, y(i)) = 1;
end

I know there are ways I could remove the for loop in this code (vectorizing). This post contains a few, including something like:

Y = full(sparse(1:length(y), y, ones(length(y),1)));

But I had to convert y to doubles to be able to use this, and the result is actually about 3x slower than my "for" approach, using 10.000.000 as the length of y.

1. Is it likely that doing this kind of vectorization will lead to better performance for a very large y? I've read many times that vectorizing calculations leads to better performance (not only in MATLAB), but this kind of solution seems to result in more calculations.

2. Is there a way to actually improve performance over the for approach in this example? Maybe the problem here is simply that acting on doubles instead of ints isn't the best thing for comparison, but I couldn't find a way to use sparse otherwise.

Answer 1

Here is a test to comapre:

function [t,v] = testIndicatorMatrix()
    y = randi([1 10], [1e6 1], 'double');
    funcs = {
        @() func1(y);
        @() func2(y);
        @() func3(y);
        @() func4(y);
    };

    t = cellfun(@timeit, funcs, 'Uniform',true);
    v = cellfun(@feval, funcs, 'Uniform',false);
    assert(isequal(v{:}))
end

function Y = func1(y)
    m = numel(y);
    Y = zeros(m, 10);
    for i = 1:m
        Y(i, y(i)) = 1;
    end
end

function Y = func2(y)
    m = numel(y);
    Y = full(sparse(1:m, y, 1, m, 10, m));
end

function Y = func3(y)
    m = numel(y);
    Y = zeros(m,10);
    Y(sub2ind([m,10], (1:m).', y)) = 1;
end

function Y = func4(y)
    m = numel(y);
    Y = zeros(m,10);
    Y((y-1).*m + (1:m).') = 1;
end

I get:

>> testIndicatorMatrix
ans =
    0.0388
    0.1712
    0.0490
    0.0430

Such a simple for-loop can be dynamically JIT-compiled at runtime, and would run really fast (even slightly faster than vectorized code)!

Answer 2

It seems you are looking for that full numeric matrix Y as the output. So, you can try this approach -

m = numel(y);
Y1(m,10) = 0; %// Faster way to pre-allocate zeros than using function call `zeros`
  %// Source - http://undocumentedmatlab.com/blog/preallocation-performance
linear_idx = (y-1)*m+(1:m)'; %//'# since y is mentioned as a column vector, 
                              %// so directly y can be used instead of y(:)
Y1(linear_idx)=1; %// Y1 would be the desired output

---

Benchmarking

Using Amro's benchmark post and increasing the datasize a bit -

y = randi([1 10], [1.5e6 1], 'double');

And finally doing the faster pre-allocation scheme mentioned earlier of using Y(m,10)=0; instead of Y = zeros(m,10);, I got these results on my system -

>> testIndicatorMatrix
ans =
    0.1798
    0.4651
    0.1693
    0.1457

That is the vectorized approach mentioned here (the last one in the benchmark suite) is giving you more than 15% performance improvement over your for-loop code (the first one in the benchmark suite). So, if you are using large datasizes and intend to get full versions of sparse matrices, this approach would make sense (in my personal opinion).

Answer 3

Does something like this not work for you?

tic;
N = 1e6;
y = randperm( N );
Y = spalloc( N, N, N );
inds = sub2ind( size(Y), y(:), (1:N)' );
Y = sparse( 1:N, y, 1, N, N, N );
toc

The above outputs

Elapsed time is 0.144683 seconds.