Why this Octave code won't work?

Question

Let Y be a vector of length N, containing numbers from 1 to 10. As example code you can use:

Y = vec(1:10);

I am writing the code which must create an N x 10 matrix, each row consisting of all zeros except for a 1 only in the position which corresponds to the number in vector Y. Thus, 1 in Y becomes 10000000000, 3 becomes 0010000000, and so on.

This approach works:

cell2mat(arrayfun(@(x)eye(10)(x,:), Y, 'UniformOutput', false))

My next idea was to "optimize", so eye(10) is not generated N times, and I wrote this:

theEye = eye(10);
cell2mat(arrayfun(@(x)theEye(x,:), Y, 'UniformOutput', false))

However, now Octave is giving me error:

error: can't perform indexing operations for diagonal matrix type
error: evaluating argument list element number 1

Why do I get this error? What is wrong?

Bonus questions — do you see a better way to do what I am doing? Is my attempt to optimize making things easier for Octave?

Answer 1

I ran this code in Octave and eye creates a matrix of a class (or whatever this is) known as a Diagonal Matrix:

octave:3> theEye = eye(10);
octave:4> theEye
theEye =

Diagonal Matrix

   1   0   0   0   0   0   0   0   0   0
   0   1   0   0   0   0   0   0   0   0
   0   0   1   0   0   0   0   0   0   0
   0   0   0   1   0   0   0   0   0   0
   0   0   0   0   1   0   0   0   0   0
   0   0   0   0   0   1   0   0   0   0
   0   0   0   0   0   0   1   0   0   0
   0   0   0   0   0   0   0   1   0   0
   0   0   0   0   0   0   0   0   1   0
   0   0   0   0   0   0   0   0   0   1

In fact, the documentation for Octave says that if the matrix is diagonal, a special object is created to handle the diagonal matrices instead of a standard matrix: https://www.gnu.org/software/octave/doc/interpreter/Creating-Diagonal-Matrices.html

What's interesting is that we can slice into this matrix outside of the arrayfun call, regardless of it being in a separate class.

octave:1> theEye = eye(10);
octave:2> theEye(1,:)
ans =

Diagonal Matrix

   1   0   0   0   0   0   0   0   0   0

However, as soon as we put this into an arrayfun call, it decides to crap out:

octave:5> arrayfun(@(x)theEye(x,:), 1:3, 'uni', 0)
error: can't perform indexing operations for diagonal matrix type

This to me doesn't make any sense, especially since we can slice into it outside of arrayfun. One may suspect that it has something to do with arrayfun and since you are specifying UniformOutput to be false, a cell array of elements is returned per element in Y and perhaps something is going wrong when storing these slices into each cell array element.

However, this doesn't seem to be the culprit either. I took the first three rows of theEye, placed them into a cell array and merged them together using cell2mat:

octave:6> cell2mat({theEye(1,:); theEye(2,:); theEye(3,:)})
ans =

   1   0   0   0   0   0   0   0   0   0
   0   1   0   0   0   0   0   0   0   0
   0   0   1   0   0   0   0   0   0   0

As such, I suspect that it may be some sort of internal bug (if you could call it that...). Thanks to user carandraug (see comment above), this is indeed a bug and it has been reported: https://savannah.gnu.org/bugs/?47510. What may also provide insight is that this code runs as expected in MATLAB.

In any case, one thing you can take away from this is that I would seriously refrain from using cell2mat. Just use straight up indexing:

Y = vec(1:10);
theEye = eye(10);
out = theEye(Y,:);

This would index into theEye and extract out the relevant rows stored in Y and create a matrix where each row is zero except for the corresponding value seen in each element Y.

Also, have a look at this post for a similar example: Replace specific columns in a matrix with a constant column vector

However, it is defined over the columns instead of the rows, but it's very similar to what you want to achieve.

Answer 2

Another approach; We start with the data:

>> len = 10;                % max number
>> vec = randi(len, [1 7])  % vector of numbers
vec =
     1    10     9     5     7     3     6

Now we build the indicator matrix:

>> I = full(sparse(1:numel(vec), vec, 1, numel(vec), len))
I =
     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     1     0     0     0     0     0
     0     0     0     0     0     0     1     0     0     0
     0     0     1     0     0     0     0     0     0     0
     0     0     0     0     0     1     0     0     0     0