Programming Efficiency

Documentation MathematicaBuilt-in Functions Advanced Documentation Linear Algebra Performance

This section discusses techniques for writing efficient Mathematica programs that solve linear algebra problems.

If you want to study the performance of your Mathematica code there are a number of timing related functions you can use.

Timing Mathematica operations.

The Timing function is particularly useful because it returns the CPU time needed for the computation.

In[1]:=

Out[2]=

Another useful function is TimeConstrained.

Time-constrained calculation.

This helps you to develop and test your code, shutting it down if it takes more than a preset limit.

In[3]:=

Out[4]=

Vectorizing Loops

A common operation involves some type of iteration over the elements of a matrix. Because the matrices may be large it is obviously important to do this efficiently.

One of the most important ways to write efficient code is to avoid explicit part references, especially in inner loops. This is a major source of inefficiency in user code and can be avoided by using functions that work on all the elements of a matrix or a vector at once. These vectorized operations allow Mathematica to use highly optimized code that runs faster.

List Creation

This problem involves creating a list and then initializing it with certain values. For example, computing

vi Sin[i]

One way to do this would be to create the list and then to use a For loop that iterates over its elements to set them to a value.

    n=10;
    v=Table[0,{n}];
    For[i=1,i≤n,i++,v[[i]]=Sin[2.0 Pi i]];

Slow way to initialize a list.

It is much faster to create the list and initialize it all at once; it is also neater.

n =10;
v=Table[ Sin[2.0 Pi i],{i,1,n}]

Faster way to initialize a list.

An even faster way is to use the optimized function Range to create the list and then to use vectorized operations.

n =10;
v=Sin[2.0Pi Range[1.,n]]

An even faster way to initialize a list.

Commands such as Table and Range are very optimized for building lists and matrices. They are discussed in the section Building Matrices.

List Updating

This problem involves updating the elements of a list; for example, computing the Sin of every element.

vi Sin[vi]

One way to do this involves iterating over the elements of the list replacing each with the updated version.

Do[v[[i]]=Sin[v[[i]]],{i,n}];

Slow way to update a list.

It is much faster to compute the new values by using a vectorized operation. The code is also neater and tidier.

v=Sin[v]

Faster way to update a list.

Another form of updating requires another argument that is not constant across the list.

One way to do this involves iterating over the elements of the list, dividing each by the square of the index.

n=Length[v];
Do[ v[[i]]=v[[i]]/i^2 ,{i,n}];

Slow way to update a list.

It is faster to compute the list of updates and then carry out a vectorized division.

d=Range[n]2;
v=v/d

Faster way to update a list.

The use of listable operations is discussed in the section Listability.

Using Built-in Support

Mathematica has a very wide range of functions, including many for list processing that have important applications for writing matrix computations. Before you start to work it is worth looking to see if any of the built-in functions will do what you want. For example, if what you are doing involves a convolution or correlation then you should use the appropriate functions; they will almost certainly be faster.

Another important point to remember is that Mathematica can represent a wide range of mathematical objects. For example, Mathematica can work with the actual linear equations that matrices are used to represent. Sometimes the equational form can be more expressive and convenient and this can lead to greater efficiency. The topic is discussed in the section Converting Equations to Sparse Arrays.

One area of built-in functions with many applications for linear algebra computations are those for structural manipulation of lists. Mathematica functions such as Part, Take, Drop, Transpose, PadLeft, PadRight, RotateLeft, and RotateRight are all heavily optimized and have many applications. These are discussed in more detail in the section Structural Operations. It is usually good if you can write your computations to use these functions.

A worked example will now be given that discusses different techniques to extend a matrix by tiling copies. For example, the input matrix

should be padded out to return the following matrix.

Three different ways to solve the problem will be demonstrated.

A Slow Way

The slowest way to do this is to write this in a step by step approach directly manipulating every element.

In[1]:=

In[2]:=

Out[3]//MatrixForm=

In[4]:=

Out[5]=

This way suffers from several speed deficiencies, the most severe being the inner loops over the part indices. These would be much faster if they were replaced by vectorized code. It was explained earlier that explicit part references in an inner loop can lead to inefficient code.

One thing to note about the function is that it will work when given a sparse matrix as input, but it will return a dense matrix as output.

In[6]:=

Out[7]=

A Faster Way

A much faster way to replicate the matrix is to write it to use Mathematica vectorizing functions.

In[1]:=

In[2]:=

Out[3]=

In[4]:=

Out[5]=

This method is much faster; it uses the function Part, which is very flexible and heavily optimized. It is described further in the section Getting Pieces of a Matrix.

This function works for a sparse matrix as input, and it returns the result as a sparse array which maintains equivalent sparsity.

In[6]:=

Out[7]=

Also Fast but Neater

Another way to extend the matrix, which involves no programming, is to use the built-in function PadRight.

In[1]:=

Out[2]=

In[3]:=

Out[4]=

This method is also neater and is probably the fastest for small input matrices. For large matrices, it is similar to the speed of the Part example shown above. The use of PadRight is described in the section Extending Matrices.

This function works for a sparse matrix as input, and it returns the result as a sparse array which maintains equivalent sparsity.

In[5]:=

Out[6]=

Matrix Contents

Matrices in Mathematica can be constructed from all the different types of object that Mathematica holds. They can contain machine-precision floating point numbers, arbitrary-precision floating point numbers, complex floating point numbers, integers, rational numbers, and general symbolic quantities. The different types are discussed in the section Matrix Types.

That Mathematica can work with these different types of matrix is a great strength of the system. However, it has a number of implications for the efficiency of numerical floating point computation: first from a mix of symbolic and numerical entries, secondly from a mix of different types of numbers, and thirdly from integer matrices.

Mixed Symbolic/Numerical Matrices

If you mix symbolic and numerical elements in a matrix, the linear algebra functions will treat the matrix as symbolic as discussed in the section Mixed Mode Matrices. It should be noted that techniques for working with symbolic matrices are typically slower than the heavily optimized techniques available for numerical matrices.

An example here compares matrix multiplication of two 200200 matrices. The second computation, which involves a symbolic matrix, is much slower.

In[1]:=

Out[2]=

In[3]:=

Out[4]=

The solution is to make sure that your matrices only contain numbers.

Mixed Numerical Type Matrices

If you mix different types of numbers in your matrices the linear algebra functions will treat the matrix as numerical and will coerce the matrix to the lowest precision. This is discussed in the section Mixed Mode Matrices. For linear algebra functions this will not have as severe an impact on performance as when working with symbolic matrices (this was demonstrated in the previous section).

This example compares matrix multiplication of two 200200 matrices. The second matrix, which contains an integer, is slower for matrix/matrix multiplication. If the operation was something more costly, this difference might not be significant.

In[1]:=

Out[2]=

In[3]:=

Out[4]=

The cost of mixed numerical matrices comes because Mathematica cannot use its efficient storage techniques, as discussed in the section Packed Arrays. This is demonstrated in the following example. Here, a vector of numbers is allocated and this is updated by some process. Because the updated values will be real numbers and the initial vector contained integers, the process cannot use packed arrays.

In[5]:=

Out[7]=

In[8]:=

Out[8]=

In this version the vector is initialized with real numbers. Now the computation can use packed arrays throughout and it runs faster.

In[9]:=

Out[11]=

In[12]:=

Out[12]=

The solution to this problem is to make sure you initialize with real numbers if you are going to work with real numbers.

Integer Matrices

A final type of problem involves working with purely integer matrices. Mathematica makes a distinction between working with integer and approximate real matrices; this is discussed in the section Matrices in Mathematica. Because an integer matrix will use symbolic computation techniques, it will be slower. Of course, if Mathematica did not use integer matrix techniques with matrices of integers it would not be so general.

In[1]:=

Out[2]=

In[3]:=

Out[5]=

The difference between the two computations can be compared. The eigenvalues of the integer matrix are returned with the exact values; in this case this involves radicals. The eigenvalues of the real matrix are returned with approximate real values.

In[6]:=

Out[6]=

In[7]:=

Out[7]=

Of course, it is a strength of Mathematica that it can return the exact result for a computation involving an integer matrix. But it should be noted that many purely numerical systems would return the approximate result.

The solution is that to work with matrices of approximate real numbers, you should start with approximate real numbers.