Implicit expansion/broadcasting of operands

Development
1

Hello Scilab users,

I would like to inform you that we have on gitlab a merge request proposed by @antoine.elias after some discussions we had inside the Scilab team: expand on math operators (!1606) · Merge requests · scilab / scilab · GitLab. This MR introduces a new behavior of arithmetical operators which is supported by Matlab since R2016b. To make a small resume, the new behavior allows singleton dimensions of operands to be changed so that the operand are considered as if they had the same size. You already know this behavior as it is already supported in Scilab, but only in the case when one of the operands is a scalar, as in the following

--> 1 + [1,2,3]

 ans = [1x3 double]

   2.   3.   4.

First, why is it worth generalizing this mechanism ? Let me consider a classical example, where you have to center and normalize columns of a matrix, you can do it by copying the row vector of means to form a matrix of same size as X

X = rand(3,5);
X = X - mean(X,1)([1 1 1],:);

or use a loop

X = rand(3,5);
m = mean(X,1);
for i = 1:3
    X(i,:) = X(i,:) - m;
end

but in both case the solution is inefficient (memory waste or explicit Scilab loop) when the size of arrays becomes large. If you further have to normalize columns, you face the same alternative (hell or … hell !)

X = X ./ stdev(X,1)([1 1 1],:);

or

s = stdev(X,1);
for i = 1:3
    X(i,:) = X(i,:) ./ s;
end

Dotted multiplication, division and power operators

I will now explain what proposes this MR in its actual state by considering first the dotted multiplication and right or left division operators .*, ./,.\,.^. These operators are given the feature of implicitely expand/broadcast operand when needed, without making actual copies in memory. This mechanism works only by increasing any singleton dimension so that operands get the same size and dimension. For example normalizing the columns of matrix X above is done by

X = rand(3,5);
X = X ./ stdev(X,1)

 X = [3x5 double]

   0.541873    1.7981974   8.1790804   0.2152067   2.709122 
   1.9386255   3.6221267   6.6003891   1.7652658   0.7404133
   0.000567    3.4207685   8.4531256   2.0847624   2.0299554

and yields the same result as

X = X ./ std(X,1)([1 1 1],:)

 X = [3x5 double]

   0.541873    1.7981974   8.1790804   0.2152067   2.709122 
   1.9386255   3.6221267   6.6003891   1.7652658   0.7404133
   0.000567    3.4207685   8.4531256   2.0847624   2.0299554

but without having to copy anything in memory. The new convention is that the dotted operator notation denotes not only element-wise operation but also implicit expand if dimensions of operands allow it.

Let us consider another simpler example:

[1 2 3] .* [4;5;6]

 ans = [3x3 double]

   4.   8.    12.
   5.   10.   15.
   6.   12.   18.

There is another way of constructing this rank one matrix by using the classical (not dotted) multiplication operator * but you have probably already understood that what could have been a Scilab beginner error will make sense with the new behavior of dotted multiplication. However, Scilab beginners are less prone to use dot prefixed operators, hence I don’t consider this case as a real problem.

Addition and substraction

Julia’s choice

It would have been consistent to introduce new dotted operators .+ and .- as it has been done in Julia language (https://julialang.org/ syntax is closely related to Matlab and Scilab ones but the language has an orthodox approach of everything and especially of such operators issues, which makes it very consistent). In Julia you can write

julia> [1 2 3] .+ [4;5;6]
3×3 Matrix{Int64}:
 5  6  7
 6  7  8
 7  8  9

and you are aware that something special will occur since you have used .+ instead of +.

Matlab’s R2016b choice

The choice that has been made by Matlab since R2016b, which is the same that has been made in the MR we are debating here, is to give the implicit expand feature to the classical, un-dotted addition and substraction operators + and -. Mathworks probably made this choice to stay consistent with the legacy of scalar expansion when addition/substraction is considered (I started this post with this point) and that’s where Julia approach differs, because in Julia even the scalar expand needs to use the dotted prefixed addition:

julia> 1 + [1 2 3]
ERROR: MethodError: no method matching +(::Int64, ::Matrix{Int64})
For element-wise addition, use broadcasting with dot syntax: scalar .+ array
The function `+` exists, but no method is defined for this combination of argument types.

julia> 1 .+ [1 2 3]
1×3 Matrix{Int64}:
 2  3  4

please remark the very explicit error message, which proposes to use dotted operator instead.

Now, the debate

Almost all users of Scilab have codes which rely on this implicit scalar expansion in addition, and that’s why any evolution of the behavior of + and - must keep this exception. Giving the implicit expand feature to + and - and not introducing .+ and .- (for other cases, when none of the operands is a scalar) is debatable. To make up you mind, I give you some links to related discussions:

https://www.mathworks.com/matlabcentral/answers/383533-how-do-i-avoid-getting-fooled-by-implicit-expansion
https://www.mathworks.com/matlabcentral/answers/318495-deactivation-switch-for-implicit-expansion

https://blogs.mathworks.com/loren/2016/10/24/matlab-arithmetic-expands-in-r2016b/
https://blogs.mathworks.com/loren/2016/11/10/more_thoughts_about_implicit_expansion/

If you want you can also download the Windows build of the MR
https://gitlab.com/scilab/scilab/-/jobs/12647188700/artifacts/browse
and experiment the proposed features, which are a great improvement for Scilab (thanks @antoine.elias). But the question is not about the features, but rather of their materialization at the Scilab user level (for all kind of users).

Sincerely yours,

S.

3 Likes
2

Hello,

I am very much in favor of this type of change. I have some ray-tracing code that could be significantly simplified—and, I imagine, made much faster—by this update. Typical examples include additions or multiplications by constant vectors:

function y=add_const(x, cst)
    y=zeros(x);
    y(1,:)=x(1,:)+cst(1);
    y(2,:)=x(2,:)+cst(2);
    y(3,:)=x(3,:)+cst(3);
endfunction

function y=mul_const(x, cst)
    y=zeros(x);
    y(1,:)=x(1,:)*cst(1);
    y(2,:)=x(2,:)*cst(2);
    y(3,:)=x(3,:)*cst(3);
endfunction
3

Hello Philippe,

On your use case, with a matrix with n=1e6 columns, the speedup is 200x compared to your solution based on loops, but only 4 compared to solution based on explicit expansion (on a MacBook pro M4):

n=1e6;
m=3;
x=rand(m,n);
y=zeros(m,n);
c=rand(m,1);

NR=5;

// with implicit expansion
tic
for i=1:NR
    y=x+c;
end
t(1)=toc();

// with loops
tic
for i=1:NR
    for j=1:n
        y(:,j)=x(:,j)+c;
    end
end
t(2)=toc();
printf("ratio (loops, %d lines) = %g\n",m,t(2)/t(1))

// with explicit expansion
j=ones(1,n);
tic;
for i=1:NR
    y=x+c(:,j);
end
toc
t(3)=toc()
printf("ratio (explicit expansion, %d lines)= %g\n",m,t(3)/t(1))

However, the first ratio decreases with the number of lines, although the second is quite stable:

ratio (loops, 3 lines) = 205.833
ratio (explicit expansion, 3 lines)= 4.66898

ratio (loops, 10 lines) = 89.4144
ratio (explicit expansion, 10 lines)= 6.34653

ratio (loops, 100 lines) = 8.64493
ratio (explicit expansion, 100 lines)= 4.10969

ratio (loops, 1000 lines) = 2.94557
ratio (explicit expansion, 1000 lines)= 5.76218