Matlab: Odd linear indexing into array with singleton dimensions

Question

I'm having trouble understanding the circumstances under which linear indexing applies.

For a 2D array, it seems intuitive:

>> clear x
>> x=[1 2;3 4]
x =
     1     2
     3     4
>> x([1 2])
ans =
     1     3
>> x([1;2])
ans =
     1
     3

For a 3D matrix, it seems intuitive:

>> clear y
>> y(1:2,1:2,1)=x
y =
     1     2
     3     4
>> y(1:2,1:2,2)=x+10
y(:,:,1) =
     1     2
     3     4
y(:,:,2) =
    11    12
    13    14
>> y([1 2])
ans =
     1     3
>> y([1;2])
ans =
     1
     3

For a 3D matrix in which the 1st two dimensions are singletons, it's not what I would expect:

>> clear y
>> y(1,1,1)=1
y =
     1
>> y(1,1,2)=2
y(:,:,1) =
     1
y(:,:,2) =
     2
>> y([1 2])
ans(:,:,1) =
     1
ans(:,:,2) =
     2
>> y([1;2])
ans(:,:,1) =
     1
ans(:,:,2) =
     2

I would have expected exactly the same as the 3D matrix without singleton dimensions.

Is there a rule that can be relied on to predict the behaviour of linear indexing?

Answer 1

The rule for linear indexing of an array x with an array ind is as follows (taken from this great post by Loren Shure):

• If at least one of x and ind is not a vector (that is, if x or ind have more than one non-singleton dimension) the output has the same shape (size) as ind:

  >> x = rand(3,4);
  >> ind = [9 3];
  >> x(ind)
  ans =
     0.276922984960890   0.743132468124916
  
  >> x(ind.')
  ans =
     0.276922984960890
     0.743132468124916
  
  >> x = rand(2,3,4);
  >> ind = [1 2; 5 6];
  >> x(ind)
  ans =
     0.814723686393179   0.905791937075619
     0.632359246225410   0.097540404999410
  

• If both x and ind are vectors, the output is a vector that has the same orientation as x (i.e. the output dimension which is non-singleton is that in x), and the same number of elements as ind:

  >> x = 10:10:70;
  >> ind = [1 3 5];
  >> x(ind)
  ans =
       10    30    50
  
  >> x(ind.')
  ans =
       10    30    50
  
  >> x = reshape(10:10:70,1,1,[]); % 1×1×7
  >> ind = reshape([2 3 4],1,1,1,1,[]); % 1×1×1×1×3
  >> x(ind)
  ans(:,:,1) =
      20
  ans(:,:,2) =
      30
  ans(:,:,3) =
      40
  

---

For completeness, if two or more indexing arrays are applied (so this is not linear indexing anymore), the shape of the output is determined by which dimensions of the original array are being indexed and by the number of elements of each indexing array. It is independent of the shape of the indexing arrays, which are simply read in column-major order as usual:

>> x = [10 20 30 40; 50 60 70 80];
>> x(2, [1 2; 3 4])
ans =
    50    70    60    80

>> x(2, [1 2 3 4])
ans =
    50    60    70    80

>> x(2, reshape([1 2 3 4],1,1,1,[]))
ans =
    50    60    70    80

If there are fewer indices than the number of dimensions of x, the trailing dimensions of x are implicitly collapsed into one before the indexing is applied:

>> x = [10 20 30 40; 50 60 70 80];
>> x(:,:,2) = x+100;
>> x([1 2], [1; 5; 7; 8])
ans =
    10   110   130   140
    50   150   170   180