Last week I received email from a blog reader asking about extending the meaning of the arithmetic operators to do more than scalar expansion, and instead extend the expansion to singleton dimensions. In this post, I discuss the evolution of scalar expansion in MATLAB, talk about possible future designs, and open up the discussion to your input. In that spirit, please add your comments on this topic here.
Contents
Scalar Expansion for Arithmetic
In very early versions of MATLAB, you could rely on scalars expanding when you were performing element-wise arithmetic. In the following code, for example,
A = magic(3);
B = 2;
C = A + B;
the scalar B is added to each element of A.
Scalar Expansion for Assignment
It wasn't until MATLAB version 5 and the introduction of N-dimensional arrays, that scalar expansion extended to assignment, e.g.,
D(3:5,1:2) = 17
D =
0 0
0 0
17 17
17 17
17 17
replacing the need for forming a right-hand side exactly the same size as the left-hand side.
What about More than Scalar Expansion?
Judging from the number of newsgroup posts (see a few at these links: 1 2 3 4 5 ), there is a lot of interest in having way to express calculations concisely that don't involve scalars or arrays that match completely but for which there is a reasonable interpretation, for example, calculating A - mean(A) where A is a matrix.
How to Add Constant Values to Each Column of a Matrix
Let's suppose we want to add a constant value to each column of a matrix. Here are some ways we might do it today, all of which require at least some temporary extra memory from MATLAB. If the matrix is large, then the temporary of the second variable will also be large, and may consume a substantial additional amount of memory.
row = 1:4; mat = [magic(3) (10:12)']
mat =
8 1 6 10
3 5 7 11
4 9 2 12
Loop over columns of matrix, adding a scalar value to each one
out = mat; for c = 1:length(row) out(:,c) = out(:,c) + row(c); end
Expand the rows into a larger matrix the size of mat 3 different ways.
The first way is via indexing into the row vector, a technique known as Tony's trick because a customer, Tony, showed it to me at the ICASSP conference in Albuquerque in 1990
out1 = mat + row(ones(1,size(mat,1)),:);
The second way it to use repmat which essentially uses Tony's trick in most cases.
out2 = mat + repmat(row, size(mat,1),1);
The third variant is to use an outer product to duplicate the rows.
out3 = mat + ones(size(mat,1),1) * row;
Let's make sure our answers are all the same.
isequal(out,out1,out2,out3)
ans =
1
Doug Schwarz, a longtime MATLAB user, has developed a library of generalized operators which could also be used to solve this class of problem (and more).
Four Choices
There are four alternatives I can think of for dealing with expansion beyond scalars.
- Keep the status quo since there are already many ways to accomplish what users want.
- Create new functions to flexibly deal with combining different size arrays for elementwise (and perhaps matrix) operations.
- Change the basic operators in MATLAB such as + to accomodate commensurate but not identically sized arrays.
- Invent new operator notation to accomodate the sort of binary elementwise operations we've discussed so far.
Pros and Cons
Each of the above ideas has pros and cons. We've talked to some users in the past about this topic and their reactions generally fell into three buckets:
- loved the idea of changing the meaning of the elementwise operators so, for example, row+column would do an "outer" sum.
- hated the idea of changing the meaning of existing operators, for reasons of both readability and backward compatibility.
- indifferent (this was, by far, the smallest group!).
Now It's Your Turn
Now it's your turn to chime in. Please leave your thoughts on this topic here.
Regarding “hated the idea of changing the meaning of existing operators, for reasons of both readability and backward compatibility”:
If you stick to only defining new operators for cases which return an error now can not lead to any real backward compatibility issues, I think.
Taking your example, if you define “A-mean(A)” for all matrices in a way that is consistent with what happens now for vectors, I think most users will be quite happy.
Personally, I always liked the expressiveness of MATLAB, even when it is not too intuitive. I mean, creating a reandom permutaion matrix by “a=eye(5); a(randperm(5),:)” is not exactly something I would’ve thought up on myself, but it’s nice, still. I think that what is called “discoverability” is the problem here.
I would suggest defining new operators… users must explicitly say that they want this or that type of expansion. To much freedom is not good, the errors are very hard to find then.
MATLAB must be able to give warnings or error messages. If you redefine the “-” operator to encompass even the case of
A-mean(A), what if the user makes a mistake and writes A-mean(A(:))? MATLAB would have no objections. If instead, the user must write e.g. A -| mean(A) with -| meaning a “vector expansion” then the mistake in A -| mean(A(:)) would be easilly discovered.
I would love to see choices 3 and 4 implemented:
“3 Change the basic operators in MATLAB such as + to accomodate commensurate but not identically sized arrays.
4 Invent new operator notation to accomodate the sort of binary elementwise operations we’ve discussed so far”
Regarding 3 I would love to see also the “*” operator to accomodate contraction and outer products of general N-dimentional tensors. For example if we define a 3rd order tensor
A = randn(3,4,5);
% And a vector (I know, the last singleton dim is dropped)
x = randn(1,4,1);
% Then the result after contraction (summation) over the 2nd dimention should be a 3×5 tensor
B = A*x;
This notation is difficult to generalize to arbitrary tensors. One could imagine that the product A*B where
A = rand(2,3,1,1);
B = rand(1,1,4,5);
Would be the outer product meaning an outer product but this is not very clean.
Regarding choice 4 then I nice thing would be new notation such as *4 for product with contraction over dim 4. Essentially take the excellent Matlab Tensor toolbox at:
http://csmr.ca.sandia.gov/~tgkolda/TensorToolbox/
and implement it using custom operators such as *k etc.
Personally I’m indifferent to options 2-4 as long as something is done that is memory efficient. I do not think we should keep things status quo since the current solutions except for Doug’s require substantial amounts of memory when dealing with large matrices.
I would also like to see this expanded into indexing as well. If I have a 5 x 4 x 10 matrix and I have indices calculated for the 5 x 4 x 1 case, I would like the index to be smart about being reused. For example: x = rand([5 4 10]); indx = rand(5,4)>.5; x(indx) = NaN; Basically reuse the indx for each “page” in the third dimension.
Thanks for such a great topic.
Stephen
I agree with Petr Posik. The user should explicitly have to request this functionality (e.g. via new/extended operators), otherwise unintended side effects could make debugging very difficult in some cases.
Generalized iterators and nice method dispatch handling in matlab could solve part of the problem.
M = magic(3);
m = [1 2 3]
M.each_col(m, @(col, mi) col + mi );
M.each_row(m, @(row, mi) row + mi );
In R, you can do sloppy arithmetic with different length arrays. The shorter array is recycled to match the longer one. For example:
In Matlab the code
>> a=1:5;
>> b=6:8;
>> a+b
generates an error
??? Error using ==> plus
Matrix dimensions must agree.
but in R we get the following output.
[1] 7 9 11 10 12
Warning message:
longer object length
is not a multiple of shorter object length in: a + b
The sloppy arithmetic allows you to avoid using repmat, or the other techniques mentioned, and a warning prevents debugging nightmares.
The shorter-object-recycling even extends to matrices.
Much as I love Matlab, I think this is one thing R does better at present. Bring on the updated operators!
Since you asked…
At the very least, conformable arrays should multiply by their inner dimension.
>> r=randn(2,3,4);
>> s=randn(4,3,2);
>> r*s
??? Error using ==> mtimes
Input arguments must be 2-D.
I also like A - mean(A) for when such things make sense
Also… xcorr (and any “companions” it may have) should allow correlating the columns of 2 (or N) dimensional array with a 1-dimensional array
>> xcorr(randn(7),randn(7,1))
??? Error using ==> xcorr
When B is a vector, A must be a vector.
nice blog, loren
- an additional (set of) operator(s) may obfuscate the syntax
- a self-recognizing (fuzzy/sloppy) syntax may obfuscate the issue
- i’m in favor of a solution along the lines of the newly introduced family of construct-handlers such as ARRAYFUN, STRUCTFUN, CELLFUN: (eg) GOPFUN with similar syntax, great flexibility (user definable functions/options), easy internal error checking, mex-speed
us
Hi Loren,
Wonderful topic and one that is important to my work.
I would like to second Ur’s suggestion. I agree with all that he wrote. The flexibility of passing function handles and anonymous functions is also a great idea. I am often forced to write loops in matlab since I can use one of those nice solutions that reformat the matrices to be the same size due to memory limitations.
Stephen
Urs said it exactly right with his typical insight and elegance. I endorse his ideas. They are what makes sense.
One thing I’ve seen elsewhere that I miss in Matlab is named wildcards for indices, rather than just :.
Then most weird array contractions (and a couple of other things) can be implemented, possibly with the help the sum operator.
For example,
Cartesian product:
A(:_a, :_b) = B(:_a)*C(:_b)
Vector expansion:
A(:_a, :_b) = B(:_a,:_b)*C(:_b)
Multiplying two tensors across a given dimension:
A(:_a,:_b,:_d)=sum(B(:_a,:_c!,:_d), C(:_c!,:_b))
here :_a etc are the different index wildcards that show what index goes where, and the exclamation sign shows that you have to sum over that index. (In true physicist fashion, once ! is seen in such an expression, the summing might even be automatic, without needing to write sum())
This is sufficiently explicit that matlab does not have to guess, yet sufficiently flexible to handle all the tasks discussed.
Couldn’t disagree more strongly with Richie, who thought sloppy arithmetic is a good thing
In R, you can do sloppy
arithmetic with different length arrays.
The shorter array is recycled to match the
longer one.
You aren’t saved by “R”s warnings here. Since it will only issue them if the lengths are factors of the total size.
I agree with Urs
Back when this issue was being discussed on the newsgroup I liked the idea of defining new operators to specify functionality like that in my genops ones, perhaps |+ or !+ for generalized addition.
However, Urs’ idea has merit also and I would be interested in seeing that fleshed out a bit more. One advantage of it is that we could implement it today for all versions of MATLAB as no changes to the parser would be required.
In spite of the fact that I wrote my genops code so that the current operators are overloaded I am mildly opposed to having MATLAB work that way “out of the box”. I think this new syntax is different enough to warrant being used only when explicitly requested.
Thanks for bringing this up again, Loren!
Hi,
I support Egor Kraev’s idea:
- intuitive syntax
- explicit syntax (”sloppy” arithmetic may hide bugs)
- very flexible
- Thanks to ! standing for “summed index”, can be made memory efficient
- and I am confident that TMW can make if FAST.
- challenge: must work (efficiently) with function calls:
a(:_i,:_j) = NonLinFun(b(:_i),c(:_j)) where NonLinFun has been programmed to operate on SCALAR b and c.
A very important issue, that needs a carefully designed solution, soon.
Regards, Philippe
I want to thank everyone for your thoughtful contributions on this topic. Several of us constantly monitor all the comments and have found some excellent ideas to consider as we plan for MATLAB’s future.
Even though I am posting my thanks now, please don’t let that stop you from adding to this list if you want.
I attempted to implement what Urs suggested. See the file dotdot.m (# 10295) on the File Exchange. It’s an M-file so it doesn’t avoid the memory allocation that could be eliminated if it were a built in function. In any case it automates the generalized expansion and uses function handles. I especially like it when you have it return a function handle to do what you want, then you can call it with new data over and over and it takes care of it no matter what size the arguments are.
Duane
This is perhap the single function that could be most valuable as an addition to Matlab. I’ve argued for its presence before. I really like the ability to use a general operator in Duane’s implementation. As distributed, compiled code, it would remove the need for repmat.
Put it in Matlab. Please.
John
I just submitted by version of Urs’ idea to the file exchange. It does essentially the same thing as Duane’s (but I started working on it before yours was published, Duane — really!), but it manages to avoid copying any data. To do that I had to use eval (gasp!!) because my initial approach using subsref and subsasgn was much slower. Also, I stole Duane’s idea of returning a function handle for a single input. Thanks, Duane!
I still like the idea of having this functionality built into MATLAB because I think the appearance of
y = x |- mean(x);
is better than
y = genop(@minus,x,mean(x));
though both are better than
y = x - repmat(mean(x),size(x,1),1);
Doug
I did a quick comparison of Doug’s matlab only code vs c code versions of the genops for @plus and I’m seeing the matlab code version take 35 times longer to run then the c code version.
I appreciate having both around, but I’d like to have function handles be passed to the c version. Unfortunately I know making genop have the same flexibility of cellfun, arrayfun, etc… is well beyond what I can handle in mex functions. If the source code for the ones already in matlab were available maybe we can make it a community project to create a gopfun mex function.
Stephen
I’ve recently switched from Mathematica to Matlab as my primary programming tool, primarily because Matlab was orders of magnitude faster for my purposes (self organizing neural networks). However, I do think that Matlab developers ought to consider the way Mathematica “vectorizes”. Almost every function when given a N-D matrix performs its operations row-wise and via the various “functional programming” functions (Map, MapAt, MapThread, Apply, etc) which allow the coder to control exactly how, and at what level, functions operate on data structures. I’ve found it takes many more lines of code to vectorize in Matlab. Perhaps it’s due to my relative lack of familiarity with Matlab, but I do think Matlab might benefit by incorporating some of the vectorization methods from Mathematica.
How to calculate a cummulative sum of a matrix over over a given varible. Example, how I calculate year totals from monthly total. The data set look like
Country month Temp
1 1 10
. . .
1 12 12
2 1 60
. . .
2 12 72
Peter-
I don’t know what you are trying to calculate. Are the columns year, month, and number. You’d probably be best off reshaping your array into a 3-dimensional array and summing along the correct dimension(s) then.
–Loren
Regarding “How to Add [or multiply] Constant Values to Each Column of a Matrix.” This bothered me for a long time. I used to use what Loren calls “Tony’s Method” (both versions) for a while but got frustrated with it because of the performance hit for large matrices. Though I found the repmat method to be faster than looping, I couldn’t get over the fact that I was doing so many extra additions (or multiplications). Finally, I wrote a few mex functions (e.g., rowmult(cVec,Matrix), colmult(rVec,Matrix), etc.) that calculate these directly. These are really basic operations that I use so often that I’m really shocked that they aren’t built in to Matlab. While it might be nice to eventually come up with some operators to represent these functions, I think it’s more important that efficient functions implementing these operations be included even if only in functional form.
Loren,
Ah ha! I can’t believe I’ve been using MATLAB in anger, daily, for years now - and hadn’t figured out the outer product.
God bless you, this just saved me 586,000 seconds… And then the rest!
Tom
Tom-
See the function bsxfun for other scalar expansion functionality.
–Loren
Hi, Loren,
How to apply a function to each row of a matrix with out using the for loop?
Thanks,
Everett
Everett-
Many MATLAB functions already work on columns of an array or can work on a selected dimension (see help for max, for example).
–Loren
You can also check out my _apply_ function at the FEX.
http://www.mathworks.com/matlabcentral/fileexchange/16467
While it doesn’t speed anything up, it makes the code somewhat cleaner. I wish the mathworks would add something like this to base MATLAB.
Stephen