MATLAB is an interpreted scripting language designed specifically for vector and matrix calculations common in science and engineering. MATLAB handles vectors and matrices with an intuitive syntax, offers a wide variety of built-in, pre-compiled routines for numerical computations, and has other useful features such as sophisticated graphing and plotting routines.

If used properly, and for its intended purpose (vector and matrix calculations), MATLAB is a very useful tool. Unfortunately, not all calculations can be done efficiently in MATLAB. Some operations are difficult or impossible to vectorize, and MATLAB has a finite number of built-in functions. Even if the code can be vectorized, doing so is often more difficult than writing the code in a compiled language.

Many people use MATLAB where it is ill-suited for the simple reason that it is the language they are most familiar with. MATLAB programs that use explicit loops to process data one element at a time will run orders of magnitude slower than equivalent compiled programs or a MATLAB program utilizing MATLAB's vector capabilities and/or pre-compiled functions.

To illustrate this point, the following program shows three ways of accomplishing the same basic vector initialization in MATLAB. The first uses the MATLAB colon operator to iterate through values from 1 to 100,000 using a compiled loop within the interpreter. The second uses the pre-compiled linspace() function (which is probably written in C) to accomplish the same iteration. The third uses an explicit interpreted loop, which must be performed by the MATLAB interpreter. Each toc statement reports the time elapsed since the previous tic statement.

fprintf('Colon operator:\n');
tic
list1 = [1:100000];
toc

fprintf('linspace function:\n');
tic
list2 = linspace(1,100000);
toc

fprintf('Explicit loop:\n');
tic
for c = 1:100000
    list3(c) = c;
end
toc
            

Below is the output produced by MATLAB running on a 2GHz Core Duo system:

Colon operator:
Elapsed time is 0.002485 seconds.
linspace function:
Elapsed time is 0.008148 seconds.
Explicit loop:
Elapsed time is 81.661672 seconds.
            

Below is the output produced by Octave running on the same Core Duo system:

Colon operator:
Elapsed time is 0.0000169277 seconds.
linspace function:
Elapsed time is 0.000270128 seconds.
Explicit loop:
Elapsed time is 0.994528 seconds.
            

Note that all three methods do in fact use loops, but the first two use compiled loops while the third uses an interpreted loop. Note also that C code below, using an explicit but compiled loop, executes in about 0.002 seconds, comparable to the MATLAB code without an explicit loop.

size_t  c;
double  list[LIST_SIZE];

for (c=0; c < 100000; ++c)
    list[c] = c;
            

This example illustrates that if MATLAB has no built-in feature such as the vector operation or pre-compiled routine such as linspace() to perform the computations you require, you will need to write code in a compiled language in order to achieve good performance.

One more potential barrier, especially in parallel computing, is the fact that MATLAB is an expensive commercial product and requires a license server in order to run on a cluster. The need to allocate funds, procure licenses, and set up a license server inevitably delays results. In addition, you will not be able to run MATLAB while the license server or the network is down.

One possible solution to this issue is Octave, a free and open-source tool very similar to MATLAB. The Octave interpreter is nearly 100% compatible with MATLAB. In fact, the Octave project considers any differences between MATLAB and Octave to be bugs.

Octave has a graphical user interface (GUI) based on the QT cross-platform GUI tools. The Octave GUI is very similar to MATLAB's GUI, but faster. (The Octave GUI is written in C++, whereas MATLAB's is written in Java.)

The general consensus at the time of this writing is that MATLAB and Octave perform about the same for properly vectorized code, while MATLAB performs better on "bad" code (using explicit loops). Octave is faster at some algorithms and MATLAB is faster at others. In fact, MATLAB and Octave utilize many of the same open source C, C++, and Fortran libraries for many of their common functions. See Table 14.2, “Selection Sort of 100,000 Integers on AMD Opteron” for some performance comparison.

Now let's consider a simple real-world example, namely approximating the area under a curve. There are a variety of approaches, such as the trapezoid rule. This example uses rectangular slices with the height equal to the function value at the center of the slice, which provides a balanced estimate and is very easy to program.

Each of the following examples was run using 10,000,000 slices. The first example program illustrates a typical approach to solving the problem in MATLAB:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Program description:
%
%   Compute the area under the curve of the square root function.
%
%   Demonstrates the value of MATLAB's just-in-time compiler.
%   Unlike the more complex selection sort example, MATLAB is able
%   to compile the explicit loop below at run-time, making it almost
%   as fast as a vector operation.
%
%   History:
%   2013-07-04  Jason Bacon Begin
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic
first_x = 0;
last_x = 1;
slices = 10000000;
slice_width = (last_x - first_x) / slices;

% Use x in the middle of each slice
x = first_x + slice_width / 2;
estimated_area = 0;
while x < last_x
    slice_area = sqrt(x) * slice_width;
    estimated_area = estimated_area + slice_area;
    x = x + slice_width;
end

fprintf('Final x = %0.16f  Estimated area = %0.16f\n', x - slice_width, estimated_area);
toc

Final x = 0.9999999497501700  Estimated area = 0.6666666666390187
Elapsed time is 0.538943 seconds.
            

The example above uses an explicit loop, which could be too slow for a large number of slices, because the MATLAB interpreter has to execute the statements in the loop many times. Below is a vectorized version of the same program:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Program description:
%
%   Compute the area under the curve of the square root function.
%
%   Demonstrates the use of vector operations.  Vector operations are
%   generally much faster than explicit loops, but require much more
%   memory in cases where vectors or matrices are not otherwise needed.
%   This may cause the program to actually run slower than one without
%   vectors, or even fail if the vectors used are too large to fit in
%   memory.
%
%   History:
%   2013-07-04  Jason Bacon Begin
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic
first_x = 0;
last_x = 1;
slices = 10000000;
slice_width = (last_x - first_x) / slices;

% Use x in the middle of each slice
x = linspace(first_x + slice_width / 2, last_x - slice_width / 2, slices);
slice_area = sqrt(x) * slice_width;
estimated_area = sum(slice_area);
fprintf('Final x = %0.16f  Estimated area = %0.16f\n', x(slices), estimated_area);
toc

Final x = 0.9999999500000000  Estimated area = 0.6666666666685899
Elapsed time is 0.434779 seconds.
            

The vectorized version uses no explicit loops, so the MATLAB interpreter will not be a bottleneck. Interestingly, however, it does not run much faster than the non-vectorized code. This is partly because the program is so simple that MATLAB's just-in-time (JIT) compiler can compile the explicit loop for us. More complex programs will not benefit from the JIT compiler as much, however, so explicit loops should be avoided where speed is essential. Speed is also limited by the increase in memory use due to the use of large arrays of 10,000,000 elements. This memory use overwhelms the CPU cache (fast memory that holds only a few megabytes), forcing the program to use the larger, slower main memory.

While using vector operations will generally speed up MATLAB code immensely, the down side is generally an increase in memory use. If we were to increase the number of slices used by this program to 1 billion, the program would now use two vectors of 8 billion bytes (8 gigabytes) each. (Each number in MATLAB program typically requires 8 bytes of memory). If your computer does not have 16 gigabytes of available RAM, the program simply won't run. Even if it does, it will not be able to utilize the CPU cache well, and will take much more time per slice than it does for smaller slice counts.

When programming in MATLAB, we are often faced with a choice between using explicit loops, which are very slow, and using huge amounts of memory, which also slows down the program and potentially prevents it from running at all due to memory exhaustion.

A C version of the program is shown below:

/***************************************************************************
 *  Description:
 *      Compute the area under the curve if the square root function.
 *
 *  Arguments:
 *      None.
 *
 *  Returns:
 *      Standard codes defined by sysexit.h
 *
 *  History: 
 *  Date        Name        Modification
 *  2013-07-04  Jason Bacon Begin
 ***************************************************************************/

#include <stdio.h>
#include <math.h>
#include <sysexits.h>

int     main(int argc,char *argv[])

{
    double  first_x,
            last_x,
            x,
            slice_width,
            estimated_area;
    unsigned long   slices;
            
    first_x = 0.0;
    last_x = 1.0;
    slices = 10000000;
    slice_width = (last_x - first_x) / slices;
    
    // Use x in the middle of each slice
    x = first_x + slice_width / 2;
    estimated_area = 0.0;
    while (x < last_x)
    {
        estimated_area += sqrt(x) * slice_width;
        // Fast, but causes x to drift due to accumulated round off when
        // adding imperfect x + imperfect slice_width
        x = x + slice_width;
    }
    printf("Final x = %0.16f  Estimated area = %0.16f\n",
        x - slice_width, estimated_area);

    return EX_OK;
}

Final x = 0.9999999497501700  Estimated area = 0.6666666666390187
        0.11 real         0.10 user         0.00 sys
            

Sine C is a compiled language, there is no disadvantage to using explicit loops, and no need to use arrays. As a result, this program is much faster than either MATLAB script, and will never be limited by available memory.

Note that this program, as well as the MATLAB program with an explicit loop, has a potential issue with drift caused by round-off error, which is inherent in all computer floating point systems. (Floating point is a storage format similar to scientific notation, which is used by computers to approximate the real number set.)

If the width of a slice cannot be represented perfectly in the computer's floating point format, then each successive value of x will be off by a little more than the previous. For example, representing 1/10 in computer floating point is like trying to represent 1/3 in decimal. It requires an infinite number of digits to represent accurately, so a computer can only approximate it to about 16 decimal digits.

After a large number of slices are processed, the value of x could drift far from the center of the theoretical slice we were aiming for. The MATLAB linspace() function, used in the vectorized example, avoids this issue. The C program below demonstrates how to avoid this problem by basing each value of x on the initial value, instead of the previous value:

/***************************************************************************
 *  Description:
 *      Compute the area under the curve if the square root function.
 *
 *  Arguments:
 *      None.
 *
 *  Returns:
 *      Standard codes defined by sysexit.h
 *
 *  History: 
 *  Date        Name        Modification
 *  2013-07-04  Jason Bacon Begin
 ***************************************************************************/

#include <stdio.h>
#include <math.h>
#include <sysexits.h>

int     main(int argc,char *argv[])

{
    double  first_x,
            last_x,
            base_x,
            x,
            slice_width,
            estimated_area;
    unsigned long
            slices,
            current_slice;

    first_x = 0.0;
    last_x = 1.0;
    slices = 10000000;
    slice_width = (last_x - first_x) / slices;
    
    current_slice = 0;
    estimated_area = 0.0;
    base_x = first_x + slice_width / 2; // Start in middle of first slice
    while (current_slice <= slices)
    {
        // Compute each x directly from a base value to avoid accumulating
        // round-off error.  This is slightly slower than x = x + slice_width
        // due to the multiplication, but produces more accurate results.
        x = base_x + (slice_width * current_slice++);
        estimated_area += sqrt(x) / slice_width;
    }
    printf("Final x = %0.16f  Estimated area = %0.16f\n",
        x - slice_width, estimated_area);

    return EX_OK;
}

This version produces a final value of x that is exactly what it should be.

Final x = 0.9999999500000000  Estimated area = 0.6666667666686176
        0.11 real         0.11 user         0.00 sys
            

Because this version of the program uses an additional multiply operation to compute each value of x, it is slightly slower. It does, however, produce more accurate results, especially for high slice counts.

Perl is an interpreted language with many of the same convenient syntactic features of the Unix shells, plus a wealth of built-in functions such as you would find in the standard C libraries. Perl is designed primarily for text-processing, and is commonly used on WEB servers to process data from WEB forms. It supports a rich set of features for handling regular expressions and is preferred by some people over traditional Unix tools like grep, awk, and sed for many text processing tasks.

Python is another interpreted scripting language, but is popular enough in scientific programming to be worth special mention. Compiled libraries such as SciPy (sigh-pie) and NumPy (num-pie), written in C, C++ and Fortran, make it possible to write Python scripts to do some heavy computation, with the same caveats as MATLAB or Octave. As long as most of the calculations are performed inside the compiled library routines, performance should be acceptable. Iterative code written in Python to do intensive computation, however, will result in run times orders of magnitude slower than equivalent C, C++, or Fortran code.

Performance of raw Python code can be improved using Cython, a compiler for building Python modules from slightly modified Python code, or using Numba, a just-in-time (JIT) compiler for normal Python code. Performance and memory use will not match that of a true compiled language, but will be vastly improved over purely interpreted code. ( See Table 14.2, “Selection Sort of 100,000 Integers on AMD Opteron”. )

Python is becoming a popular alternative to MATLAB for many tasks. While it is not at all similar to MATLAB, it does offer many of the same capabilities. For example, one of the strengths of MATLAB is it's plotting capabilities, which are often used to visualize data. The Matplotlib Python library allows Python scripts to easily create plots, graphs, and charts as well.

One of Python's primary advantages over MATLAB is that it is free, open source, and highly portable. Python scripts can be run on most Unix-compatible systems as well as Windows, at no monetary cost and with minimal installation effort.

The Python community does a good job of encouraging adherence to standards through tools like Distutils and the Python Package Index (PyPI), a website and repository where Python packages are quality-checked before inclusion on the site. This facilitates the development and sharing of software written in Python in much the same way as language-independent package managers. In fact, Python projects that have met the standards to be included on PyPI are almost invariably trivial to add to other package systems like FreeBSD ports and pkgsrc.

R is an open source statistics package with a built-in scripting language. It is a very popular substitute for commercial statistics packages such as SPSS and SAS. Like MATLAB/Octave, the R scripting language is a vector based language and can perform many common vector operations (operations on whole arrays) quickly, because the R interpreter uses compiled loops behind the scenes.

Explicit loops in R scripts are interpreted and extremely slow, and seen in Table 14.2, “Selection Sort of 100,000 Integers on AMD Opteron”. Like MATLAB/Octave, R code must therefore be vectorized as much as possible to achieve best performance. Again, this will likely increase memory use, an unavoidable trade off when using interpreted languages.

The R CRAN project provides add-on R packages that can be easily installed from within R by any user. Many of these packages integrate compiled C, C++, or Fortran code in order to maximize speed and flexibility. Bioconductor is another repository similar to R CRAN that provides many packages for bioinformatics analyses.

R has many frustrating quirks, such as the inability to use certain features within a loop, different behavior when code is run from a script than when entered from the keyboard, etc. With its wealth of statistics libraries, R is an excellent tool for doing quick-and-dirty standard statistical analyses. It is not so convenient to use for developing applications, however, as running those applications will usually require the user to know R programming.

Some developers write much of their code in an interpreted scripting language and only the portions that require speed in a compiled language. If you are using Unix shell scripts, there is no integration necessary. You simply write complete programs in the compiled language and run them from your scripts.

There are also tools available for integrating compiled subprograms into many interpreted languages. The MEX system allows you to write subprograms in C, C++, or Fortran and integrate them into a MATLAB script. Octave also supports MEX as well as its own system MKOCTFILE.

Cython is a Python-to-C translator and module generator. It allows you to write code in an extended Python syntax and compile it for better performance. The Cython compiler generates C or C++ source code that can be compiled to a Python module, which can be imported into Python programs. You may not get the same performance you would from hand-written C, but it will be orders of magnitude faster than interpreted Python.

It's important to consider whether you want to deal with the complexity of integrating multiple languages. Savvy programmers developing complex applications may prefer this approach while less savvy users or those writing simpler programs may find it easier write the entire program in C or Fortran.