Archive for the ‘matlab’ Category

January 28th, 2014

What is Software Carpentry?

Software Carpentry boot-camps aim to ensure that researchers have a working knowledge of several useful technologies from  the world of software development.  Concepts such as version control, unit-testing, task-automation and modular programming  are bread and butter to full time software-developers but are often completely unknown to many researchers; researchers who are nevertheless expected to develop computer code as part of their research output.

An in-depth education in all of these technologies can take a lot of time; time that many researchers simply can’t spare. Fortunately, however, it is possible to learn just enough to completely transform the quality of your workflow in a relatively short amount of time.  Software Carpentry boot camps aim to lay the foundations in around two days.

Software Carpentry…but in MATLAB

Traditionally, these bootcamps are taught using open-source languages such as Python or R but many proprietary languages are also used in academia such as  MATLAB, IDL, Mathematica and Maple to name but a few. In the most recent version of MATLAB, Mathworks introduced a unit testing framework and so now seemed like a good time to try out Software Carpentry in MATLAB.

On Tuesday 14th January I hosted the first ever MATLAB-based Software Carpentry bootcamp in conjunction with The Software Sustainability Institute (of which I am a 2013 fellow) and Mathworks. Held at The University of Manchester, this two-day event gave free software development instruction to researchers from a wide variety of disciplines and, if the feedback forms are to be believed, was a resounding success.

Shoaib Sufi and Aleksandra Pawlik of the Software Sustainability Institute taught material on bash shell scripting and git version control respectively with Ken Deeley of Mathworks providing MATLAB instruction.  I, along with Mathworks’ Jos Martin and Juan Martinez, acted as the less-than-glamourous but ever-helpful classroom assistants.

Fun and games with BYOD – Licensing

For the decade or so that I’ve been teaching programming, I’ve done so in fully equipped teaching labs containing row upon row of identical computers where each have all the required software pre-installed.  For this event, we decided to try a Bring Your Own Device (BYOD) approach..i.e. students bring their own, personal laptops and we provided a list of required software that needed to be installed for the event.  Since Manchester’s MATLAB site license is network-based only and does not allow students to install on personally owned equipment, Mathworks kindly supplied standalone trial licenses for all course attendees.

BYOD has a number of benefits for the student and, in my mind at least, the most important of these benefits is that students are left with a fully working development environment after the course.  This means that they can start applying their new skills immediately after leaving the course which hopefully leads to them being used ‘in anger’ on their research projects.  Since the students were only supplied with trial licenses, this was not to be the case with MATLAB.  They will have to switch to using their on-campus machines or purchase their own copy of ‘standalone’ MATLAB in order to continue working.

Availability of software is not something that usually concerns an organiser of a Software Carpentry boot camp since the likes of Python and R are available everywhere for free.  When using a proprietary language such as MATLAB, however, it’s very much of an issue and I’d advise any future organiser of a MATLAB boot camp to carefully consider this before proceeding.  Of course this isn’t just true of MATLAB, it’s true of any proprietary language one may choose.  It will also be less of a problem if your institution has an all you can eat, ‘Total Academic Headcount’  unlimited site license.

Fun and games with BYOD – switches and glitches

All three major desktop operating systems were represented in the laptops of the 20 or so students — something that occasionally made for fun times.  Here is a list of some of the minor issues that arose over the two days

  • Teaching bash scripting to Windows users always makes me wince a little.  Environments such as Cygwin and git bash are great and a little bit of bash never hurt anyone but I can’t help but wonder if we should be teaching a native scripting language instead such as PowerShell.  After all, Linux users would be surprised if they came to a scripting seminar and we made them use pash!  Of course, if the student ever gets access to a HPC system, it is highly probable that it will be running some variant of Linux and so perhaps everyone should learn at least a little bash.
  • One Mac user had some fancy graphical overlay program which jazzed up his desktop.  Looked great but it turned out that it sometimes caused MATLAB to crash.
  • Some of the MATLAB commands that interacted with the file system worked slightly differently across the three operating systems…something we hadn’t appreciated ahead of time. This caused delays while we figured out a platform-independent way to proceed.
  • Some of the Linux machines exposed a bug in TLS (thread local storage) in Intel’s MKL which caused errors in MATLAB.
  • MATLAB’s keyboard shortcuts are different across operating systems.  It is possible to change the behaviour but MATLABers who’ve been around a while didn’t want to. This sometimes led to confusion if the instructor mentioned an explicit keyboard shortcut that happens to be different on the students machine.

None of these issues were particularly major but the need to consider and resolve them did take up the time of instructors and demonstrators.  In a lab-setting, this extra work would not be necessary.  Obviously we’ll fix some of the above issues in the next iteration of the course but I believe that BYOD sessions will always require a higher number of demonstrators/glamorous assistants than more traditional lab-sessions simply because of the inevitable variation of software and hardware.

No toolboxes….almost!

I can see my Mathworks friends rolling their eyes already but I have to get this off my chest.  Regular readers of this blog will know that I have got some issues with Mathworks’ toolbox system in MATLAB.  I *really* wanted a software carpentry course that only included pure MATLAB–no toolboxes at all.  In the event, we decided on a set of  course materials that required the use of the statistics toolbox because it made certain things so much easier.

For example, MATLAB uses NaNs to represent missing data — a design choice that quickly leads to the desire for functions such as nanmean and nanmedian.  Unfortunately, despite their simplicity, these functions are in the statistics toolbox and so using them leads to less accessible software (since your users need to have that extra toolbox).  Of course you could code up your own versions, or use free implementations easily enough but then you are not using idiomatic MATLAB.  All very frustrating.  This wasn’t the only reason why we added the statistics toolbox to the list of requirements of course but hopefully makes my point.

OK, moan over, I’m done….for now.

Course Material

The main course page is at http://apawlik.github.io/2014-01-14-manchester/ and the course material for version control using git and shell scripting using bash are already available.  We hope to have the MATLAB material (which was developed by Mathworks) available in the near future once we’ve sorted out some legal issues.

The tutorials were very interactive with the instructors demonstrating commands while students followed along on their own machines–barely a slide to be seen, just how I like it.  There were regular exercises with a high ratio of demonstrators to students.  This was very much a hands-on course which, in my opinion, is the best way to learn programming.

One of the advantages of instructor-led courses is that students can go off-piste if they wish and learn all sorts of extra material. One student, for example, asked senior Mathworks developer, Jos Martin, for a quick code-review of her research simulation at the end of the first day.  Listening to Jos’ critique of the code was instructive for several of us.  Jos is also the lead developer of MATLAB’s parallel computing toolbox — something I took advantage of by asking him questions relevant to a code-optimization project I’m currently working on.

It’s very difficult to replicate this kind of interaction in online courses!

Useful teaching technology

The instructors used a couple of pieces of software that I believe significantly enhanced the teaching experience and I plan to use them in future courses of my own.

  • ZoomIt - This is a free presentation tool that allows you to zoom into any area of the screen and subsequently annotate it.  This might not sound like much but it can really improve a presentation.  Here’s a video of some of its functionality.
  • Etherpad – Etherpad is a web based collaborative document editor which we used to post code snippets, links and anything else that anyone felt was useful on the day. We also used it as a chat room which was sometimes useful.
  • PostIt notes – Very low tech but very effective.  Every student had red and green sticky notes.  If a student had no sticky note on the back of their laptop, it meant they were working and would prefer to be left alone.  Green means that they had completed the exercise and red (or in our case, orange) means ‘I want help’.

Feedback and the future

Both Mathworks and the Software Sustainability Institute asked students to fill out their own course feedback forms.  Since the students had had such a great experience, they were all happy to do so but whittling those two feedback forms down to one that satisfies both institutions is something we should definitely aim for!  The feedback was extremely positive with almost everyone agreeing that they had got a lot out of the course.

Of course, it wasn’t perfect and there are several things we could do to make it better for next time:-

  • Some of the ‘Fun and Games with BYOD’ described earlier could have been avoided by modifying the course material a little.
  • The individual sections on bash, git and MATLAB were great but I felt that they needed to be tied together better.  They felt too much like self-contained mini-courses and it wouldn’t take much extra effort on our part to link them together a little more.
  • The course was free and places were strictly limited. A couple of people cancelled less than 24 hours before the first day which didn’t give us enough time to offer the places out to others. Worse still, a few more people simply didn’t bother to turn up and didn’t offer any explanation.  I find that this often happens with free courses and I’ve yet to figure out a way to improve the situation.
  • We need to get all of the course material on line!

All told, the course was a great success and I hope to be running more of them in the future. A huge thanks to the instructors, Shoaib Sufi, Alexandra Pawlik and Ken Deeley who did the lion’s share of the hard work on the day but also to Mathworks’ Jos Martin and Juan Martinez who did a great job of helping out students with the exercises, answering questions and generally ensuring that everything ran as smoothly as possible.  I’d also like to thank Mathworks’ Stasi Revel and Tanya Morton who helped out with sorting out trial licenses and, finally, thanks to Software Carpentry’s Greg Wilson who gave support and advice in the weeks leading to the event.

December 17th, 2013

In a recent Stack Overflow query, someone asked if you could switch off the balancing step when calculating eigenvalues in Python. In the document A case where balancing is harmful, David S. Watkins describes the balancing step as ‘the input matrix A is replaced by a rescaled matrix A* = D-1AD, where D is a diagonal matrix chosen so that, for each i, the ith row and the ith column of A* have roughly the same norm.’  

Such balancing is usually very useful and so is performed by default by software such as  MATLAB or Numpy.  There are times, however, when one would like to switch it off.

In MATLAB, this is easy and the following is taken from the online MATLAB documentation

A = [ 3.0     -2.0      -0.9     2*eps;
     -2.0      4.0       1.0    -eps;
     -eps/4    eps/2    -1.0     0;
     -0.5     -0.5       0.1     1.0];

[VN,DN] = eig(A,'nobalance')
VN =

    0.6153   -0.4176   -0.0000   -0.1528
   -0.7881   -0.3261         0    0.1345
   -0.0000   -0.0000   -0.0000   -0.9781
    0.0189    0.8481   -1.0000    0.0443

DN =
    5.5616         0         0         0
         0    1.4384         0         0
         0         0    1.0000         0
         0         0         0   -1.0000

At the time of writing, it is not possible to directly do this in Numpy (as far as I know at least). Numpy’s eig command currently uses the LAPACK routine DGEEV to do the heavy lifting for double precision matrices.  We can see this by looking at the source code of numpy.linalg.eig where the relevant subsection is

lapack_routine = lapack_lite.dgeev
        wr = zeros((n,), t)
        wi = zeros((n,), t)
        vr = zeros((n, n), t)
        lwork = 1
        work = zeros((lwork,), t)
        results = lapack_routine(_N, _V, n, a, n, wr, wi,
                                  dummy, 1, vr, n, work, -1, 0)
        lwork = int(work[0])
        work = zeros((lwork,), t)
        results = lapack_routine(_N, _V, n, a, n, wr, wi,
                                  dummy, 1, vr, n, work, lwork, 0)

My plan was to figure out how to tell DGEEV not to perform the balancing step and I’d be done. Sadly, however, it turns out that this is not possible. Taking a look at the reference implementation of DGEEV, we can see that the balancing step is always performed and is not user controllable–here’s the relevant bit of Fortran

*     Balance the matrix
*     (Workspace: need N)
*
      IBAL = 1
      CALL DGEBAL( 'B', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )

So, using DGEEV is a dead-end unless we are willing to modifiy and recompile the lapack source — something that’s rarely a good idea in my experience. There is another LAPACK routine that is of use, however, in the form of DGEEVX that allows us to control balancing.  Unfortunately, this routine is not part of the numpy.linalg.lapack_lite interface provided by Numpy and I’ve yet to figure out how to add extra routines to it.

I’ve also discovered that this functionality is an open feature request in Numpy.

Enter the NAG Library

My University has a site license for the commercial Numerical Algorithms Group (NAG) library.  Among other things, NAG offers an interface to all of LAPACK along with an interface to Python.  So, I go through the installation and do

import numpy as np
from ctypes import *
from nag4py.util import Nag_RowMajor,Nag_NoBalancing,Nag_NotLeftVecs,Nag_RightVecs,Nag_RCondEigVecs,Integer,NagError,INIT_FAIL
from nag4py.f08 import nag_dgeevx

eps = np.spacing(1)
np.set_printoptions(precision=4,suppress=True) 

def unbalanced_eig(A):
    """
    Compute the eigenvalues and right eigenvectors of a square array using DGEEVX via the NAG library.
    Requires the NAG C library and NAG's Python wrappers http://www.nag.co.uk/python.asp
    The balancing step that's performed in DGEEV is not performed here.
    As such, this function is the same as the MATLAB command eig(A,'nobalance')

    Parameters
    ----------
    A : (M, M) Numpy array
        A square array of real elements.

        On exit: 
        A is overwritten and contains the real Schur form of the balanced version of the input matrix .

    Returns
    -------
    w : (M,) ndarray
        The eigenvalues
    v : (M, M) ndarray
        The eigenvectors

    Author: Mike Croucher (www.walkingrandomly.com)
    Testing has been mimimal
    """

    order = Nag_RowMajor
    balanc = Nag_NoBalancing
    jobvl = Nag_NotLeftVecs
    jobvr = Nag_RightVecs
    sense = Nag_RCondEigVecs

    n = A.shape[0]
    pda = n
    pdvl = 1

    wr = np.zeros(n)
    wi = np.zeros(n)

    vl=np.zeros(1);
    pdvr = n
    vr = np.zeros(pdvr*n)

    ilo=c_long(0)
    ihi=c_long(0)

    scale = np.zeros(n)
    abnrm = c_double(0)
    rconde = np.zeros(n)
    rcondv = np.zeros(n)

    fail = NagError()
    INIT_FAIL(fail)

    nag_dgeevx(order,balanc,jobvl,jobvr,sense, 
               n, A.ctypes.data_as(POINTER(c_double)), pda, wr.ctypes.data_as(POINTER(c_double)),
               wi.ctypes.data_as(POINTER(c_double)),vl.ctypes.data_as(POINTER(c_double)),pdvl, 
               vr.ctypes.data_as(POINTER(c_double)),pdvr,ilo,ihi, scale.ctypes.data_as(POINTER(c_double)),
               abnrm, rconde.ctypes.data_as(POINTER(c_double)),rcondv.ctypes.data_as(POINTER(c_double)),fail)

    if all(wi == 0.0):
            w = wr
            v = vr.reshape(n,n)
    else:
            w = wr+1j*wi
            v = array(vr, w.dtype).reshape(n,n)

    return(w,v)

Define a test matrix:

A = np.array([[3.0,-2.0,-0.9,2*eps],
          [-2.0,4.0,1.0,-eps],
          [-eps/4,eps/2,-1.0,0],
          [-0.5,-0.5,0.1,1.0]])

Do the calculation

(w,v) = unbalanced_eig(A)

which gives

(array([ 5.5616,  1.4384,  1.    , -1.    ]),
 array([[ 0.6153, -0.4176, -0.    , -0.1528],
       [-0.7881, -0.3261,  0.    ,  0.1345],
       [-0.    , -0.    , -0.    , -0.9781],
       [ 0.0189,  0.8481, -1.    ,  0.0443]]))

This is exactly what you get by running the MATLAB command eig(A,’nobalance’).

Note that unbalanced_eig(A) changes the input matrix A to

array([[ 5.5616, -0.0662,  0.0571,  1.3399],
       [ 0.    ,  1.4384,  0.7017, -0.1561],
       [ 0.    ,  0.    ,  1.    , -0.0132],
       [ 0.    ,  0.    ,  0.    , -1.    ]])

According to the NAG documentation, this is the real Schur form of the balanced version of the input matrix.  I can’t see how to ask NAG to not do this. I guess that if it’s not what you want unbalanced_eig() to do,  you’ll need to pass a copy of the input matrix to NAG.

The IPython notebook

The code for this article is available as an IPython Notebook

The future

This blog post was written using Numpy version 1.7.1. There is an enhancement request for the functionality discussed in this article open in Numpy’s git repo and so I expect this article to become redundant pretty soon.

December 6th, 2013

A question I get asked a lot is ‘How can I do nonlinear least squares curve fitting in X?’ where X might be MATLAB, Mathematica or a whole host of alternatives.  Since this is such a common query, I thought I’d write up how to do it for a very simple problem in several systems that I’m interested in

This is the MATLAB version. For other versions,see the list below

The problem

xdata = -2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9
ydata = 0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001

and you’d like to fit the function

 F(p_1,p_2,x) = p_1 \cos(p_2 x)+p_2 \sin(p_1 x)

using nonlinear least squares.  You’re starting guesses for the parameters are p1=1 and P2=0.2

For now, we are primarily interested in the following results:

  • The fit parameters
  • Sum of squared residuals

Future updates of these posts will show how to get other results such as confidence intervals. Let me know what you are most interested in.

How you proceed depends on which toolboxes you have. Contrary to popular belief, you don’t need the Curve Fitting toolbox to do curve fitting…particularly when the fit in question is as basic as this. Out of the 90+ toolboxes sold by The Mathworks, I’ve only been able to look through the subset I have access to so I may have missed some alternative solutions.

Pure MATLAB solution (No toolboxes)

In order to perform nonlinear least squares curve fitting, you need to minimise the squares of the residuals. This means you need a minimisation routine.  Basic MATLAB comes with the fminsearch function which is based on the Nelder-Mead simplex method.  For this particular problem, it works OK but will not be suitable for more complex fitting problems.  Here’s the code

format compact
format long
xdata = [-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9];
ydata = [0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001];

%Function to calculate the sum of residuals for a given p1 and p2
fun = @(p) sum((ydata - (p(1)*cos(p(2)*xdata)+p(2)*sin(p(1)*xdata))).^2);

%starting guess
pguess = [1,0.2];

%optimise
[p,fminres] = fminsearch(fun,pguess)

This gives the following results

p =
1.881831115804464 0.700242006994123
fminres =
0.053812720914713

All we get here are the parameters and the sum of squares of the residuals.  If you want more information such as 95% confidence intervals, you’ll have a lot more hand-coding to do. Although fminsearch works fine in this instance, it soon runs out of steam for more complex problems.

MATLAB with optimisation toolbox

With respect to this problem, the optimisation toolbox gives you two main advantages over pure MATLAB.  The first is that better optimisation routines are available so more complex problems (such as those with constraints) can be solved and in less time.  The second is the provision of the lsqcurvefit function which is specifically designed to solve curve fitting problems.  Here’s the code

format long
format compact
xdata = [-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9];
ydata = [0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001];

%Note that we don't have to explicitly calculate residuals
fun = @(p,xdata) p(1)*cos(p(2)*xdata)+p(2)*sin(p(1)*xdata);

%starting guess
pguess = [1,0.2];

[p,fminres] = lsqcurvefit(fun,pguess,xdata,ydata)

This gives the results

p =
   1.881848414551983   0.700229137656802
fminres =
   0.053812696487326

MATLAB with statistics toolbox

There are two interfaces I know of in the stats toolbox and both of them give a lot of information about the fit. The problem set up is the same in both cases

%set up for both fit commands in the stats toolbox
xdata = [-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9];
ydata = [0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001];

fun = @(p,xdata) p(1)*cos(p(2)*xdata)+p(2)*sin(p(1)*xdata);

pguess = [1,0.2];

Method 1 makes use of NonLinearModel.fit

mdl = NonLinearModel.fit(xdata,ydata,fun,pguess)

The returned object is a NonLinearModel object

class(mdl)
ans =
NonLinearModel

which contains all sorts of useful stuff

mdl = 
Nonlinear regression model:
    y ~ p1*cos(p2*xdata) + p2*sin(p1*xdata)

Estimated Coefficients:
          Estimate             SE                 
    p1      1.8818508110535      0.027430139389359
    p2    0.700229815076442    0.00915260662357553

          tStat               pValue              
    p1    68.6052223191956    2.26832562501304e-12
    p2    76.5060538352836    9.49546284187105e-13

Number of observations: 10, Error degrees of freedom: 8
Root Mean Squared Error: 0.082
R-Squared: 0.996,  Adjusted R-Squared 0.995
F-statistic vs. zero model: 1.43e+03, p-value = 6.04e-11

If you don’t need such heavyweight infrastructure, you can make use of the statistic toolbox’s nlinfit function

[p,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(xdata,ydata,fun,pguess);

Along with our parameters (p) this also provides the residuals (R), Jacobian (J), Estimated variance-covariance matrix (CovB), Mean Squared Error (MSE) and a structure containing information about the error model fit (ErrorModelInfo).

Both nlinfit and NonLinearModel.fit use the same minimisation algorithm which is based on Levenberg-Marquardt

MATLAB with Symbolic Toolbox

MATLAB’s symbolic toolbox provides a completely separate computer algebra system called Mupad which can handle nonlinear least squares fitting via its stats::reg function.  Here’s how to solve our problem in this environment.

First, you need to start up the mupad application.  You can do this by entering

mupad

into MATLAB. Once you are in mupad, the code looks like this

xdata := [-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9]:
ydata := [0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001]:
stats::reg(xdata,ydata,p1*cos(p2*x)+p2*sin(p1*x),[x],[p1,p2],StartingValues=[1,0.2])

The result returned is

[[1.88185085, 0.7002298172], 0.05381269642]

These are our fitted parameters,p1 and p2, along with the sum of squared residuals. The documentation tells us that the optimisation algorithm is Levenberg-Marquardt– this is rather better than the simplex algorithm used by basic MATLAB’s fminsearch.

MATLAB with the NAG Toolbox

The NAG Toolbox for MATLAB is a commercial product offered by the UK based Numerical Algorithms Group.  Their main products are their C and Fortran libraries but they also have a comprehensive MATLAB toolbox that contains something like 1500+ functions.  My University has a site license for pretty much everything they put out and we make great use of it all.  One of the benefits of the NAG toolbox over those offered by The Mathworks is speed.  NAG is often (but not always) faster since its based on highly optimized, compiled Fortran.  One of the problems with the NAG toolbox is that it is difficult to use compared to Mathworks toolboxes.

In an earlier blog post, I showed how to create wrappers for the NAG toolbox to create an easy to use interface for basic nonlinear curve fitting.  Here’s how to solve our problem using those wrappers.

format long
format compact
xdata = [-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9];
ydata = [0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001];

%Note that we don't have to explicitly calculate residuals
fun = @(p,xdata) p(1)*cos(p(2)*xdata)+p(2)*sin(p(1)*xdata);
start = [1,0.2];
[p,fminres]=nag_lsqcurvefit(fun,start,xdata,ydata)

which gives

Warning: nag_opt_lsq_uncon_mod_func_easy (e04fy) returned a warning indicator (5) 
p =
   1.881850904268710   0.700229557886739
fminres =
   0.053812696425390

For convenience, here’s the two files you’ll need to run the above (you’ll also need the NAG Toolbox for MATLAB of course)

MATLAB with curve fitting toolbox

One would expect the curve fitting toolbox to be able to fit such a simple curve and one would be right :)

xdata = [-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9];
ydata = [0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001];

opt = fitoptions('Method','NonlinearLeastSquares',...
               'Startpoint',[1,0.2]);
f = fittype('p1*cos(p2*x)+p2*sin(p1*x)','options',opt);

fitobject = fit(xdata',ydata',f)

Note that, unlike every other Mathworks method shown here, xdata and ydata have to be column vectors. The result looks like this

fitobject = 

     General model:
     fitobject(x) = p1*cos(p2*x)+p2*sin(p1*x)
     Coefficients (with 95% confidence bounds):
       p1 =       1.882  (1.819, 1.945)
       p2 =      0.7002  (0.6791, 0.7213)

fitobject is of type cfit:

 class(fitobject)

ans =
cfit

In this case it contains two fields, p1 and p2, which are the parameters we are looking for

>> fieldnames(fitobject)
ans = 
    'p1'
    'p2'
>> fitobject.p1
ans =
   1.881848414551983
>> fitobject.p2
ans =
    0.700229137656802

For maximum information, call the fit command like this:

[fitobject,gof,output] = fit(xdata',ydata',f)

fitobject = 

     General model:
     fitobject(x) = p1*cos(p2*x)+p2*sin(p1*x)
     Coefficients (with 95% confidence bounds):
       p1 =       1.882  (1.819, 1.945)
       p2 =      0.7002  (0.6791, 0.7213)

gof = 

           sse: 0.053812696487326
       rsquare: 0.995722238905101
           dfe: 8
    adjrsquare: 0.995187518768239
          rmse: 0.082015773244637

output = 

           numobs: 10
         numparam: 2
        residuals: [10x1 double]
         Jacobian: [10x2 double]
         exitflag: 3
    firstorderopt: 3.582047395989108e-05
       iterations: 6
        funcCount: 21
     cgiterations: 0
        algorithm: 'trust-region-reflective'
          message: [1x86 char]
October 10th, 2013

From the wikipedia page on Division by Zero: “The IEEE 754 standard specifies that every floating point arithmetic operation, including division by zero, has a well-defined result”.

MATLAB supports this fully:

>> 1/0
ans =
   Inf
>> 1/(-0)
ans =
  -Inf
>> 0/0
ans =
   NaN

Julia is almost there, but doesn’t handled the signed 0 correctly (This is using Version 0.2.0-prerelease+3768 on Windows)

julia> 1/0
Inf

julia> 1/(-0)
Inf

julia> 0/0
NaN

Python throws an exception. (Python 2.7.5 using IPython shell)

In [4]: 1.0/0.0
---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
 in ()
----> 1 1.0/0.0

ZeroDivisionError: float division by zero

In [5]: 1.0/(-0.0)
---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
 in ()
----> 1 1.0/(-0.0)

ZeroDivisionError: float division by zero

In [6]: 0.0/0.0
---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
 in ()
----> 1 0.0/0.0

ZeroDivisionError: float division by zero

Update:
Julia does do things correctly, provided I make it clear that I am working with floating point numbers:

julia> 1.0/0.0
Inf

julia> 1.0/(-0.0)
-Inf

julia> 0.0/0.0
NaN
September 25th, 2013

I support scientific applications at The University of Manchester (see my LinkedIn profile if you’re interested in the details) and part of my job involves working on code written by researchers in a variety of languages.  When I say ‘variety’ I really mean it – MATLAB, Mathematica, Python, C, Fortran, Julia, Maple, Visual Basic and PowerShell are some languages I’ve worked with this month for instance.

Having to juggle the semantics of so many languages in my head sometimes leads to momentary confusion when working on someone’s program.  For example, I’ve been doing a lot of Python work recently but this morning I was hacking away on someone’s MATLAB code.  Buried deep within the program, it would have been very sensible to be able to do the equivalent of this:

a=rand(3,3)

a =
    0.8147    0.9134    0.2785
    0.9058    0.6324    0.5469
    0.1270    0.0975    0.9575

>> [x,y,z]=a(:,1)

Indexing cannot yield multiple results.

That is, I want to be able to take the first column of the matrix a and broadcast it out to the variables x,y and z. The code I’m working on uses MUCH bigger matrices and this kind of assignment is occasionally useful since the variable names x,y,z have slightly more meaning than a(1,3), a(2,3), a(3,3).

The only concise way I’ve been able to do something like this using native MATLAB commands is to first convert to a cell. In MATLAB 2013a for instance:

>> temp=num2cell(a(:,1));
>> [x y z] = temp{:}

x =
    0.8147

y =
    0.9058

z =
    0.1270

This works but I think it looks ugly and introduces conversion overheads. The problem I had for a short time is that I subconsciously expected multiple assignment to ‘Just Work’ in MATLAB since the concept makes sense in several other languages I use regularly.

Python:

from pylab import rand
a=rand(3,3)
[a,b,c]=a[:,0]

Mathematica:

a = RandomReal[1, {3, 3}]
{x,y,z}=a[[All,1]]

Julia:

a=rand(3,3);
(x,y,z)=a[:,1]

I’ll admit that I don’t often need this construct in MATLAB but it would definitely be occasionally useful. I wonder what other opinions there are out there? Do you think multiple assignment is useful (in any language)?

September 16th, 2013

Last week I gave a live demo of the IPython notebook to a group of numerical analysts and one of the computations we attempted to do was to solve the following linear system using Numpy’s solve command.

 \left( \begin{array}{ccc} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{array} \right) x = \left( \begin{array}{c} 15 \\ 15 \\ 15 \\ \end{array} \right)

Now, the matrix shown above is singular and so we expect that we might have problems. Before looking at how Numpy deals with this computation, lets take a look at what happens if you ask MATLAB to do it

>> A=[1 2 3;4 5 6;7 8 9];
>> b=[15;15;15];
>> x=A\b
Warning: Matrix is close to singular
or badly scaled. Results may be
inaccurate. RCOND =  1.541976e-18. 
x =
  -39.0000
   63.0000
  -24.0000

MATLAB gives us a warning that the input matrix is close to being singular (note that it didn’t actually recognize that it is singular) along with an estimate of the reciprocal of the condition number. It tells us that the results may be inaccurate and we’d do well to check. So, lets check:

>> A*x

ans =
   15.0000
   15.0000
   15.0000

>> norm(A*x-b)

ans =
2.8422e-14

We seem to have dodged the bullet since, despite the singular nature of our matrix, MATLAB has able to find a valid solution. MATLAB was right to have warned us though…in other cases we might not have been so lucky.

Let’s see how Numpy deals with this using the IPython notebook:

In [1]:
import numpy
from numpy import array
from numpy.linalg import solve
A=array([[1,2,3],[4,5,6],[7,8,9]])
b=array([15,15,15])
solve(A,b)

Out[1]:

array([-39.,  63., -24.])

It gave the same result as MATLAB [See note 1], presumably because it’s using the exact same LAPACK routine, but there was no warning of the singular nature of the matrix.  During my demo, it was generally felt by everyone in the room that a warning should have been given, particularly when working in an interactive setting.

If you look at the documentation for Numpy’s solve command you’ll see that it is supposed to throw an exception when the matrix is singular but it clearly didn’t do so here. The exception is sometimes thrown though:

In [4]:

C=array([[1,1,1],[1,1,1],[1,1,1]])
x=solve(C,b)

---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
 in ()
      1 C=array([[1,1,1],[1,1,1],[1,1,1]])
----> 2 x=solve(C,b)

C:\Python32\lib\site-packages\numpy\linalg\linalg.py in solve(a, b)
    326     results = lapack_routine(n_eq, n_rhs, a, n_eq, pivots, b, n_eq, 0)
    327     if results['info'] > 0:
--> 328         raise LinAlgError('Singular matrix')
    329     if one_eq:
    330         return wrap(b.ravel().astype(result_t))

LinAlgError: Singular matrix

It seems that Numpy is somehow checking for exact singularity but this will rarely be detected due to rounding errors. Those I’ve spoken to consider that MATLAB’s approach of estimating the condition number and warning when that is high would be better behavior since it alerts the user to the fact that the matrix is badly conditioned.

Thanks to Nick Higham and David Silvester for useful discussions regarding this post.

Notes
[1] – The results really are identical which you can see by rerunning the calculation after evaluating format long in MATLAB and numpy.set_printoptions(precision=15) in Python

August 5th, 2013

While working on someone’s MATLAB code today there came a point when it was necessary to generate a vector of powers.  For example, [a a^2 a^3....a^10000] where a=0.999

a=0.9999;
y=a.^(1:10000);

This isn’t the only way one could form such a vector and I was curious whether or not an alternative method might be faster. On my current machine we have:

>> tic;y=a.^(1:10000);toc
Elapsed time is 0.001526 seconds.
>> tic;y=a.^(1:10000);toc
Elapsed time is 0.001529 seconds.
>> tic;y=a.^(1:10000);toc
Elapsed time is 0.001716 seconds.

Let’s look at the last result in the vector y

>> y(end)
ans =
   0.367861046432970

So, 0.0015-ish seconds to beat.

>> tic;y1=cumprod(ones(1,10000)*a);toc
Elapsed time is 0.000075 seconds.
>> tic;y1=cumprod(ones(1,10000)*a);toc
Elapsed time is 0.000075 seconds.
>> tic;y1=cumprod(ones(1,10000)*a);toc
Elapsed time is 0.000075 seconds.

soundly beaten! More than a factor of 20 in fact. Let’s check out that last result

>> y1(end)
ans =
   0.367861046432969

Only a difference in the 15th decimal place–I’m happy with that. What I’m wondering now, however, is will my faster method ever cause me grief?

This is only an academic exercise since this is not exactly a hot spot in the code!

May 3rd, 2013

I recently picked up a few control theory books from the University library to support a project I am involved with right now and was interested in the seemingly total dominance of MATLAB in this subject area.  Since I’m not an expert in control systems, I’m not sure if this is because MATLAB is genuinely the best tool for the job or if it’s simply because it’s been around for a very long time and so has become entrenched.  Comments from anyone who works in relevant fields would be most welcome.

On its own, MATLAB is insufficient to teach introductory control systems courses — you also need the control systems toolbox as a bare minimum but most books and courses also seem to require Simulink and the symbolic math toolbox.  All of these are included in the student edition of MATLAB which is very reasonably priced.

If you are not a registered student, however, and don’t work for someone who can provide you with MATLAB it’s going to be very expensive!  As far as I can tell, your only option would be to purchase commercial licenses which are very expensive (as in thousands of dollars/pounds for MATLAB and a few toolboxes).

What else is out there?

I have a strong interest in mathematical software and so I know that there are several products that have support for control theory. Here are some that I know of and have access to myself

  • Mathematica – Its symbolic math support far exceeds that of MATLAB and it is on an equal footing numerically but its control systems support is much more recent and I don’t know of a textbook that utilizes it.  One benefit of Mathematica is that it doesn’t separate functionality out into toolboxes – everything is just built in.  Another benefit to tinkerers is the home edition which gives you the full product at a much lower price than commercial licenses.
  • Maple – This also has very strong symbolic and numeric math support.  It also comes with some Control Systems support built in.  Like Mathematica, it has a home edition for non-commercial tinkering and learning.
  • Labview - A graphical programming language that I’m only just starting to get used to.  It has lots of users and advocates in my employers electrical and mechanical engineering departments.  There is no support for symbolic computing as far as I know.
  • Python – Python is a superb general purpose scripting language that’s also completely free.  Numerics are taken care of by Numpy, symbolics by Sympy and there is a control theory module, the development of which is coordinated by Richard Murray of Caltech (The same Richard Murray that co-wrote the book Feedback Systems: An Introduction for Scientists and Engineers).
  • Octave – Octave is a free implementation of the MATLAB .m language.  It also has a free control package.
  • Scilab – Scilab is a free numerical environment that also has a free control package.

I haven’t mentioned Simulink alternatives in this post since I’ve discussed them before.

Questions

Some questions that arise are

  • Are there any other alternatives to those listed above?
  • Do these alternatives have sufficient functionality to support undergraduate courses in control systems and control theory?
  • What would be the best language to use if you were teaching control systems as a Massively Open Online Course (MOOC)?
  • Does it matter to employers which computational language you learned your control systems in as an undergraduate?

I find that the final point is very divisive among people I discuss it with.  On the one hand you have those who say ‘It’s the concepts that matter, the language you choose to implement them in is much less important’ and on the other hand you have those who say ‘It’s gotta be MATLAB, my father used MATLAB and his grandfather before him. Industry uses MATLAB, I only know MATLAB, we must teach MATLAB.’

April 23rd, 2013

I was recently working on some MATLAB code with Manchester University’s David McCormick.  Buried deep within this code was a function that was called many,many times…taking up a significant amount of overall run time.  We managed to speed up an important part of this function by almost a factor of two (on his machine) simply by inserting two brackets….a new personal record in overall application performance improvement per number of keystrokes.

The code in question is hugely complex, but the trick we used is really very simple.  Consider the following MATLAB code

>> a=rand(4000);
>> c=12.3;
>> tic;res1=c*a*a';toc
Elapsed time is 1.472930 seconds.

With the insertion of just two brackets, this runs quite a bit faster on my Ivy Bridge quad-core desktop.

>> tic;res2=c*(a*a');toc
Elapsed time is 0.907086 seconds.

So, what’s going on? Well, we think that in the first version of the code, MATLAB first calculates c*a to form a temporary matrix (let’s call it temp here) and then goes on to find temp*a’.  However, in the second version, we think that MATLAB calculates a*a’ first and in doing so it takes advantage of the fact that the result of multiplying a matrix by its transpose will be symmetric which is where we get the speedup.

Another demonstration of this phenomena can be seen as follows

>> a=rand(4000);
>> b=rand(4000);
>> tic;a*a';toc 
Elapsed time is 0.887524 seconds.
>> tic;a*b;toc  
Elapsed time is 1.473208 seconds.
>> tic;b*b';toc
Elapsed time is 0.966085 seconds.

Note that the symmetric matrix-matrix multiplications are faster than the general, non-symmetric one.

March 10th, 2013

Ever since I took a look at GPU accelerating simple Monte Carlo Simulations using MATLAB, I’ve been disappointed with the performance of its GPU random number generator. In MATLAB 2012a, for example, it’s not much faster than the CPU implementation on my GPU hardware.  Consider the following code

function gpuRandTest2012a(n)

mydev=gpuDevice();
disp('CPU - Mersenne Twister');
tic
CPU = rand(n);
toc

sg = parallel.gpu.RandStream('mrg32k3a','Seed',1);
parallel.gpu.RandStream.setGlobalStream(sg);
disp('GPU - mrg32k3a');
tic
Rg = parallel.gpu.GPUArray.rand(n);
wait(mydev);
toc

Running this on MATLAB 2012a on my laptop gives me the following typical times (If you try this out yourself, the first run will always be slower for various reasons I’ll not go into here)

>> gpuRandTest2012a(10000)
CPU - Mersenne Twister
Elapsed time is 1.330505 seconds.
GPU - mrg32k3a
Elapsed time is 1.006842 seconds.

Running the same code on MATLAB 2012b, however, gives a very pleasant surprise with typical run times looking like this

CPU - Mersenne Twister
Elapsed time is 1.590764 seconds.
GPU - mrg32k3a
Elapsed time is 0.185686 seconds.

So, generation of random numbers using the GPU is now over 7 times faster than CPU generation on my laptop hardware–a significant improvment on the previous implementation.

New generators in 2012b

The MATLAB developers went a little further in 2012b though.  Not only have they significantly improved performance of the mrg32k3a combined multiple recursive generator, they have also implemented two new GPU random number generators based on the Random123 library.  Here are the timings for the generation of 100 million random numbers in MATLAB 2012b

CPU - Mersenne Twister
Elapsed time is 1.370252 seconds.
GPU - mrg32k3a
Elapsed time is 0.186152 seconds.
GPU - Threefry4x64-20
Elapsed time is 0.145144 seconds.
GPU - Philox4x32-10
Elapsed time is 0.129030 seconds.

Bear in mind that I am running this on the relatively weak GPU of my laptop!  If anyone runs it on something stronger, I’d love to hear of your results.

  • Laptop model: Dell XPS L702X
  • CPU: Intel Core i7-2630QM @2Ghz software overclockable to 2.9Ghz. 4 physical cores but total 8 virtual cores due to Hyperthreading.
  • GPU: GeForce GT 555M with 144 CUDA Cores.  Graphics clock: 590Mhz.  Processor Clock:1180 Mhz. 3072 Mb DDR3 Memeory
  • RAM: 8 Gb
  • OS: Windows 7 Home Premium 64 bit.
  • MATLAB: 2012a/2012b