# Optimize inner loops

I'm Mathematica newbie so please be gentle :) I have this, heavily non-optimized part of code, which I would like to speed up. I have put all matrices to be RandomReal, but in my code they take specific values. Also, matrices mat3, mat4, mat6 and mat7 consist of random simple trigonometric functions. Any kind of help is much appreciated. I have read about Map, Table, Do, Functional Programming... but I don't know how to apply it.

Code is part of Finite element-Finite strip calculation. I don't mind losing computing time during "real" calculation, but this part of code is just placing elements of matrices into other places with just few calculations. Matrices mat2, mat5, mat9, mat10 and mat11 need to be stored for later data processing.

Bonus question: is it possible to Compile this part of code?

Hope I gave optimal amount of data about my problem

All the best, Aleksandar

limit1 = 10;
limit2 = 20;
limit3 = 10;
limit4 = 15;
mat1 = RandomReal[{-100, 100}, {limit3, 2}];
mat2 = RandomReal[{-100, 100}, {limit1, limit2, limit3, 2}];
mat3 = Table[Sin[m*\[Pi]*y] + Cos[s*\[Pi]*y], {m, limit2}, {s, limit3}];
mat4 = Table[Sin[m*s*\[Pi]*y] + Cos[s*\[Pi]*y], {m, limit2}, {s, limit3}];
mat5 = RandomReal[{-100, 100}, {limit1, limit2, limit3, 6}];
mat6 = Table[Sin[m*m*\[Pi]*y] + Cos[s*\[Pi]*y], {m, limit2}, {s, limit3}];
mat7 = Table[Sin[m*\[Pi]*y] + Cos[s*s*\[Pi]*y], {m, limit2}, {s, limit3}];
mat8 = RandomReal[{-100, 100}, {limit2, limit3, limit4, 2}];
mat9 = RandomReal[{-100, 100}, {limit1, limit2, limit3, limit4}];
mat10 = RandomReal[{-100, 100}, {limit1, limit2, limit3, limit4}];
mat11 = RandomReal[{-100, 100}, {limit1, limit2, limit3, limit4}];

For[n = 1, n < limit1 + 1, n++,
For[i = 1, i < limit2 + 1, i++,

For[j = 1, j < limit3 + 1, j++,

y = (mat1[[j, 1]] + mat1[[j, 2]])/2;

mat2[[n, i, j, 1]] = mat3[[i, j]];

mat2[[n, i, j, 2]] = mat4[[i, j]];

mat5[[n, i, j, All]] = 1/2 (mat6[[i, j]] + mat7[[i, j]]);

Clear[y];

For[k = 1, k < limit4 + 1, k++,

zz = 1/2 (mat8[[i, j, k, 1]] + mat8[[i, j, k, 2]]);

mat9[[n, i, j, k]] = {mat5[[n, i, j, 1]] + zz*mat5[[n, i, j, 4]],
mat5[[n, i, j, 2]] + zz*mat5[[n, i, j, 5]],
mat5[[n, i, j, 3]] + zz*mat5[[n, i, j, 6]]};
mat10[[n, i, j, k]] = mat11[[n, i, j, k]].mat9[[n, i, j, k]];

Clear[zz]
]

]]
]

-
@NasserM.Abbasi He sets mat3,mat4,mat6 and mat7 using unevaluated y. Then when he sets y and sets mat 2 using mat3 and mat4, y evaluates using his definition. So y is being used, but it's not obvious from the code. – jVincent Nov 20 '12 at 11:45
i suspect there is a typo in the last line mat10[[n, i, j, k]] = mat11[[n, i, j, k]].mat9[[n, i, j, k]];, did you mean mat10[[n, i, j, k]] = mat11[[n, i, j, k]]*mat9[[n, i, j, k]]; ? – Thies Heidecke Nov 20 '12 at 21:43
not actually, mat11[[n, i, j, k]] and mat9[[n, i, j, k]] are matrix and vector, respectively :) – sasaborg Nov 20 '12 at 23:19

The code as it is now looks very much FORTRAN style, which is fine. But Mathematica offers you a wide range of ways to make your code more readable, faster and easier to spot potential bugs. So let's go through through some of the possible ways to improve your code:

Variable Naming

I know that in languages like C and FORTRAN it's common to give variables alphabetically and rather short names, but in my experience it helps readability of your code for your future you (when you maybe haven't looked at your code for a few weeks) and other persons that try to understand your code to try naming variables

• verbose
• selfexplaining

. So e.g. when i read limit1, limit2, etc. i have no clue what those numbers represent apart from being a kind of boundary probably. Later in your code i see that you use them as boundaries for your loop variables, so why not name them maxn, maxi, etc. for example? now i can anticipate that they are connected to loop variables, even before i saw the loops. The same goes for mat1, mat2, etc... It makes reasoning about code hard, because it's just abstract, and most people are better at problem solving when they can visualize stuff.

Making functions explicit

In your code in mat(3|4|6|7) you introduce the variable y

mat3 = Table[Sin[i*\[Pi]*y] + Cos[j*\[Pi]*y], {i, maxi}, {j, maxj}];


and use it later inexplicitly by changing the value of y and then using the changed value of mat3, e.g.

y = (mat1[[j, 1]] + mat1[[j, 2]])/2;
mat2[[n, i, j, 1]] = mat3[[i, j]];


While this is possible it makes it kind of hard to see the dependence of mat2 on y, it's even obscuring it almost.

So why not make the dependence explicit and save ourselves a lot of possible headaches later ;P In your example we could express the dependence like this:

mat3[y_] = Table[Sin[i*\[Pi]*y] + Cos[j*\[Pi]*y], {i, maxi}, {j, maxj}];
mat2[[n, i, j, 1]] = mat3[ (mat1[[j, 1]] + mat1[[j, 2]])/2 ][[i, j]];


which makes it obvious and explicit that mat3 is dependent on y, and likewise mat2 is dependent on mat1.

It's still kind of horrible, because we compute the whole matrix mat3 for every element of mat2 and throw away the rest. There are at least two ways out. One is making the dependence on i and j explicit, too:

mat3elem[i_,j_,y_] = Sin[i*\[Pi]*y] + Cos[j*\[Pi]*y];
mat2[[n, i, j, 1]] = mat3elem[i, j, (mat1[[j, 1]] + mat1[[j, 2]])/2]];


another way would be to get rid of the i- and j-loops and use vectorised operations that work on the whole matrix like halirutan stated in his answer.

Getting rid of loops by using Functional Programming constructs

In procedural languages like C and FORTRAN loops are used for various purposes, e.g.

traversing arrays and doing the same thing with every element, e.g. assigning values to an array based on the index or another array

n = 10;
array = ConstantArray[0, n];
For[i = 1, i <= n, i++, array[[i]] = i^2]

(* {1,4,9,16,25,36,49,64,81,100} *)


filtering arrays, e.g.

inputarray = {1, 3, 4, 5, 10, 4, 7, 9, 2, 13};
n = Length[inputarray];
outputarray = {};
For[i = 1, i < n, i++
If[ PrimeQ[inputarray[[i]]],
AppendTo[outputarray, inputarray[[i]]]
]
];
outputarray

(* {3,5,7,2,13} *)


repetition, e.g. iterating a function over some varaible to reach some fixpoint or get a list of intermediate values

x = 100.0;
For[i = 1; rootapprox = x/2, i < 7, i++,
rootapprox = (rootapprox + x/rootapprox)/2;
Print[rootapprox]
]


This variety of stuff makes loops on one hand very versatile and powerful, on the other side very unexplanatory constructs, that take extra time to see what they are used for.

Let's look at alternatives for each of the above cases:

traversing arrays

halirutan already showed vectorisation through Listable functions and Part as a way to apply operations over whole columns, rows or matrices, which is a great starter i would go for in most circumstances because it's concise and selfexplanatory. Another way to apply a function is by mapping it over a list of values, which for the above example could look like this:

n = 10;
array = Map[ Function[x, x^2], Range[n] ];

(* {1,4,9,16,25,36,49,64,81,100} *)


or shorthand

n = 10;
array = #^2 &  /@  Range[10]

(* {1,4,9,16,25,36,49,64,81,100} *)


At first it may look kind of arbitrary which style to use, but the functional style has one major advantage: it eliminates the need for intermediate state in form of variables that we only use for temporary results. In our case the array variable. We just say what our input is ( Range[10]=={1,2,3,...,9,10}) and what we want to do to each element (square it Function[x,x^2]==(#^2 &)) and get the result by mapping the function over our input values (via Map[] or /@').

filtering arrays

a way of doing filtering of arrays in a functional style is

inputarray = {1, 3, 4, 5, 10, 4, 7, 9, 2, 13};
Select[inputarray, PrimeQ]

(* {3,5,7,2,13} *)


by which we get rid of the need of explicitly indexing every element in the array and thus even the output variable outputarray altogether.

repetition

repetition can be done in a functional style via recursion, of which one often needed simple case can be achieved via Nest, which just describes the case, where we feed the output of a computation to itself as the new input. When we are just interested in the end result we can use

x = 100.0;
Nest[(# + x/#)/2 &, x/2, 6]

(* 10. *)


or if we are interested in the intermediate results, too

x = 100.0;
NestList[(# + x/#)/2 &, x/2, 6]

(* {50.,26.,14.9231,10.8121,10.0305,10.,10.} *)


.

Putting it all together

i'll add examples how you could apply these ideas to your code later. For now maybe you got some ideas to try out yourself.

-
Great answer..! – cormullion Nov 20 '12 at 17:07
Thanks, now I'm in lot of thinking :) This is what I needed, and comparison with FORTRAN is great (I was there). I never had a chance to learn functional programming so I will try now :) I will let you know about results. – sasaborg Nov 20 '12 at 23:25

This is awful. It is one very typical example of "how to use Mathematica the wrong way*. OK, enough complaining. Let me give you one hint. Lets say you have a 500x500 and a 1000x1000 matrix and you want to copy the smaller one in the upper left corner of the larger one. We do this step 100 times.

In your style this would go like

m1 = RandomReal[{0, 1}, {1000, 1000}];
m2 = RandomReal[{0, 1}, {500, 500}];

AbsoluteTiming[
Do[
For[j = 1, j <= 500, j++,
For[i = 1, i <= 500, i++,
m1[[j, i]] = m2[[j, i]]
]
], {100}]
]
(* {65.111817, Null} *)


This takes over 1 Minute here. What goes faster (in fact 2500x faster) is to do it in one chunk

AbsoluteTiming[
Do[
m1[[1 ;; 500, 1 ;; 500]] = m2,
{100}]
]
(* {0.025199, Null} *)


You can generalize this to any dimension. Here I use vectors and a matrix:

v1 = {1, 1, 1, 1};
v2 = {2, 2, 2, 2};
m = RandomInteger[{0, 10}, {4, 4}];


Setting the 2nd row

m[[2]] = v1 - 1/2 v2;
m // MatrixForm


Setting the 3rd column

m[[All, 3]] = v1 + v2;


Setting a submatrix

m[[1 ;; 2, 3 ;; 4]] = 9;
m // MatrixForm
`

and so on and so forth..

-
I was aware of the simple constructs as "Do" that you presented. But I don't know how to apply this to may problem (because I have some calculations inside the loop). Anyway, thanks for answer. – sasaborg Nov 20 '12 at 12:57
@sasaborg My examples are directly usable in the snip of code you gave. You can replace most loops with simple vector-assignments and vector-operations. – halirutan Nov 20 '12 at 13:10
Ok, I'll try to implement it and let you know how did it go. – sasaborg Nov 20 '12 at 13:33