Contents

Problem Statement

As part of program, suppose we need to see if some value lies between two others. let's say A is the value we are checking at the limits are low=3 and high=7.

low = 3;
high = 7;

Mathematically you might write this as

  low < A < high

Let's try that for A values both inside the range and outside to start. And let me place the expression into an anonymous function so I don't have to keep repeating it.

myExpr = @(x) low < x < high;
inResult = myExpr(pi)
outResult = myExpr(17)

inResult =
     1
outResult =
     1

WHAT???

It can't be true that both pi and 17 lie between 3 and 7. So what's going on? Let's dissect the expression.

First Part of Expression: low < A

Let's look at the first part of the expression.

step1In = low < pi
step1Out = low < 17

step1In =
     1
step1Out =
     1

and we see that this is true for both of our inputs.

Second Part of Expression: previous output < high

The second part of our expression uses the output from the first expression and continues from there.

step2In = step1In < high
step2Out = step1Out < high

step2In =
     1
step2Out =
     1

and we see that we get ones, or true, for both of these. What's going on?

Look at the Types

whos step1*

  Name          Size            Bytes  Class      Attributes

  step1In       1x1                 1  logical              
  step1Out      1x1                 1  logical              

The results from step 1 are logical - are these numbers greater than low? And the answers for both of our values is yes, or true, represented in MATLAB as logical values. When we take these values as inputs in the second step, what happens is the true values are interpreted as numeric inputs with value 1. And then we ask if 1 is less than high. Which it is in both of these cases!

How to Get the Expected Answer

How do we get the expected answer, and it's easy. We simply combine two logical expressions, but in a different way than above.

myExprCorrect = @(x) (low < x) & (x < high)
inResult1 = myExprCorrect(pi)
outResult1 = myExprCorrect(17)

myExprCorrect = 
    @(x)(low<x)&(x<high)
inResult1 =
     1
outResult1 =
     0

What we did is checked first to see if the number was greater than low and separately checked the same number with high. After getting two logical answers, we combine them. They must both be true for numbers that lie between low and high and hence the result should yield true only under those conditions.

For what it's worth, I always use parentheses to group my expressions to make them very readable for me so I don't need to wonder later what I intended to be testing.

Have Compound Expressions Caused You Problems?

Let me know here if you've had trouble with expressions like the one in this post.

Get the MATLAB code

Published with MATLAB® 7.14

Contents

Visualizing Quiksort

With much amusement, many of us from MathWorks recently viewed this video depicting the quiksort algorithm in dance. I don't think you'll see a MATLAB animation quite like this!

Visualizing Other Sort Algorithms

If you are interested in similar visual depictions of other sorting algorithms, check out this link.

A Place Where Math Meets the Humanities

In addition to these animations showing off mathematical algorithms, I have also enjoyed listening to and trying my hand at change ringing. I find similar points of intersection between my love of math and my interest in other art forms as well. Also check out the 17 regular tilings of the plane.

Do you have similar links between some interesting math and art? Let me know here.

Get the MATLAB code

Published with MATLAB® 7.14

Contents

Running the Blackjack Simulation in MATLAB

The calculation of the result of playing a single hand of blackjack was written by Cleve Moler, and is available on the file exchange. During this article, we're going to use a slightly modified version of that program where we've extracted just the calculation of a single hand of blackjack so that we can invoke it in a variety of different ways. There's a discussion of the algorithm in chapter 9 of Cleve's book "Numerical Computing with MATLAB".

Cleve describes the amount won playing a single hand of blackjack as follows:

"[win] zero for a 'push', win $15 for a blackjack, win or lose $10 on a hand that has not been split or doubled, win or lose $20 on hands that have been split or doubled once, and win or lose $30 or $40 on hands that have been doubled after a split. The $30 and $40 payoffs occur rarely (and may not be allowed at some casinos)."

We can calculate the winnings of playing several individual hands of blackjack simply by invoking blackjack_kernel a few times.

winnings = zeros(1,10);
for ii=1:10
    winnings(ii) = blackjack_kernel;
end
winnings

winnings =

    15   -10     0     0   -10    10   -10   -10   -10    10

Running Blackjack as a Script Directly in MATLAB

The script blackjack_script simulates a number of players playing many hands of blackjack, and computes for each player their net profit after those hands. Let's run the script and then plot the results as a histogram showing the distribution of profit and loss among the players.

type blackjack_script

blackjack_script;
hist(results);
title('Profit and loss playing blackjack');
xlabel('Profit'); ylabel('Count');

numHands   = 500;
numPlayers = 50;
results    = zeros(1, numPlayers);
tic
for h = 1:numHands
    thisHand = zeros(1, numPlayers);
    for p = 1:numPlayers
        thisHand(p) = blackjack_kernel;
    end
    results = results + thisHand;
end
t = toc;
fprintf('Time to simulate %d players playing %d hands of blackjack: %.3g seconds\n', ...
        numPlayers, numHands, t);

Time to simulate 50 players playing 500 hands of blackjack: 8.04 seconds

Submitting Batch Jobs to Run Blackjack on the "local" Cluster

Now that our script works correctly, we can submit the script for execution on a cluster. In this case, we're going to use the "local" cluster since it's always available, but you might submit the script for execution on a remote cluster if you have one available. When you do submit to a remote cluster, you're free to close down your MATLAB session and collect the results later. If you're using the local cluster, you can carry on using your MATLAB session while the batch job runs in the background; also, you can have several batch jobs running simultaneously.

We invoke the batch command with the name of the script we wish to run, and here we're also specifying which 'Profile' we wish to use. Profiles identify connections you have made to different clusters that might be available, the 'local' profile is always available. Other cluster types require an MDCS installation. If you only have a single Profile, you don't need to specify the Profile argument. The return from the batch command is a job object. The job object allows us to track the progress of execution on the cluster, and when it has finished, we can access the results.

As well as using the batch command directly from the MATLAB command line, you can right-click on a script file in MATLAB's Current Folder browser, and select "Run Script as Batch Job". This uses your default cluster profile.

In this example, the first thing that we do with the job is call wait so that we only proceed when the job has finished executing. Alternatively, you can carry on working in MATLAB and check the state of the job by looking at job.State, loading the results when the state is 'finished'.

job = batch('blackjack_script', 'Profile', 'local');

% Wait for the job to complete execution.
wait(job);

% Load the results of executing the job into the MATLAB workspace
load(job);

% And display the diary output from the job
disp('Diary output from the batch job:')
diary(job);

% We can plot the results as we did previously.
hist(results);
title('Profit and loss playing blackjack');
xlabel('Profit'); ylabel('Count');

Diary output from the batch job:
Time to simulate 50 players playing 500 hands of blackjack: 7.93 seconds

Adding parfor Loops to the Script

Running a batch job on the local machine doesn't make things go any quicker, although running a job on a remote cluster might give you access to a more powerful machine. Even on the local machine, we can still use parfor to gain a speedup. Each iteration of the outer for loop in blackjack_script is independent, and the calculation of results is a "reduction" or summary operation that parfor can understand. So, we can recast the script into a variant using parfor simply by changing the outer loop:

type blackjack_script_parfor

numHands   = 500;
numPlayers = 50;
results    = zeros(1, numPlayers);
tic
parfor h = 1:numHands
    thisHand = zeros(1, numPlayers);
    for p = 1:numPlayers
        thisHand(p) = blackjack_kernel;
    end
    results = results + thisHand;
end
t = toc;
fprintf('Time to simulate %d players playing %d hands of blackjack: %.3g seconds\n', ...
        numPlayers, numHands, t);

Running the parfor script with an open MATLAB pool

We can check that the parfor loop is working correctly by opening an interactive matlabpool session. An interactive matlabpool session launches workers based on the Profile you specify, and makes them available to your MATLAB session for running the body of parfor loops. Normally, it's a good idea to leave matlabpool open while you're working, but we're going to close it just after running the script to free up workers so that we can run a batch job.

% Open matlabpool using the 'local' profile
matlabpool open local

% Run our script
blackjack_script_parfor

% Close the pool again to free up workers for when we run the batch job
matlabpool close

Starting matlabpool using the 'local' profile ... connected to 6 labs.
Time to simulate 50 players playing 500 hands of blackjack: 1.87 seconds
Sending a stop signal to all the labs ... stopped.

Running a Batch Job Containing a parfor Loop

To run a batch job where a parfor loop is used, we need to specify an additional argument to the batch command - 'Matlabpool'. The value of the option specifies how many additional workers to use when running the script. The machine I'm using has 6 cores, so I can run 1 worker to control the script execution, and 5 additional workers to act as a matlabpool. These additional workers operate on the body of the parfor loop in parallel to calculate the result more quickly.

job = batch('blackjack_script_parfor', ... % Run the script containing parfor
            'Profile', 'local', ...        % using the 'local' profile
            'Matlabpool', 5);              % with a Matlabpool of size 5

% Wait for the job to complete execution.
wait(job);

% Load the results of executing the job into the MATLAB workspace
load(job);

% And display the diary output from the job
disp('Diary output from the batch job using parfor:')
diary(job);

% We can plot the results as we did previously.
hist(results);
title('Profit and loss playing blackjack');
xlabel('Profit'); ylabel('Count');

Diary output from the batch job using parfor:
Time to simulate 50 players playing 500 hands of blackjack: 1.92 seconds

Conclusions

We have seen how using the batch command, we can send a script to be executed on a cluster. We've used the 'local' cluster here, which is a part of Parallel Computing Toolbox, but we could have just as easily used a remote cluster if a MATLAB Distributed Computing Server installation was available to us. Using a remote cluster allows even more compute resources to be used during the execution of a job. Some for loops (where each iteration can be executed independently of the others) lend themselves to parallel execution using parfor, and we saw how we can run a batch job where our script contained a parfor loop.

There are more details about getting your code running on a cluster in the documentation. Do you have any other questions about running your code on a cluster? Let us know here.

Get the MATLAB code

Published with MATLAB® 7.14

Contents

New at MATLAB Central

In addition to the usual services on MATLAB Central,

there are three new additions to MATLAB Central collection. I'll list them here and delve into more details in later posts.

Which Community Service Do You Use Most?

I'm wondering which service you use most. Let me know here.

Get the MATLAB code

Published with MATLAB® 7.14

Contents

How objects spend their time

Let’s start with some basics – some of the places where objects spend time and how to minimize it. Objects will spend time in four basic places – object construction, property access, method invocation, and object deletion.

Object Construction

Object construction time is mostly spent copying the default values of properties from the class definition to the object and then calling the object’s constructor function(s). Generally speaking, there isn’t much to consider about default values in terms of performance since the expressions that create the default values are executed once when the class is first used, but then each object is just given a copy of the values from the class. Generally more superclasses will mean more constructor function calls for each object creation so this is a factor to consider in performance critical code.

Property Access

Property access is one of the most important factors in object performance. While changes over the past several releases have made property performance in R2012a more uniform and closer to struct performance, there is still some additional overhead for properties and the potential for aliasing with handle objects means that handle objects don’t get the same level of optimization that structs and value objects get from the MATLAB JIT. The simpler a property can be, the faster it will be. For example, making a property observable by listeners (using SetObservable/GetObservable) will turn off many optimizations and make property access slower. Using a set or get function will turn off most optimizations and also introduce the extra time to call the set or get function which is generally much greater than the time to just access the property. MATLAB doesn’t currently inline functions including set/get functions and so these are always executed as function calls. MATLAB optimizes property reading separately from property writing so it is important not to add a get-function just because the property needs a set-function.

Consider the following class:

type SimpleCylinder

classdef SimpleCylinder
    properties
        R
        Height 
    end
    
    methods
        function V = volume(C)
            V = pi .* [C.R].^2 .* [C.Height];
        end
    end
end

We can measure the time to create 1000 cylinders and compute their volumes:

tic
C1 = SimpleCylinder;
for k = 1:1000,
    C1(k).R = 1;
    C1(k).Height = k;
end
V = volume(C1);
toc

Elapsed time is 0.112309 seconds.

Now consider a slightly different version of the above class where the class checks all the property values:

type SlowCylinder

classdef SlowCylinder
    properties
        R
        Height
    end
    methods
        function V = volume(C)
            V = pi .* [C.R].^2 .* [C.Height];
        end
        
        function C = set.R(C, R)
            checkValue(R);
            C.R = R;
        end
        
        function C = set.Height(C, Height)
            checkValue(Height);
            C.Height = Height;
        end
    end
end

function checkValue(x)
    if ~isa(x, 'double') || ~isscalar(x)
        error('value must be a scalar double.');
    end
end

We can measure the same operations on this class:

tic
C2 = SlowCylinder;
for k = 1:1000,
    C2(k).R = 1;
    C2(k).Height = k;
end
A = volume(C2);
toc

Elapsed time is 0.174094 seconds.

Optimizing for Property Usage Inside the Class

If much of the performance critical code is inside methods of the class, it might make sense to consider using two property definitions for properties accessed in such performance critical code. One property is private to the class and doesn’t define any set or get functions. A second dependent property is public and passes through to the private property but adds error checking in its set-function. This allows the class to check values coming from outside the class, but not check values inside the class. Set functions always execute except when setting the default value during object creation and this allows the class to use its public interface if it is more convenient to do so. Set functions may do convenient transformations or other work in addition to just checking that the input value is legal for the property. However, if a performance-critical method doesn’t need this work, it can be helpful to use two properties.

For example, consider a new version of the cylinder class that checks its inputs but is designed to keep loops inside the class methods and use unchecked properties inside those methods.

type NewCylinder

classdef NewCylinder
    properties(Dependent)
        R
        Height
    end
    properties(Access=private)
        R_
        Height_
    end
    methods
        function C = NewCylinder(R, Height)
            if nargin > 0
                if ~isa(R, 'double') || ~isa(Height, 'double')
                    error('R and Height must be double.');
                end
                
                if ~isequal(size(R), size(Height))
                    error('Dimensions of R and Height must match.');
                end
                for k = numel(R):-1:1
                    C(k).R_ = R(k);
                    C(k).Height_ = Height(k);
                end
            end            
        end
        
        function V = volume(C)
            V = pi .* [C.R_].^2 .* [C.Height_];
        end
        
        function C = set.R(C, R)
            checkValue(R);
            C.R_ = R;
        end

        function R = get.R(C)
            R = C.R_;
        end
        
        function C = set.Height(C, Height)
            checkValue(Height);
            C.Height_ = Height;
        end
        
        function Height = get.Height(C)
            Height = C.Height_;
        end        
    end
end

function checkValue(x)
    if ~isa(x, 'double') || ~isscalar(x)
        error('value must be a scalar double.');
    end
end

Here we measure the same operations as above.

tic
C3 = NewCylinder(ones(1,1000), 1:1000);
A = volume(C3);
toc

Elapsed time is 0.006654 seconds.

Method Invocation

Method invocation using function call notation e.g. f(obj, data) is generally faster than using obj.f(data). Method invocation, like function calls on structs, cells, and function handles will not benefit from JIT optimization of the function call and can be many times slower than function calls on purely numeric arguments. Because of the overhead for calling a method, it is always better to have a loop inside of a method rather than outside of a method. Inside the method, if there is a loop, it will be faster if the loop just does indexing operations on the object and makes calls to functions that are passed numbers and strings from the object rather than method or function calls that take the whole object. If function calls on the object can be factored outside of loops, that will generally improve performance.

Calling a method on an object:

C4 = NewCylinder(10, 20);
tic
for k = 1:1000
    volume(C4);
end
toc

Elapsed time is 0.013509 seconds.

Calling a method on the object vector:

C5 = NewCylinder(ones(1,1000), 1:1000);
tic
volume(C5);
toc

Elapsed time is 0.001903 seconds.

Calling a function on a similar struct and struct array First calling the function inside a loop:

CS1 = struct('R', 10, 'Height', 20);
tic
for k = 1:1000
    cylinderVolume(CS1);
end

Next, we call the function on a struct array:

toc
CS2 = struct('R', num2cell(ones(1,1000)), ...
             'Height', num2cell(1:1000));
tic
cylinderVolume(CS2);
toc

Elapsed time is 0.008510 seconds.
Elapsed time is 0.000705 seconds.

Deleting Handle Objects

MATLAB automatically deletes handle objects when they are no longer in use. MATLAB doesn't use garbage collection to clean up objects periodically but instead destroys objects when they first become unreachable by any program. This means that MATLAB destructors (the delete method) are called more deterministically than in environments using garbage collection, but it also means that MATLAB has to do more work whenever a program potentially changes the reachability of a handle object. For example, when a variable that contains a handle goes out of scope, MATLAB has to determine whether or not that was the last reference to that variable. This is not as simple as checking a reference count since MATLAB has to account for cycles of objects. Changes in R2011b and R2012a have made this process much faster and more uniform. However, there is one aspect of object destruction that we are still working on and that has to do with recursive destruction. As of R2012a, if a MATLAB object is destroyed, any handle objects referenced by its properties will also be destroyed if no longer reachable and this can in turn lead to destroying objects in properties of those objects and so on. This can lead to very deep recursion for something like a very long linked list. Too much recursion can cause MATLAB to run out of system stack space and crash. To avoid such an issue, you can explicitly destroy elements in a list rather than letting MATLAB discover that the whole list can be destroyed.

Consider a doubly linked list of nodes using this node class:

type dlnode

classdef dlnode < handle
    properties
        Data
    end
    properties(SetAccess = private)
        Next
        Prev
    end
    
    methods
        function node = dlnode(Data)
            node.Data = Data;
        end

        function delete(node)
            disconnect(node);
        end

        function disconnect(node)
            prev = node.Prev;
            next = node.Next;
            if ~isempty(prev)
                prev.Next = next;
            end
            if ~isempty(next)
                next.Prev = prev;
            end
            node.Next = [];
            node.Prev = [];
        end
        
        function insertAfter(newNode, nodeBefore)
            disconnect(newNode);
            newNode.Next = nodeBefore.Next;
            newNode.Prev = nodeBefore;
            if ~isempty(nodeBefore.Next)
                nodeBefore.Next.Prev = newNode;
            end
            nodeBefore.Next = newNode;
        end
        
        function insertBefore(newNode, nodeAfter)
            disconnect(newNode);
            newNode.Next = nodeAfter;
            newNode.Prev = nodeAfter.Prev;
            if ~isempty(nodeAfter.Prev)
                nodeAfter.Prev.Next = newNode;
            end
            nodeAfter.Prev = newNode;
        end       
    end
end

    
    

Create a list of 1000 elements:

top = dlnode(0);
tic
for i = 1:1000
    insertBefore(dlnode(i), top);
    top = top.Prev;
end
toc

Elapsed time is 0.123879 seconds.

Destroy the list explicitly to avoid exhausting the system stack:

tic
while ~isempty(top)
    oldTop = top;
    top = top.Next;
    disconnect(oldTop);
end
toc

Elapsed time is 0.113519 seconds.

Measure time for varying lengths of lists. We expect to see time vary linearly with the number of nodes.

N = [500 2000 5000 10000];
% Create a list of 10000 elements:
CreateTime = [];
TearDownTime = [];
for n = N
    top = dlnode(0);
    tic
    for i = 1:n
        insertBefore(dlnode(i), top);
        top = top.Prev;
    end
    CreateTime = [CreateTime;toc];
    tic
    while ~isempty(top)
        oldTop = top;
        top = top.Next;
        disconnect(oldTop);
    end
    TearDownTime = [TearDownTime; toc];
end
subplot(2,1,1);
plot(N, CreateTime);
title('List Creation Time vs. List Length');
subplot(2,1,2);
plot(N, TearDownTime);
title('List Destruction Time vs. List Length');

A Look to the Future

We continue to look for opportunities to improve MATLAB object performance and examples from you are very helpful for learning what changes will make an impact on real applications. If you have examples or scenarios you want us to look at, please let me know. Also, if you have your own ideas or best practices, it would be great to share them as well. You can post ideas and comments here.

Get the MATLAB code

Published with MATLAB® 7.14

Contents

The 12a release of Statistics Toolbox includes new functions for

These regression techniques aren’t new to Statistics Toolbox. What is new is that MathWorks addded a wide set of support functions that simplify common analysis tasks like plotting, outlier detection, generating predictions, performing stepwise regression, applying robust regression...

We'll start with a simple example using linear regression.

Create a dataset

I'm going to generate a basic dataset in which the relationship between X and Y is modeled by a straight line (Y = mX + B) and add in some normally distributed noise. Next, I'll generate a scatter plot showing the relationship between X and Y

clear all
clc
rng(1998);

X = linspace(1,100,50);
X = X';
Y = 7*X + 50 + 30*randn(50,1);
New_X = 100 * rand(10,1);

scatter(X,Y, '.')

Use linear regression to model Y as a function of X

Looking at the data, we can see a clear linear relationship between X and Y. I'm going to use the new LinearModel function to model Y as a function of X. The output from this function will be stored as an object named "myFit" which is displayed on the screen.

myFit = LinearModel.fit(X,Y)

myFit = 
Linear regression model:
    y ~ 1 + x1

Estimated Coefficients:
                   Estimate    SE         tStat     pValue    
    (Intercept)      51.6       8.2404    6.2618    9.9714e-08
    x1             7.0154      0.14131    49.644    6.4611e-43

Number of observations: 50, Error degrees of freedom: 48
Root Mean Squared Error: 29.1
R-squared: 0.981,  Adjusted R-Squared 0.98
F-statistic vs. constant model: 2.46e+03, p-value = 6.46e-43

The first line shows the linear regression model. When perform a regression we need to specify a model that describes the relationship between our variables. By default, LinearModel assumes that you want to model the relationship as a straight line with an intercept term. The expression "y ~ 1 + x1" describes this model. Formally, this expression translates as "Y is modeled as a linear function which includes an intercept and a variable". Once again note that we are representing a model of the form Y = mX + B...

The next block of text includes estimates for the coefficients, along with basic information regarding the reliability of those estimates.

Finally, we have basic information about the goodness-of-fit including the R-square, the adjusted R-square and the Root Mean Squared Error.

Use myFit for analysis

Earlier, I mentioned that the new regression functions include a wide variety of support functions that automate different analysis tasks. Let's look at some them.

First, lets generate a plot that we can use to evaluate the quality of the resulting fit. We'll do so by applying the standard MATLAB plot command to "myFit".

plot(myFit)

Notice that this simple command creates a plot with a wealth of information including

MATLAB has also automatically labelled our axes and added a legend.

Look for patterns in the residuals

Alternatively, lets assume that we wanted to see whether there was any pattern to the residuals. (A noticeable pattern to the residuals might suggest that our model is to simple and that it failed to capture a real work trend in the data set. This technique can also be used to check and see whether the noise component is constant across the dataset).

Here, I'll pass "myFit" to the new "plotResiduals" method and tell plotResiduals to plot the residuals versus the fitted values.

figure
plotResiduals(myFit, 'fitted')

My plot looks like random noise - which in this case is a very good thing.

Look for autocorrelation in the residuals

Autocorrelation in my data set could also throw off the quality of my fit. The following command will modify the residual plot by plotting residuals versus lagged residuals.

figure
plotResiduals(myFit, 'lagged')

Look for outliers

Suppose that I wanted to check for outliers... We also have a plot for that.

figure
plotDiagnostics(myFit, 'cookd')

Cook's Distance is a metric that is commonly used to see whether a dataset contains any outliers. For any given data point, Cook's Distance is calculated by performing a brand new regression that excludes that data point. Cook's distance measures how much the shape of the curve changes between the two fits. If the curve moves by a large amount, that data point has a great deal of influence on the model and might very well be an outlier.

In this example, non of our data points look as if they are outliers.

Use the resulting model for prediction

Last, but not least, lets assume that we wanted to use our model for prediction. This is as easy as applying the "predict" method.

Predictions = predict(myFit, New_X)

Predictions =
       310.51
       596.17
       64.866
       490.79
       306.06
       308.87
       621.55
       366.15
       543.32
       317.33

Discover the full set of methods available with regression objects

I hope that you agree that all these built in plots and analysis routines represent a significant improvement in usability. However, if you're anything like me, you immediate reaction is going to be "Great, you've built a lot of nice stuff, however, how do you expect me to find out about this?"

What I'd like to do now is show you a couple of simple tricks that you can use to discover all the new cabilities that we've added. The first trick is to recognize that "myFit" is an object and that objects have methods associated with them. All of the commands that we've sued so far like "plot", "plotResiduals", and "predict" are methods for the LinearModel object.

Any time that I'm working with one of the built in objects that ship with MATLAB my first plot is to inspect the full set of methods that ship with that object. This is as easy as typing methods(myFit) at the command line. I can use this to immediately discover all the built in capabilities that ship with the object. If one of those options catches my eye, I can use the help system to get more information.

methods(myFit)

Methods for class LinearModel:

addTerms              plot                  plotSlice             
anova                 plotAdded             predict               
coefCI                plotAdjustedResponse  random                
coefTest              plotDiagnostics       removeTerms           
disp                  plotEffects           step                  
dwtest                plotInteraction       
feval                 plotResiduals         

Static methods:

fit                   stepwise              

Discover the full set of information included in the regression object

Here's another really useful trick to learn about the new regression objects. You can use the MATLAB variable editor to walk through the object and see all the information that is availabe.

You should have an object named "myFit" in the MATLAB workspace. Double clicking on the object will open the object in the Variable Editor.

Formulas

At the start of this blog there was some brief introduction to "formulas". I'd like to conclude this talk by providing a bit more information about formulas. Regression analysis requires the ability to specify a model that describes the relationship between your predictors and your response variables.

Let's change our initial example such that we're working with a high order polynomial rather than a straight line. I'm also going to change this from a curve fitting problem to a surface fitting problem.

X1 = 100 * randn(100,1);
X2 = 100 * rand(100,1);
X = [X1, X2];
Y = 3*X1.^2 + 5*X1.*X2 + 7* X2.^2 + 9*X1 + 11*X2 + 30 + 100*randn(100,1);
myFit2 = LinearModel.fit(X,Y)

myFit2 = 
Linear regression model:
    y ~ 1 + x1 + x2

Estimated Coefficients:
                   Estimate    SE       tStat     pValue   
    (Intercept)     25771      11230    2.2949     0.023893
    x1             178.75      54.95     3.253    0.0015717
    x2             613.39      187.2    3.2767    0.0014579

Number of observations: 100, Error degrees of freedom: 97
Root Mean Squared Error: 5.42e+04
R-squared: 0.211,  Adjusted R-Squared 0.195
F-statistic vs. constant model: 13, p-value = 1.03e-05

Let's take a look at the output from this example. We can see, almost immediately, that something has gone wrong with our fit.

If we look at the line that describes the linear regression model we can see what went wrong. By default, LinearModel will fit a plane to the dataset. Here, the relationship between X and Y is modelled as a high order polynomial. We need to pass this additional piece of information to "LinearModel".

Modeling a high order polynomial (Option 1)

Here are a couple different ways that I can use LinearModel to model a high order polynomial. The first option is to write out the formula by hand.

myFit2 = LinearModel.fit(X,Y, 'y ~ 1 + x1^2 + x2^2 + x1:x2 + x1 + x2')

myFit2 = 
Linear regression model:
    y ~ 1 + x1*x2 + x1^2 + x2^2

Estimated Coefficients:
                   Estimate    SE           tStat      pValue     
    (Intercept)    32.317         38.071    0.84886        0.39811
    x1             9.1735         0.2119     43.291     6.9346e-64
    x2              12.19         1.6568     7.3577     6.9388e-11
    x1:x2          4.9985      0.0033943     1472.6    7.0691e-207
    x1^2           2.9992      0.0006268       4785    5.5265e-255
    x2^2           6.9849       0.015281     457.11    3.9814e-159

Number of observations: 100, Error degrees of freedom: 94
Root Mean Squared Error: 101
R-squared: 1,  Adjusted R-Squared 1
F-statistic vs. constant model: 7.02e+06, p-value = 3.23e-260

Modeling a high order polynomial (Option 2)

Alternatively, I can simple use the the string "poly22" to indicate a second order polynomial for both X1 and X2 an automatically generate all the appropriate terms and cross terms.

myFit2 = LinearModel.fit(X, Y, 'poly22')

myFit2 = 
Linear regression model:
    y ~ 1 + x1^2 + x1*x2 + x2^2

Estimated Coefficients:
                   Estimate    SE           tStat      pValue     
    (Intercept)    32.317         38.071    0.84886        0.39811
    x1             9.1735         0.2119     43.291     6.9346e-64
    x2              12.19         1.6568     7.3577     6.9388e-11
    x1^2           2.9992      0.0006268       4785    5.5265e-255
    x1:x2          4.9985      0.0033943     1472.6    7.0691e-207
    x2^2           6.9849       0.015281     457.11    3.9814e-159

Number of observations: 100, Error degrees of freedom: 94
Root Mean Squared Error: 101
R-squared: 1,  Adjusted R-Squared 1
F-statistic vs. constant model: 7.02e+06, p-value = 3.23e-260

Nonlinear regression: Generate our data

Let's consider a nonlinear regression example. This time around, we'll work with a sin curve. The equation for a sin curve is governed by four key parameters.

We'll start by generating a our dataset

X = linspace(0, 6*pi, 90);
X = X';
Y = 10 + 3*(sin(1*X + 5)) + .2*randn(90,1);

Nonlinear regression: Generate a fit

Next we'll using the NonLinearModel function to perform a nonlinear regression. Here we need to

myFit3 = NonLinearModel.fit(X,Y, 'y ~ b0 + b1*sin(b2*x + b3)', [11, 2.5, 1.1, 5.5])

myFit3 = 
Nonlinear regression model:
    y ~ b0 + b1*sin(b2*x + b3)

Estimated Coefficients:
          Estimate    SE           tStat     pValue     
    b0     9.9751       0.02469    404.02    9.0543e-143
    b1     2.9629      0.033784    87.703     6.4929e-86
    b2    0.99787     0.0021908    455.47     3.029e-147
    b3     5.0141      0.023119    216.88    1.4749e-119

Number of observations: 90, Error degrees of freedom: 86
Root Mean Squared Error: 0.227
R-Squared: 0.989,  Adjusted R-Squared 0.989
F-statistic vs. constant model: 2.58e+03, p-value = 4.36e-84

Nonlinear regression: Work with the resulting model

Once again, the output from the regression analysis is an object which we can used for analysis. For example:

figure
scatter(X,Y)
hold on
plot(X, myFit3.Fitted, 'r')

Nonlinear regression: Alternative ways to specify the regression model

One last point to be aware of: The syntax that I have been using to specify regression models is based on "Wilkinson's notation". This is a standard syntax that is commonly used in Statistics. If you prefer, you also have the option to specify your model using anonymous functions.

For example, that command could have been written as

myFit4 = NonLinearModel.fit(X,Y, @(b,x)(b(1) + b(2)*sin(b(3)*x + b(4))), [11, 2.5, 1.1, 5.5])

myFit4 = 
Nonlinear regression model:
    y ~ (b1 + b2*sin(b3*x + b4))

Estimated Coefficients:
          Estimate    SE           tStat     pValue     
    b1     9.9751       0.02469    404.02    9.0543e-143
    b2     2.9629      0.033784    87.703     6.4929e-86
    b3    0.99787     0.0021908    455.47     3.029e-147
    b4     5.0141      0.023119    216.88    1.4749e-119

Number of observations: 90, Error degrees of freedom: 86
Root Mean Squared Error: 0.227
R-Squared: 0.989,  Adjusted R-Squared 0.989
F-statistic vs. constant model: 2.58e+03, p-value = 4.36e-84

Conclusion

I've been working with Statistics Toolbox for close to five years now. Other than the introduction of dataset arrays a few years back, nothing has gotten me nearly as excited as the release of these new regression capabilities. So, it seemed fitting to conclude with an example that combines dataset arrays and the new regression objects.

Tell me what you think is going on in the following example. (Extra credit if you can work in the expression "dummy variable")

load carsmall

ds = dataset(MPG,Weight);
ds.Year = ordinal(Model_Year);

mdl = LinearModel.fit(ds,'MPG ~ Year + Weight^2')

mdl = 
Linear regression model:
    MPG ~ 1 + Weight + Year + Weight^2

Estimated Coefficients:
                   Estimate      SE            tStat      pValue    
    (Intercept)        54.206        4.7117     11.505    2.6648e-19
    Weight          -0.016404     0.0031249    -5.2493    1.0283e-06
    Year_76            2.0887       0.71491     2.9215     0.0044137
    Year_82            8.1864       0.81531     10.041    2.6364e-16
    Weight^2       1.5573e-06    4.9454e-07      3.149     0.0022303

Number of observations: 94, Error degrees of freedom: 89
Root Mean Squared Error: 2.78
R-squared: 0.885,  Adjusted R-Squared 0.88
F-statistic vs. constant model: 172, p-value = 5.52e-41

Post your answers here.

Get the MATLAB code

Published with MATLAB® 7.14

Contents

Conference Information

You'll find information about the conference here as well as a link for registration. You can also see me extend a video welcome to each of you.

Brief Conference Agenda

In addition to a keynote address from Jim Tung, and the following tracks:

The conference will have venues to meet MathWorks folks and exhibitors.

Join Us!

On behalf of Jim, MathWorks, and myself, we look forward to seeing you at the MATLAB Virtual Conference!

Get the MATLAB code

Published with MATLAB® 7.14

Contents

Setting up the Arduino and Breadboard

Depending on the platform you are working on, Arduino's website has a detailed set of instructions for setting up the Arduino IDE and connecting the Arduino board. Once the setup was complete, I worked on the breadboard to have the setup similar to CIRC-08 from the Sparkfun's Arduino Inventor Kit, except that I used an RGB LED instead of Green LED and hence instead of using just pin 9 - I used pins 9, 10 and 11.

At this stage the Arduino board and rest of the components were ready.

Setting up MATLAB and Support Package

MATLAB support package for Arduino comes with a server program adiosrv.pde that can be downloaded on the Arduino board. Once downloaded, the adiosrv.pde will:

The steps I performed to write/download the server program to Arduino board were as following:

I was now ready to code in MATLAB.

Writing MATLAB Code

MATLAB support package for Arduino provides a very easy to use function set to perform basic operations like analog read, analog write, digital read and digital write.

For the purpose of this example, I followed the following sequence:

In MATLAB the above sequence looks like:

 % Connect to the board
 a = arduino('COM8');

 a.pinMode(9,'output');
 a.pinMode(10,'output');
 a.pinMode(11,'output');

 % Start the timer now
 tic;

 while toc/60 < 1 % Run for 1 minute.
 % Read analog input from analog pin 0
 % The value returned is between 0 and 1023
     sensorValue = a.analogRead(0);

 % For potentiometer value from 0 to 511, we will fade LED from red to
 % green keeping blue at constant 0.
    if 0 <= sensorValue && sensorValue < 512
       greenIntensity = floor(sensorValue/2.0);
       redIntensity   = 255-greenIntensity;
       blueIntensity  = 0;
    else
 % For potentiometer value from 511 to 1023, we will fade LED from
 % green to blue keeping red at constant 0.
       blueIntensity   = floor(sensorValue/2.0) - 256;
       greenIntensity  = 255 - blueIntensity;
       redIntensity    = 0;
    end

 % Write the intensity to analog output pins
    a.analogWrite(9,redIntensity);
    a.analogWrite(10,greenIntensity);
    a.analogWrite(11,blueIntensity);
 end

 % Close session
 delete(a);
 clear a;

Conclusion

So why, you might ask, do I need MATLAB when I can just code it in Arduino IDE? The simple answer is - interactive development and debugging capabilities.

While working on this project, there were so many instances where I wanted to see the values returned from the potentiometer as I was turning the knob or quickly set the RGB value of the LED. Using the Arduino IDE, I would have updated the code, re-compiled it and uploaded the code to the board before I can see the updates in effect. With MATLAB support package, I was able to just tweak the settings on the fly without any extra step.

Make sure you check out the webinar: Learning Basic Mechatronics Concepts Using the Arduino Board and MATLAB, to learn more about Analog and Digital I/O as well as DC, Servo and Stepper motor control.

Have you used/wanted to use Arduino for any cool projects? If so, please feel free to share your project details in here.

Get the MATLAB code

Published with MATLAB® 7.13

Contents

Errors are Different!

By default, standalone applications display both ordinary output and and errors via the same channel. This can be convenient, as it makes capturing all the output relatively simple. With this simplicity comes a hidden danger: error messages may get lost if the application prints out a lot of messages. In this final installment of my series of posts about redirecting program output, I'm going to show you how to treat errors differently than other types of messages. One reason for doing so is to make error messages more obvious, so they are easier to find. When I'm looking through a long listing of program output, I'm often scanning primarily for error messages, because errors indicate some intervention is necessary. I'd like the errors to be emphasized or decorated so I'll spot them quickly.

Applications and shared libraries created by MATLAB Compiler split the output from MATLAB functions, into two streams, one for errors and one for all other types of output, including warnings. Since MATLAB Compiler-generated components allow you to install a separate handler for each of these streams, it is conceptually simple to emphasize error messages.

A Brief Review of the Application

monitor is a C program that calls a single MATLAB function, sendmessages. The MATLAB function displays three kinds of output: some informational text, a warning message and an error message. monitor.c relies on three monitoring functions, StartMonitoring, MonitorHandler and StopMonitoring, to redirect its output to a user-visible window. Separate, platform-specific files contain different sets of monitoring functions.

Download the source code for the monitor program from MATLAB Central. See the monitorREADME in the ZIP file for a description of each of the downloaded files and complete directions for building the application.

Splitting the Error and Output Streams

The generated library initialization function libmsgfcnInitializeWithHandlers takes two input arguments: an error message handler and a print message handler. The example monitor.c I discussed in last week's post installs the same handler, the function MonitorHandler, for both the error and the print output streams:

 if (!libmsgfcnInitializeWithHandlers(MonitorHandler, MonitorHandler))

Change the first parameter to a specialized error handler, the ErrorHandler function, and the monitor program decorates the error messages it displays, making them easier to see.

 if (!libmsgfcnInitializeWithHandlers(ErrorHandler, MonitorHandler))

Emphasizing Error Messages

I've chosen to emphasize error messages in two different ways. On UNIX the monitor window displays error messages in red. The UNIX ErrorHandler function uses X terminal escape codes to set the color of the displayed text. ErrorHandler sends the sequence ESC [ 1;31m to set the text to red, outputs the error message, and then restores the foreground color to black with the sequence ESC [ 1;30m.

int ErrorHandler(const char *s)
{
    const int csi = 0x9b;
    const char *redOn = "1;31m";
    const char *redOff = "1;30m";

    time_t now = time(NULL);
    char nowIsTheHour[128];
    fwrite(&csi, sizeof(int), 1, tmpLog);
    fprintf(tmpLog, "%s", redOn);
    sprintf(nowIsTheHour, "%s", ctime(&now));
    fprintf(tmpLog, nowIsTheHour);
    fprintf(tmpLog, "%s\n", s);
    fwrite(&csi, sizeof(int), 1, tmpLog);
    fprintf(tmpLog, "%s", redOff);
    return strlen(s)+strlen(nowIsTheHour)+1;
}

The UNIX monitoring mechanism sends its output to a temporary file, the contents of which are displayed in an X terminal running tail -f. Refer to the previous post in this series for complete details of the implementation.

The monitoring mechanism in the Windows program monitorWinDialog uses a dialog box containing a ListBox. The DialogProc function (the function that handles events directed at the dialog box), prefixes a line of five asterixes to error messages before displaying them:

case ERROR_MESSAGE:
    strcpy(buffer, "***** ");
    strcat(buffer, (const char *)lParam);
    ListBox_AddString(GetDlgItem(box, IDC_LIST1),
                      (const char *)buffer);
    break;

The error handler function ErrorHandler in the Microsoft Windows event log version, monitorWinEvents calls ReportEvent with an event type of EVENTLOG_ERROR_TYPE.

if (!ReportEvent(hEventLog, EVENTLOG_ERROR_TYPE, EXAMPLE_CATEGORY,
                 INFORMATIONAL_MESSAGE, NULL, 1, 0,
                 (LPCSTR*)pInsertStrings, NULL))

In response, the event log viewer displays a stop-sign icon next to the event. The print handler MonitorHandler creates events with an event type of EVENTLOG_INFORMATION_TYPE, causing the event logger to display a "talk bubble" icon next to the event.

Seeing the Differences

Once you've changed monitor.c to install separate handlers on the error and print message streams, build and run the application.

First use MATLAB Compiler to create the shared library containing sendmessages:

>> mcc -W lib:libmsgfcn -T link:lib sendmessages.m

Then, use mbuild to build the main program and link it with the shared library:

(UNIX): mbuild -g monitor.c monitorUNIX.c -L. -lmsgfcn

(Windows) -- two commands because there are two programs:

mbuild -g monitorWinDialog.c monitor.c WinLogger.rc libmsgfcn.lib

mbuild -g monitorWinEvents.c monitor.c libmsgfcn.lib

The first command builds the dialog box version (monitorWinDialog.exe) and the second command builds the event log version (monitorWinEvents.exe).

Now run it. The program takes three input arguments: the informational message, the warning and the error message.

(UNIX): ./monitor information "look out!" oops

(Windows): monitorWinDialog.exe information "look out!" oops

It should be easy to pick out the error message in the program output.

That's all I've got to say about managing errors and output in a deployed application. Do your applications do things differently? Are any of these ideas going to be useful to you? What would you like to see the Compiler do differently? Let me know.

Get the MATLAB code

Published with MATLAB® 7.13

Contents

Overview

In this post, we first will introduce the basics of using the GPU with MATLAB and then move onto solving a 2nd-order wave equation using this GPU functionality. This blog post is inspired by a recent MATLAB Digest article on GPU Computing that I coauthored with one of our developers, Jill Reese. Since the original demo was made, the GPU functions available in MATLAB have grown. If you compare the code below to the code in the paper- they are slightly different, reflecting these new capabilites. Additionally, small changes were made to enable easier explanation of the code in this blog format. The GPU functionality shown in this post requires the Parallel Computing Toolbox.

GPU Background

Originally used to accelerate graphics rendering, GPUs are now increasingly applied to scientific calculations. Unlike a traditional CPU, which includes no more than a handful of cores, a GPU has a massively parallel array of integer and floating-point processors, as well as dedicated, high-speed memory. A typical GPU comprises hundreds of these smaller processors. These processors can be used to greatly speed-up particular types of applications.

A good rule of thumb is that your problem may be a good fit for the GPU if it is:

Applications that do not satisfy these criteria might actually run more slowly on a GPU than on a CPU.

Learning About Our GPU

With that background, we can now start working with the GPU in MATLAB. Let's start first by querying our GPU to see just what we are working with:

gpuDevice

ans = 
  parallel.gpu.CUDADevice handle
  Package: parallel.gpu

  Properties:
                      Name: 'Tesla C2050 / C2070'
                     Index: 1
         ComputeCapability: '2.0'
            SupportsDouble: 1
             DriverVersion: 4
        MaxThreadsPerBlock: 1024
          MaxShmemPerBlock: 49152
        MaxThreadBlockSize: [1024 1024 64]
               MaxGridSize: [65535 65535]
                 SIMDWidth: 32
               TotalMemory: 3.1820e+009
                FreeMemory: 2.6005e+009
       MultiprocessorCount: 14
              ClockRateKHz: 1147000
               ComputeMode: 'Default'
      GPUOverlapsTransfers: 1
    KernelExecutionTimeout: 1
          CanMapHostMemory: 1
           DeviceSupported: 1
            DeviceSelected: 1

We are running on a Tesla C2050. Currently, Parallel Computing Toolbox supports NVDIA GPUs with Compute Capability 1.3 or greater. Here is a link to see if your card is supported.

A Simple Example using Overloaded Functions

Over 100 operations (e.g. fft, ifft, eig) are now available as built-in MATLAB functions that can be executed directly on the GPU by providing an input argument of the type GPUArray. These GPU-enabled functions are overloaded—in other words, they operate differently depending on the data type of the arguments passed to them. Here is a list of all the overloaded functions.

Let's create a GPUArray and perform a fft using the GPU. However, let's first do this on the CPU so that we can see the difference in code and performance

A1 = rand(3000,3000);
tic;
B1 = fft(A1);
time1 = toc;

To perform the same operation on the GPU, we first use gpuArray to transfer data from the MATLAB workspace to device memory. Then we can run fft, which is one of the overloaded functions on that data:

A2 = gpuArray(A1);
tic;
B2 = fft(A2);
time2 = toc;

The second case is executed on the GPU rather than the CPU since its input (a GPUArray) is held on the GPU. The result, B2, is stored on the GPU. However, it is still visible in the MATLAB workspace. By running class(B2), we can see that it is also a GPUArray.

class(B2)

ans =
parallel.gpu.GPUArray

To bring the data back to the CPU, we run gather.

B2 = gather(B2);
class(B2)

ans =
double

B2 is now on the CPU and has a class of double.

Speedup For Our Simple Example

In this simple example, we can calculate the speedup of performing the fft on our data.

speedUp = time1/time2;
disp(speedUp)

    3.6122

It looks like our fft is running ~3.5x faster on the GPU. This is pretty good, especially since fft is multi-threaded in core MATLAB. However, to perform a realistic comparison, we should really include the time spent transferring the vector to and from the GPU. If we do that - we find that our acceleration is greatly reduced. Let's see what happens if we include the time to do our data transfer.

tic;
A3 = gpuArray(A1);
B3 = fft(A3);
B3 = gather(B3);
time3 = toc;

speedUp = time1/time3;
disp(speedUp)

    0.5676

Understanding and Limiting Your Data Transfer Overhead

Data transfer overhead can become so significant that it degrades the application's overall performance, especially if you repeatedly exchange data between the CPU and GPU to execute relatively few computationally intensive operations. However not all hope is lost! By limiting the data transfer between the GPU and the CPU, we can still acheive speedup.

Instead of creating the data on the CPU and transferring it to the GPU - we can just create on the GPU directly. There are several constructors available as well as constructor-like functions such as meshgrid. Let's see how that effects our time. To be accurate, we should also retime the original serial code including the time it takes to generate the random matrix on the CPU.

tic;
A4 = rand(3000,3000);
B4 = fft(A4);
time4 = toc;

tic;
A5 = parallel.gpu.GPUArray.rand(3000,3000);
B5 = fft(A5);
B5 = gather(B5);
time5 = toc;

speedUp = time4/time5;
disp(speedUp);

    1.4100

This is better, although we still do see the influence of gathering the data back from the GPU. In this simple demo, the effect is exaggerated because we are running a single, already fast operation on the GPU. It is more efficient to perform several operations on the data while it is on the GPU, bringing the data back to the CPU only when required.

Solving the Wave Equation

To put the above example into context, let's use the same GPU functionality on a more "real-world" problem. To do this, we are going to solve a second-order wave equation:

The algorithm we use to solve the wave equation combines a spectral method in space and a second-order central finite difference method in time. Specifically, we apply the Chebyshev spectral method, which uses Chebyshev polynomials as the basis functions. Our example is largely based on an example in Trefethen's book: "Spectral Methods in MATLAB"

Code Changes to Run Algorithm on GPU

When accelerating our alogrithm, we focus on speeding up the code within the main time stepping while-loop. The operations in that part of the code (e.g. fft and ifft, matrix multiplication) are all overloaded functions that work with the GPU. As a result, we do not need to change the algorithm in any way to execute it on a GPU. We get to simply transfer the data to the GPU and back to the CPU when finished.

type('stepSolution.m')

function [vv,vvold] = stepsolution(vv,vvold,ii,N,dt,W1T,W2T,W3T,W4T,...
                                        WuxxT1,WuxxT2,WuyyT1,WuyyT2)
    V = [vv(ii,:) vv(ii,N:-1:2)];
    U = real(fft(V.')).';
    
    W1test = (U.*W1T).';
    W2test = (U.*W2T).';
    W1 = (real(ifft(W1test))).';
    W2 = (real(ifft(W2test))).';
    
    % Calculating 2nd derivative in x
    uxx(ii,ii) = W2(:,ii).* WuxxT1 - W1(:,ii).*WuxxT2;
    uxx([1,N+1],[1,N+1]) = 0;
    
    V = [vv(:,ii); vv((N:-1:2),ii)];
    U = real(fft(V));
    
    W1 = real(ifft(U.*W3T));
    W2 = real(ifft(U.*W4T));
    
    % Calculating 2nd derivative in y
    uyy(ii,ii) = W2(ii,:).* WuyyT1 - W1(ii,:).*WuyyT2;
    uyy([1,N+1],[1,N+1]) = 0;
    
    % Computing new value using 2nd order central finite difference in
    % time
    vvnew = 2*vv - vvold + dt*dt*(uxx+uyy);
    vvold = vv; vv = vvnew;

end

Code Changes to Transfer Data to the GPU and Back Again

The variables dt, x, index1, and index2 are calculated on the CPU and transferred over to the GPU using gpuArray.

For example:

N = 256;
index1 = 1i*[0:N-1 0 1-N:-1];

index1 = gpuArray(index1);

For the rest of the variables, we create them directly on the GPU using the overloaded versions of the transpose operator ('), repmat ,and meshgrid.

For example:

W1T = repmat(index1,N-1,1);

When, we are done with all our calculations on the GPU, we use gather once to bring data back from the GPU. Note, that we don't have to transfer our data back to the CPU between time steps. This all comes together to look like this:

type WaveGPU.m

function vvg = WaveGPU(N,Nsteps)
%% Solving 2nd Order Wave Equation Using Spectral Methods
% This example solves a 2nd order wave equation: utt = uxx + uyy, with u =
% 0 on the boundaries. It uses a 2nd order central finite difference in
% time and a Chebysehv spectral method in space (using FFT). 
%
% The code has been modified from an example in Spectral Methods in MATLAB
% by Trefethen, Lloyd N.

% Points in X and Y
x = cos(pi*(0:N)/N); % using Chebyshev points

% Send x to the GPU
x = gpuArray(x);
y = x';

% Calculating time step
dt = 6/N^2;

% Setting up grid
[xx,yy] = meshgrid(x,y);

% Calculate initial values
vv = exp(-40*((xx-.4).^2 + yy.^2));
vvold = vv;

ii = 2:N;
index1 = 1i*[0:N-1 0 1-N:-1];
index2 = -[0:N 1-N:-1].^2;

% Sending data to the GPU
dt = gpuArray(dt);
index1 = gpuArray(index1);
index2 = gpuArray(index2);

% Creating weights used for spectral differentiation 
W1T = repmat(index1,N-1,1);
W2T = repmat(index2,N-1,1);
W3T = repmat(index1.',1,N-1);
W4T = repmat(index2.',1,N-1);

WuxxT1 = repmat((1./(1-x(ii).^2)),N-1,1);
WuxxT2 = repmat(x(ii)./(1-x(ii).^2).^(3/2),N-1,1);
WuyyT1 = repmat(1./(1-y(ii).^2),1,N-1);
WuyyT2 = repmat(y(ii)./(1-y(ii).^2).^(3/2),1,N-1);

% Start time-stepping
n = 0;

while n < Nsteps
    [vv,vvold] = stepsolution(vv,vvold,ii,N,dt,W1T,W2T,W3T,W4T,...
                                 WuxxT1,WuxxT2,WuyyT1,WuyyT2);
    n = n + 1;
end


% Gather vvg back from GPU memory when done
vvg = gather(vv);

Comparing CPU and GPU Execution Speeds

We ran a benchmark in which we measured the amount of time it took to execute 50 time steps for grid sizes of 64, 128, 512, 1024, and 2048 on an Intel Xeon Processor X5650 and then using an NVIDIA Tesla C2050 GPU.

For a grid size of 2048, the algorithm shows a 7.5x decrease in compute time from more than a minute on the CPU to less than 10 seconds on the GPU. The log scale plot shows that the CPU is actually faster for small grid sizes. As the technology evolves and matures, however, GPU solutions are increasingly able to handle smaller problems, a trend that we expect to continue.

Other Ways to Use the GPU with MATLAB

To accelerate an algorithm with multiple element-wise operations on a GPU, you can use arrayfun, which applies a function to each element of an GPUarray. Additionally if you have your own CUDA code, you can use the CUDAKernel interface to integrate this code with MATLAB.

Ben Tordoff's previous guest blog post Mandelbrots Sets on the GPU shows a great example using both of these capabilities.

Addendum: There has been some confusion about using assignin and evalin inside the body of parfor loops. You can use these functions if they use the 'base' workspace. However, please use caution when using them since the 'base' workspace will not be the same as the client 'base' workspace.

Your Thoughts?

Have you experimented with our new GPU functionality? Let us know what you think or if you have any questions by leaving a comment for this post here.

Get the MATLAB code

Published with MATLAB® 7.13