Can someone explain this large difference in compute time?

The two computations below should be identical, but tic/toc reports very different compute times for them in R2018a. Why?

function test

      N=3e4+1;
      F1=tril( triu(ones(N),0), 0);
      F2=diag(ones(N,1));

      x=rand(1e4,1); x(N)=0;

      tic; 
       z1=x.'*(F1*x)+7;
      toc;
      %Elapsed time is 0.231846 seconds.

      tic; 
       z2=x.'*(F2*x)+7;
      toc;
      %Elapsed time is 2.148849 seconds.

end

[OCDER]

That's very interesting. If you copy F2 as another variable, F3, you get the same speed boost as F1.

N = 10000;
x = rand(N,1); 
x(end) = 0;

F1 = tril(triu(ones(N),0), 0);
F2 = diag(ones(N,1),0);
F3 = F2;

tic 
    z1 = x.'*(F1*x)+7; 
toc %Elapsed time is 0.049871 seconds.

tic 
    z2 = x.'*(F2*x)+7; 
toc %Elapsed time is 0.138747 seconds.

tic 
    z3 = x.'*(F3*x)+7; 
toc %Elapsed time is 0.043385 seconds.

if isequal(F1, F2); disp('F1 = F2, Checked'); end
if isequal(F2, F3); disp('F2 = F3, Checked'); end
if isequal(z1, z2); disp('z1 = z2, Checked'); end
if isequal(z2, z3); disp('z2 = z3, Checked'); end

[Bruno Luong]

Prepare for this:

N = 10000;
x = rand(N,1); 
x(end) = 0;

F1 = tril(triu(ones(N),0), 0);
F2 = diag(ones(N,1),0);
F3 = F2;

tic 
    z3 = x.'*(F3*x)+7; 
toc %Elapsed time is 0.216113 seconds.

tic 
    z1 = x.'*(F1*x)+7; 
toc %Elapsed time is 0.021705 seconds.

tic 
    z2 = x.'*(F2*x)+7; 
toc %Elapsed time is 0.021853 seconds.

[Bruno Luong]

I read in a BLOG written by someone on TMW that discuss about of storing the results of a function call and returning the same results without computing again if the function is invoked with the same input arguments.

Can it explains that?

[Matt J]

You mean memoize()?

[Matt J]

More fun stuff:

    N=1e4+1;
    F1=tril( triu(ones(N),0), 0);
    F2=diag(ones(N,1),0);

    x=rand(1e4,1); x(N)=0;

    tic; 
     z1=x.'*(F1*x)+7;
    toc; %Elapsed time is 0.026374 seconds.

    tic; 
     z2=x.'*(F2*x)+7;
    toc; %Elapsed time is 0.247395 seconds.

    isequal(F1,F2);

    tic; 
     z2=x.'*(F2*x)+7;
    toc; %Elapsed time is 0.026936 seconds.

[Bruno Luong]

"You mean memoize()?"

That's it: Loren's blog.

[Matt J]

I don't see how that could be playing a role here.

[Bruno Luong]

Because the first call with F2/F3 takes more time regardless which one. They shares mxArray header.

[Matt J]

But in my original post, the first call takes the shortest time. Also, the compute time is unnaturally long for F2. Memoization should never lengthen compute time.

[Bruno Luong]

But I my comment is not applicable your OP, I comment on OCER and mine finding.

[OCDER]

Interesting finding and perhaps another hint. Just accessing F2 via a function like round, floor, ceil fixes the issue again.

N = 10000;
x = rand(N,1); 
x(end) = 0;
F1 = tril(triu(ones(N),0), 0);
F2 = diag(ones(N,1),0);
F3 = round(F2); %<=== doing this round here affects F2 too, somehow...
tic 
    z1 = x.'*(F1*x)+7; 
toc %Elapsed time is 0.048131 seconds.
tic 
    z2 = x.'*(F2*x)+7; 
toc %Elapsed time is 0.050953 seconds.  <=== time is faster
tic 
    z3 = x.'*(F3*x)+7; 
toc %Elapsed time is 0.045269 seconds.

Answer 1

[Matt J]

I think I know what's happening. Apparently, diag() doesn't return it's full result until this result is later used for actual computation. The output is stored in some sparse form until that happens. So, the tic/toc for the z1 calculation accounts for both the computation time and the memory allocation of the full matrix as well. The remarkable similarity between the total times for creation of F_i + calculation of z_i supports this:

      N=1e4+1;
      tic
      F1=diag(ones(N,1),0); 
      toc  %Elapsed time is 0.020561 seconds.

      tic
      F2=diag(ones(N,1),0)+0;
      toc  %Elapsed time is 0.266266 seconds.

      x=rand(1e4,1); x(N)=0;

      tic; 
       z1=x.'*(F1*x)+7;
      toc; %Elapsed time is 0.260155 seconds.

      tic; 
       z2=x.'*(F2*x)+7;
      toc; %Elapsed time is 0.028795 seconds.

[Matt J]

I've known for a while that this thing happens with zeros() (see below) but didn't think it extended to diag() as well.

    N=1e4+1;
    tic
    F1=zeros(N); 
    toc  %Elapsed time is 0.000103 seconds.

    tic
    F2=zeros(N)+0;
    toc  %Elapsed time is 0.257329 seconds.

    x=rand(1e4,1); x(N)=0;

    tic; 
     z1=x.'*(F1*x)+7;
    toc; %Elapsed time is 0.276728 seconds.

    tic; 
     z2=x.'*(F2*x)+7;
    toc; %Elapsed time is 0.033366 seconds.

[Bruno Luong]

+1, I think you get it.

So beside "copy-on-write" method, MATLAB does some sort of "create-on-use" on certain functions.

Answer 2

[Bruno Luong]

Interesting. Just a guess: TRIU and TRIL mark internally the resulting mxArray and MATLAB detects later the specific form and call specific faster BLAS routine when MTIMES is invokes.