Does A(m:n) trigger the explicit creation of an index array (m:n)?

I was under the impression that in an indexing expression A(m:n) when A is a numeric vector, Matlab was smart enough to optimize the operation by not explicitly constructing an intermediate index array I=(m:n), as it must in more general indexing expressions. I'm sure this was discussed in Answers or newsgroup posts from long ago (circa 2009).

However, the only way I can make sense of the different timings below is that there is substantial overhead coming from the expressions (1:end). So, am I mistaken? Do all indexing expressions A(expr) result in the explicit creation and storage of expr, even if it consists of a colon expression?

runTest()

ans = 0.0490

ans = 0.1041

ans = 0.0111

function runTest

N=1e7;
[a]=deal( rand(N,1) );

timeit(@() f1(a))

timeit(@() f2(a))

timeit(@() f3(a))
    
    function f1(a)
    
     c=a(1:N)./2;
    
    end
    
    function f2(a)
    
     a(1:N)=a(1:N)./2;
    
    end
    
    function f3(a)
    
      a(:)=a(:)/2;
    
    end

end

[Stephen23]

https://www.mathworks.com/matlabcentral/answers/159410-oop-subsref-runs-twice-with-end

[Matt J]

@Stephen23 I don't think that link is applicable to built-in classes. In any case, I've modified the example by eliminating 'end', and the results are the same.

Answer 1

[James Lebak]

For numeric indexing, MATLAB does make a special optimization for a range expression like m:n. We do additional optimizations with a literal colon, and I think this is why your results don't match your expectations.

I timed the following functions.

function f4(a, N)
c=a(1:N)./2;
end

function f5(a, N)
idx = 1:N;
c = a(idx)./2;
end

On my machine with N = 1e7 I see

f4: 0.0283

f5: 0.1015

Comparing f4 and f5, the overhead that you see is the effect of creating an explicit index vector and using it as the index. We preserve the compressed nature of the indexing vector in f4 and that's why it's faster.

An additional optimization applies to your f3 function because of the use of the literal colon. When the compiler can prove that an indexing expression will apply to every element of the array, it is able to make further optimizations. This ends up being much faster than the more generic implementation that we use for general indexing expressions, and it's why your f3 is so much faster than your f1 or f2. We can't apply this optimization to a range expression or to an index vector, because the compiler can't prove that it will select every element of the array.

[Walter Roberson]

It would not have surprised me if MATLAB was able to optimize the case of constant numeric limits in a colon expression, but not the case of variable limits.

a(1:100)     %hypothesized to be optimized

N = 100; a(1:N)    %hypothesized to not be as optimized

[Matt J]

We do additional optimizations with a literal colon, and I think this is why your results don't match your expectations.

Thanks @James Lebak, but if there are optimizations for literal colons then that makes the differences between f1 and f2 confusing. Why is f2 twice as slow as f1? The only difference between them is that there is an extra literal colon the left hand side. If such expressions are optimized, I would not expect that to be a bottleneck.

[James Lebak]

This is very much an implementation detail, but basically, the output of the right-hand side of an assignment is always a temporary. That means that when you have a simple assignment to a variable like

c = ...

it's extremely cheap, just a pointer copy.

When you add an index to the variable, either a literal-colon

a(:) = ...

or a range expression

a(1:N) = ...

you are making an extra indexed copy of the right-hand side and performing elementwise assignment into the left-hand side. That is why f2 is twice as slow as f1.

[Matt J]

the output of the right-hand side of an assignment is always a temporary

Not always, right? The documentation says some assignments are done in place. That doesn't apply to literal colons?

[James Lebak]

Sorry, that was an overstatement on my part.

The documentation you pointed to is correct about in-place assignment and re-using a variable that appears both as an input and an output to a function.

A better way to say it is that the output of indexed reference is usually a temporary value. Occasionally (for example, when every index is a literal colon) it may be a shared copy of the indexee.

Adding indexing to the left-hand side of an assignment usually produces an extra copy operation.

[Matt J]

OK, thanks. I wonder what the potential might be to generalize in-place assignment optimization to indexed expressions. Could the existing engine be modified easily at least for cases where the variable is referenced with the exact same indexing expression on both the right and left hand sides?

[James Lebak]

I found an open enhancement request for this. I will add a mention of this answers post to the request as well. I don't think it's an easy modification.

[Steven Lord]

Could the existing engine be modified easily at least for cases where the variable is referenced with the exact same indexing expression on both the right and left hand sides?

I'm not sure that would always work.

A = ones(3);
whos A

  Name      Size            Bytes  Class     Attributes

  A         3x3                72  double              

A(5) = A(5)+1i;
whos A

  Name      Size            Bytes  Class     Attributes

  A         3x3               144  double    complex   

The second whos call shows that A requires twice as much memory. [You could also see an example where it uses less, if the assignment overwrote the sole complex value in the array with a real value by subtracting the imaginary part.]

A(5) = A(5)-1i;
whos A

  Name      Size            Bytes  Class     Attributes

  A         3x3                72  double              

So it's not always as easy as just writing over the value in memory.

[Matt J]

@Steven Lord That's a fair point, but the same issue is present for non-indexed susasgn expressions. According to the documentation, a statement like,

A=A+B;

currently results in an in-place assignment. Though not mentioned in the doc, I imagine based on your post that this does not apply when A and B are not both purely real.

But if the engine is already set up to parse and deal separately with the complex-valued case, why would it become harder just because an index is appended to A?

[Matt J]

Though not mentioned in the doc, I imagine based on your post that this does not apply when A and B are not both purely real.

Based on the test below, it appears that in-place optimizations are still applied even when the in-place variable oscillates between real and complex-valued content. The speed improvement factor is noticeably reduced, though.

    N=1e8;
    a=deal( rand(N,1) );

    f1(a,2); %always real 

Elapsed time is 0.917327 seconds.

    f2(a,2); %always real (IN-PLACE)

Elapsed time is 0.194343 seconds.

    f1(a,1i);  %oscillating complex->real->complex

Elapsed time is 3.098619 seconds.

    f2(a,1i);  %oscillating complex->real->complex (IN-PLACE)

Elapsed time is 2.396123 seconds.

function a=f1(a,s)

    tic;
    for i=1:10
        tmp=a./s;
        a=tmp;
    end
    toc

end

function a=f2(a,s)

    tic;
    for i=1:10
        a=a./s;
    end
    toc

end