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

Contents

Find Documentation On-Line

If you go to MathWorks website, you can search in the box in the upper right corner for the term "documention". Here's a snippet of the screen you see next.

Find New Beta Documentation Center

Choose the second link, 'Online Documentation for MathWorks Product Family', and you reach a page that looks like this one.

See the prominent green box near the top 'Try the Documentation Center Beta' to see the new way to search for topics in our documentation and you will see a page that looks like this one.

Documentation Center Features

In Release 2011b, the documentation already has many features that allow you to pinpoint relevant information quickly.

Browse Documentation Topics by category

Search Documentation Topics

* Get search suggestions as you type, organized by content type * Refine search results by product, category, and content type * Single-click access to Functions and Blocks in a product * Expand search results to all Support content

Learning About Curve Fitting

Let's see what happens when we search for 'curve fitting' in the Doc Center. Here you see a screen with several different areas of focus.

You can see relevant documentation in MathWorks products in the main, right-hand portion of the window. In addition, the left side contains information clustered in other ways, allowing you to refine the search overall, by the type of information (functions, blocks, etc.), or by product.

Documentation Center

We hope you like the presentation of information in the Documentation Center and find it helpful. We have plans for more features for the Documentation Center. Please feel free to post your thoughts here.

Get the MATLAB code

Published with MATLAB® 7.13

Contents

My List of Best Practices

Clearly (at least to me), this is not everything you generally need to do. You still need to comment the code, add good help information and examples, etc. But these are the main coding practices and tools I always rely on.

Missing from Your List? Additions to My List?

What's on my list that you don't currently do? Do you have a major addition to my list (there can't be too many, or I won't remember to do them all!)? Let me know here.

Get the MATLAB code

Published with MATLAB® 7.13

Contents

Integrating Output with a Window System

Applications developed with the MATLAB Compiler direct their output to the console window. While this may suffice for batch or background tasks, it is often not suitable for interactive or graphical environments. In an earlier post, I demonstrated how to use message handlers to redirect program output and error messages to a log file. For many batch or background applications, simply capturing output to a log file provides sufficient information for validation and debugging. However, if you need to more closely monitor the progress of your application while it is running, you might want to redirect output to a visible window.

I've written a short C program in monitor.c to demonstrate a few applications of this idea: on UNIX, the program sends output to an X terminal; on Windows, I created two programs. The simplest of the Windows programs displays the output in a Windows dialog box. The other, slightly more complex, sends the output as events to the Windows Event Log.

In all three cases (UNIX and Windows), the monitor main program calls a MATLAB function, sendmessages, contained in a shared library generated by MATLAB Compiler.

Creating a Shared Library

The sendmessages function displays an informational message, a warning and then issues an error:

function sendmessages(info, warn, err)
    disp(info); warning(warn); error(err);

Create a C shared library containing sendmessages with MATLAB Compiler:

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

To follow the rest of this posting, download the source code for the main program that calls sendmessages 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 example application.

Build the application using mbuild:

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).

You can't run monitorWinEvents.exe on Microsoft Windows until you register the event source. (UNIX readers can skip this step.) To register the event source, follow the step-by-step instructions in the monitorREADME file in the source code distribution ZIP file.

Displaying Messages in a Visible Window

The main program in monitor.c relies on three external functions:

void StartMonitoring(void) void StopMonitoring(void) int
MonitorHandler(const char *msg)

StartMonitoring initializes the monitoring system, MonitorHandler redirects the program's output to the monitoring window and StopMonitoring releases the resources required to capture the program's output. The implementation of these functions varies by platform.

On UNIX, for example, the StartMonitor function starts an xterm which runs tail -f to monitor the log file:

void StartMonitoring() {
    if (tmpLog == NULL) {
        char logFileName[80]; char cmd[128];

        sprintf(logFileName, "/tmp/monitorLog.%d", getpid()); tmpLog =
        fopen(logFileName, "w"); sprintf(cmd, "xterm -e tail -f %s&",
        logFileName); system(cmd);
    }
}

The UNIX message handler, MonitorHandler, is identical to the message handler from the previous post's example, except that it writes messages to a temporary log file instead of a persistent one.

On Windows, the simplest version of the monitoring system sends messages to a dialog box. StartMonitoring creates the dialog box and starts it running in a separate thread. MonitorHandler calls LogMessage to add lines of text to the ListBox contained in the dialog box. StopMonitoring waits for the dialog box's thread to terminate. LogMessage sends the lines of text to the ListBox using the Windows messaging system. Since each item in a list box only displays a single line of text, LogMessage sends multi-line text as a sequence of separate messages, one message per line.

int LogMessage(const char *s, int msgType) {
    char msg[1024]; char *pos = &msg[0]; char *cr = NULL; msg[0] = '\0';
    strncat(msg, s, 1024);

    while (pos != NULL && *pos != '\0') {
       cr = strchr(pos, '\n');
        if (cr != NULL)
            *cr = '\0';
        pos = CleanText(pos); if (strlen(pos) > 0)
            SendMessage(loggerDialog, WM_COMMAND, MAKEWPARAM(msgType,
            0),
                        (LPARAM)pos);
        pos = cr+1;
    } return strlen(s);
}

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

On UNIX machines and when running the dialog box version (monitorWinDialog) on a Windows machine, you'll see a window pop up and display the messages. To view the output of the event log version on a Windows PC, you'll need to open the Event Logger and click on the Application category.

Next time: handling errors differently than information and warning messages. Until then, let me know how your application makes sure error messages get delivered to users. Do you need a permanent record, like the Windows event logger provides?

Get the MATLAB code

Published with MATLAB® 7.13

Contents

Adorned Ginger Man

Here's Rafael's adorned gingerman.

gingerR

Here's a comment from Rafael for this submission:

I decided to make a pixelated gingerbread man with the output of the chaotic map, using a lemniscate of Gerono as his green bow tie :)

The Code

I took the liberty of adding a small number of comments to Rafael's code to highlight some parts of the plot.

type gingerR

function gingerR
%function pixelatedGBM from Rafael
    rng(42) % could comment this out to get other "men"
    [x,y] = gingerbreadman;
    % scale and rotate our gingerbread man
    r = [x y] * [cosd(135) -sind(135); sind(135) cosd(135)];
    minR = min(r); maxR = max(r);
    r = (r - repmat(minR,size(r,1),1))./repmat(maxR-minR,size(r,1),1);
    r(:,2) = r(:,2).^1.5;
    
    % create a pixel representation of it
    N = 25;
    b = linspace(0,1,N);
    dif = b(2)-b(1);
    [xb,yb] = meshgrid(b,b);
    
    C = zeros(N);
    x = r(:,1); y = r(:,2);
    for i = 1:numel(xb)
        C(i) = length(find(x >= xb(i) & x < xb(i) +dif & y >= yb(i) ...
            & y < yb(i)+dif));
    end
    C = reshape(C,size(xb));
    % smooth a little for better results
    smooth = @(A,L) ((eye(size(A,1)) + ...
        L ^ 2 * diff(eye(size(A,1)),2)' * diff(eye(size(A,1)),2) + ...
        2 * L * diff(eye(size(A,1)))' * diff(eye(size(A,1)))) \ A);
    D = smooth(smooth(C,0.1)',0.1)';
    
    % let's draw it :)
    set(figure,'Position',[0 0 300 400],'Color',[1 1 1])
    movegui(gcf,'center')
    set(surf(xb,yb,zeros(size(D)),D),'ZData',xb.*0-0.01);
    view(2); shading flat; grid off; axis off equal;
    colormap([1 1 1; pink(19)])
    hold on
    t = linspace(0,2*pi,50);
    % bowtie
    set(fill((sin(t))/9+0.5,(sin(2*t))/18+0.65,[0 .8 0]),...
        'EdgeAlpha',0)
    set(fill((sin(t))/30+0.5,(cos(t))/30+0.65,[0 .9 0]),...
        'EdgeAlpha',0)
    t = linspace(0,2*pi,5);
    % buttons
    set(fill((sin(t))/20+0.5,(cos(t))/20+0.5,[.8 0 0]),'EdgeAlpha',0)
    set(fill((sin(t))/20+0.5,(cos(t))/20+0.3,[.8 0 0]),'EdgeAlpha',0)
    % sugar
    plot3(r(1:10:end,1),r(1:10:end,2),r(1:10:end,1)*0-0.005,'.',...
        'MarkerSize',3,'MarkerEdgeColor',[1 1 1]);
end

Thanks for Participating

I really appreciate all the effort the contributors made. The decision was tough. I encourage you all to grab some of the contributions from the comments in the previous post and play with them yourselves.

Let me also take this moment to wish all of you a happy, healthy holiday season and new year!

Get the MATLAB code

Published with MATLAB® 7.13

Contents

Gingerbreadman Map Equations

The equations for the gingerbreadman map look simple enough. For any given point in space: , define the next point in the sequence by . I now show the same equations in this MATLAB code, accounting for some initial values x0,y0 to seed the calculation.

dbtype 17:25 gingerbreadman

17    
18    % main calculation
19    for i = 1:n
20        if i == 1
21            x(i) = 1 - y0 + abs(x0);
22            y(i) = x0;
23        else
24            x(i) = 1 - y(i-1) + abs(x(i-1));
25            y(i) = x(i-1);

Plot Results

I want to show several results from different starting values for x0,y0. First I'll set the random number generator seed for repeatable results.

rng(42)
gingerbreadman

gingerbreadman

gingerbreadman

gingerbreadman

gingerbreadman

Code Listing

Here's the full code listing for gingerbreadman. You can see it allows you to specify the initial conditions, and return points and initial conditions should you choose. With no output arguments, gingerbreadman creates plots like you've been seeing.

type gingerbreadman

function [xout,yout, x0, y0] = gingerbreadman(x0,y0)
%  Gingerbreadman map producing a chaotic 2-D map.

%  Copyright 2011 The MathWorks, Inc.

% if not enough inputs, assign random numbers
if nargin < 2
    x0 = randn();
    y0 = randn();
end

% iteration counter
n = 10000;

x = zeros(n,1);
y = zeros(n,1);

% main calculation
for i = 1:n
    if i == 1
        x(i) = 1 - y0 + abs(x0);
        y(i) = x0;
    else
        x(i) = 1 - y(i-1) + abs(x(i-1));
        y(i) = x(i-1);
    end
end

% if output is requested, return gingerbread x,y values and
% x0, y0 initial conditions
%
% otherwise plot results
if nargout > 0
    xout = x;
    yout = y;
else
    scatter(x, y, '.');
end

References

Make the Plot Prettier!

Instead of using the plotting function scatter, post your thoughts (in code) here for a chance to win some MATLAB bling. I will look at entries up through (was: Sunday, November 27, 2011) Wednesday December 14, 2011 and announce the winner (the one whose plot I find most interesting) shortly after that.

Get the MATLAB code

Published with MATLAB® 7.13

Contents

Recreate Data Set 1 from the Previous Post

clear all
clc
rng(1998);
mu = [0 0 0 0 0 0 0 0];
i = 1:8;
matrix = abs(bsxfun(@minus,i',i));
covariance = repmat(.5,8,8).^matrix;
X = mvnrnd(mu, covariance, 20);
Beta = [3; 1.5; 0; 0; 2; 0; 0; 0];
ds = dataset(Beta);
Y = X * Beta + 3 * randn(20,1);
b = regress(Y,X);
ds.Linear = b;

Use Lasso to Fit the Model

The syntax for the lasso command is very similar to that used by linear regression. In this line of code, I am going estimate a set of coefficients B that models Y as a function of X. To avoid over fitting, I'm going to apply five-fold cross validation.

[B Stats] = lasso(X,Y, 'CV', 5);

When we perform a linear regression, we generate a single set of regression coefficients. By default lasso will create 100 different models. Each model was estimated using a slightly larger . All of the model coefficients are stored in array B. The rest of the information about the model is stored in a structure named Stats.

Let's look at the first five sets of coefficients inside of B. As you traverse the rows you can see that as increases, the value model coefficients are usually shrinking towards zero.

disp(B(:,1:5))
disp(Stats)

       3.9147       3.9146       3.9145       3.9143       3.9142
      0.13502      0.13498      0.13494      0.13488      0.13482
      0.85283      0.85273      0.85262      0.85247      0.85232
     -0.92775     -0.92723      -0.9267       -0.926     -0.92525
       3.9415       3.9409       3.9404       3.9397       3.9389
      -2.2945       -2.294      -2.2936       -2.293      -2.2924
       1.3566       1.3567       1.3568       1.3569       1.3571
     -0.14796     -0.14803      -0.1481     -0.14821     -0.14833
         Intercept: [1x100 double]
            Lambda: [1x100 double]
             Alpha: 1
                DF: [1x100 double]
               MSE: [1x100 double]
    PredictorNames: {}
                SE: [1x100 double]
      LambdaMinMSE: 0.585
         Lambda1SE: 1.6278
       IndexMinMSE: 78
          Index1SE: 89

Create a Plot Showing Mean Square Error Versus Lambda

The natural question to ask at this point in time is "OK, which of these 100 different models should I use?". We can answer that question using lassoPlot. lassoPlot generates a plot that displays the relationship between and the cross validated mean square error (MSE) of the resulting model. Each of the red dots shows the MSE for the resulting model. The vertical line segments stretching out from each dot are error bars for each estimate.

You can also see a pair of vertical dashed lines. The line on the right identifies the value that minimizes the cross validated MSE. The line on the left indicates the highest value of whose MSE is within one standard error of the minimum MSE. In general, people will chose the that minimizes the MSE. On occasion, if a more parsimonious model is considered particularly advantageous, a user might choose some other value that falls between the two line segments.

lassoPlot(B, Stats, 'PlotType', 'CV')

Use the Stats Structure to Extract a Set of Model Coefficients.

The value that minimizes the MSE is stored in the Stats structure. You can use this information to index into Beta and extract the set of coefficients that minimize the MSE.

Much as in the feature selection example, we can see that the lasso algorithm has eliminated four of the five distractors from the resulting model. This new, more parsimonious model will be significantly more accurate for prediction than a standard linear regression.

ds.Lasso = B(:,Stats.IndexMinMSE);
disp(ds)

    Beta    Linear      Lasso   
      3       3.5819      3.0591
    1.5      0.44611      0.3811
      0      0.92272    0.024131
      0     -0.84134           0
      2       4.7091      1.5654
      0      -2.5195           0
      0      0.17123      1.3499
      0     -0.42067           0

Run a Simulation

Here, once again, it's very dangerous to base any kind of analysis on a single observation. Let's use a simulation to compare the accuracy of a linear regression with the lasso. We'll start by preallocating some variables.

MSE = zeros(100,1);
mse = zeros(100,1);
Coeff_Num = zeros(100,1);
Betas = zeros(8,100);
cv_Reg_MSE = zeros(1,100);

Next, we'll generate 100 different models and estimate the number of coefficients contained in the lasso model as well as the difference in the cross validated MSE between a standard linear regression and the lasso model.

As you can see, on average, the lasso model only contains 4.5 terms (the standard linear regression model includes 8). More importantly, the cross validated MSE for the linear regression model is about 30% larger than that generated from the lasso. This is an incredibly powerful results. The lasso algorithm is every bit as easy to apply as standard linear regression, however, it offers significant improvements in predictive accuracy compared to regression.

rng(1998);

for i = 1 : 100

    X = mvnrnd(mu, covariance, 20);
    Y = X * Beta + randn(20,1);

    [B Stats] = lasso(X,Y, 'CV', 5);
    Betas(:,i) = B(:,Stats.IndexMinMSE) > 0;
    Coeff_Num(i) = sum(B(:,Stats.IndexMinMSE) > 0);
    MSE(i) = Stats.MSE(:, Stats.IndexMinMSE);

    regf = @(XTRAIN, ytrain, XTEST)(XTEST*regress(ytrain,XTRAIN));
    cv_Reg_MSE(i) = crossval('mse',X,Y,'predfun',regf, 'kfold', 5);

end

Number_Lasso_Coefficients = mean(Coeff_Num);
disp(Number_Lasso_Coefficients)

MSE_Ratio = median(cv_Reg_MSE)/median(MSE);
disp(MSE_Ratio)

         4.57
       1.2831

Choosing the Best Technique

Regularization methods and feature selection techniques both have unique strengths and weaknesses. Let's close this blog posting with some practical guidance regarding pros and cons for the various techniques.

Regularization techniques have two major advantages compared to feature selection.

However, feature selection methods also have their advantages

With this said and done, each of the three regularization techniques also offers its own unique advantages and disadvantages.

Let's assume that you are running a cancer research study.

Sequential feature selection is completely impractical with this many different variables. You can't use ridge regression because it won't force coefficients completely to zero quickly enough. At the same time, you can't use lasso since you might need to identify more than 500 different genes. The elastic net is one possible solution.

Conclusion

If you'd like more information on this topic there is a MathWork's webinar titled Computational Statistics: Feature Selection, Regularization, and Shrinkage which provides a more detailed treatment of these topics.

In closing, I'd like to ask you whether any of you have practical examples applying feature selection or regularization algorithms in your work?

If so, please post here here.

Get the MATLAB code

Published with MATLAB® 7.13

Contents

Create a Dataset

To start with, I am going to generate a dataset that is deliberately engineered to be very difficult to fit using traditional linear regression. The dataset has three characteristics that present challenges of linear regression

clear all
clc
rng(1998);

Create eight X variables. The mean of each variable will be equal to zero.

mu = [0 0 0 0 0 0 0 0];

The variables are correlated with one another. The covariance matrix is specified as

i = 1:8;
matrix = abs(bsxfun(@minus,i',i));
covariance = repmat(.5,8,8).^matrix;

Use these parameters to generate a set of multivariate normal random numbers.

X = mvnrnd(mu, covariance, 20);

Create a hyperplane that describes Y = f(X) and add in a noise vector. Note that only three of the eight predictors have coefficient values > 0. This means that five of the eight predictors have zero impact on the actual model.

Beta = [3; 1.5; 0; 0; 2; 0; 0; 0];
ds = dataset(Beta);
Y = X * Beta + 3 * randn(20,1);

For those who are interested, this example was motivated by Robert Tibshirani's classic 1996 paper on the lasso:

Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. J. Royal. Statist. Soc B., Vol. 58, No. 1, pages 267-288).

Let's model Y = f(X) using ordinary least squares regression. We'll start by using traditional linear regression to model Y = f(X) and compare the coefficients that were estimated by the regression model to the known vector Beta.

b = regress(Y,X);
ds.Linear = b;
disp(ds)

    Beta    Linear  
      3       3.5819
    1.5      0.44611
      0      0.92272
      0     -0.84134
      2       4.7091
      0      -2.5195
      0      0.17123
      0     -0.42067

As expected, the regression model is performing very poorly.

Let's see if we can improve on these results using a technique called "All possible subset regression". First, I'm going to generate 255 different regession models which encompass every possible combinate of the eight predictors. Next, I am going to measure the accuracy of the resulting regression using the cross validated mean squared error. Finally, I'm going to select the regression model that produced the most accurate results and compare this to our original regression model.

Create an index for the regression subsets.

index = dec2bin(1:255);
index = index == '1';
results = double(index);
results(:,9) = zeros(length(results),1);

Generate the regression models.

for i = 1:length(index)
    foo = index(i,:);
    rng(1981);
    regf = @(XTRAIN, YTRAIN, XTEST)(XTEST*regress(YTRAIN,XTRAIN));
    results(i,9) = crossval('mse', X(:,foo), Y,'kfold', 5, 'predfun',regf);
end

Sort the outputs by MSE. I'll show you the first few.

index = sortrows(results, 9);
index(1:10,:)

ans =
  Columns 1 through 5
            1            1            0            0            1
            1            1            0            1            1
            1            0            1            0            1
            1            1            0            0            1
            1            0            0            0            1
            1            1            0            1            1
            1            0            0            1            1
            1            0            0            0            1
            1            0            1            0            1
            1            0            0            1            1
  Columns 6 through 9
            1            0            0       10.239
            1            0            0       10.533
            1            0            0       10.615
            0            0            0       10.976
            1            0            0       11.062
            0            0            0       11.257
            1            0            0       11.641
            0            0            0       11.938
            0            0            0       12.068
            0            0            0       12.249

Compare the Results

beta = regress(Y, X(:,logical(index(1,1:8))));
Subset = zeros(8,1);
ds.Subset = Subset;
ds.Subset(logical(index(1,1:8))) = beta;
disp(ds)

    Beta    Linear      Subset 
      3       3.5819     3.4274
    1.5      0.44611     1.1328
      0      0.92272          0
      0     -0.84134          0
      2       4.7091     4.1048
      0      -2.5195    -2.0063
      0      0.17123          0
      0     -0.42067          0

The results are rather striking.

All possible subset regression appears to have generated a significantly better model. This R^2 value for this regression model isn't as good as the original linear regression; however, if we're trying to generate predictions from a new data set, we'd expect this model to perform significantly better.

Perform a Simulation

It's always dangerous to rely on the results of a single observation. If we were to change either of the noise vectors that we used to generate our original dataset we might end up with radically different results. This next block of code will repeat our original experiment a hundred times, each time adding in different noise vectors. Our output will be a bar chart that shows how frequently each of the different variables was included in the resulting model. I'm using parallel computing to speed things up though you could perform this (more slowly) in serial.

matlabpool open

sumin = zeros(1,8);
parfor j=1:100
    mu = [0 0 0 0 0 0 0 0];
    Beta = [3; 1.5; 0; 0; 2; 0; 0; 0];
    i = 1:8;
    matrix = abs(bsxfun(@minus,i',i));
    covariance = repmat(.5,8,8).^matrix;
    X = mvnrnd(mu, covariance, 20);
    Y = X * Beta + 3 * randn(size(X,1),1);

    index = dec2bin(1:255);
    index = index == '1';
    results = double(index);
    results(:,9) = zeros(length(results),1);

    for k = 1:length(index)
        foo = index(k,:);
        regf = @(XTRAIN, YTRAIN, XTEST)(XTEST*regress(YTRAIN,XTRAIN));
        results(k,9) = crossval('mse', X(:,foo), Y,'kfold', 5, ...
            'predfun',regf);
    end

% Sort the outputs by MSE
index = sortrows(results, 9);
sumin = sumin + index(1,1:8);
end  % parfor

matlabpool close

Starting matlabpool using the 'local4' configuration ... connected to 3 labs.
Sending a stop signal to all the labs ... stopped.

%Plot results.
bar(sumin)

The results clearly show that "All possible subset selection" is doing a really good job identifying the key variables that drive our model. Variables 1, 2, and 5 show up in most of the models. Individual distractors work their way into the models as well, however, their frequency is no where near as high as the real predictors.

Introducing Sequential Feature Selection

All possible subset selection suffers from one enormous limitation. It's computationally infeasible to apply this technique to datasets that contain large numbers of (potential) predictors. Our first example had eight possible predictors. Testing all possible subsets required generate (2^8)-1 = 255 different models. Suppose that we were evaluating a model with 20 possible predictors. All possible subsets would require use to test 1,048,575 different models. With 40 predictors we'd be looking at a billion different combinations.

Sequential feature selection is a more modern approach that tries to define a smart path through the search space. In turn, this should allow us to identify a good set of cofficients but ensure that the problem is still computationally feasible.

Sequential feature selection works as follows:

% This next bit of code will show you how to use sequential feature
% selection to solve a significantly larger problem.

Generate a Data Set

clear all
clc
rng(1981);

Z12 = randn(100,40);
Z1 = randn(100,1);
X = Z12 + repmat(Z1, 1,40);
B0 = zeros(10,1);
B1 = ones(10,1);
Beta = 2 * vertcat(B0,B1,B0,B1);
ds = dataset(Beta);

Y = X * Beta + 15*randn(100,1);

Use linear regression to fit the model.

b = regress(Y,X);
ds.Linear = b;

Use sequential feature selection to fit the model.

opts = statset('display','iter');

fun = @(x0,y0,x1,y1) norm(y1-x1*(x0\y0))^2;  % residual sum of squares
[in,history] = sequentialfs(fun,X,Y,'cv',5, 'options',opts)

Start forward sequential feature selection:
Initial columns included:  none
Columns that can not be included:  none
Step 1, added column 26, criterion value 868.909
Step 2, added column 31, criterion value 566.166
Step 3, added column 9, criterion value 425.096
Step 4, added column 35, criterion value 369.597
Step 5, added column 16, criterion value 314.445
Step 6, added column 5, criterion value 283.324
Step 7, added column 18, criterion value 269.421
Step 8, added column 15, criterion value 253.78
Step 9, added column 40, criterion value 242.04
Step 10, added column 27, criterion value 234.184
Step 11, added column 12, criterion value 223.271
Step 12, added column 17, criterion value 216.464
Step 13, added column 14, criterion value 212.635
Step 14, added column 22, criterion value 207.885
Step 15, added column 7, criterion value 205.253
Step 16, added column 33, criterion value 204.179
Final columns included:  5 7 9 12 14 15 16 17 18 22 26 27 31 33 35 40 
in =
  Columns 1 through 12
     0     0     0     0     1     0     1     0     1     0     0     1
  Columns 13 through 24
     0     1     1     1     1     1     0     0     0     1     0     0
  Columns 25 through 36
     0     1     1     0     0     0     1     0     1     0     1     0
  Columns 37 through 40
     0     0     0     1
history = 
      In: [16x40 logical]
    Crit: [1x16 double]

I've configured sequentialfs to show each iteration of the feature selection process. You can see each of the variables as it gets added to the model, along with the change in the cross validated mean squared error.

Generate a linear regression using the optimal set of features

beta = regress(Y,X(:,in));

ds.Sequential = zeros(length(ds),1);
ds.Sequential(in) = beta

ds = 
    Beta    Linear      Sequential
    0       -0.74167          0   
    0       -0.66373          0   
    0         1.2025          0   
    0        -1.7087          0   
    0         1.6493     2.5441   
    0         1.3721          0   
    0        -1.6411    -1.1617   
    0        -3.5482          0   
    0         4.1643     2.9544   
    0        -1.9067          0   
    2        0.28763          0   
    2         1.7454     2.9983   
    2         0.1005          0   
    2         2.2226     3.5806   
    2         4.4912     4.0527   
    2         4.7129     5.3879   
    2         0.9606     1.9982   
    2         5.9834     5.1789   
    2         3.9699          0   
    2        0.94454          0   
    0          2.861          0   
    0        -2.1236    -2.4748   
    0        -1.6324          0   
    0         1.8677          0   
    0       -0.51398          0   
    0         3.6856     4.0398   
    0        -3.3103    -4.1396   
    0       -0.98227          0   
    0         1.4756          0   
    0       0.026193          0   
    2         2.9373     3.8792   
    2        0.63056          0   
    2         2.3509     1.9472   
    2         0.8035          0   
    2         2.9242       4.35   
    2         0.4794          0   
    2        -0.6632          0   
    2         0.6066          0   
    2       -0.70758          0   
    2         4.4814     3.8164   

As you can see, sequentialfs has generated a much more parsimonious model than the ordinary least squares regression. We also have a model that should generate significantly more accurate predictions. Moreover, if we needed to (hypothetically) extract the five most important variables we could simple grab the output from Steps 1:5.

Run a Simulation

As a last step, we'll once again run a simulation to benchmark the technique. Here, once again, we can see that the feature selection technique is doing a very good job identifying the true predictors while screening out the distractors.

matlabpool open

sumin = 0;
parfor j=1:100
    Z12 = randn(100,40);
    Z1 = randn(100,1);
    X = Z12 + repmat(Z1, 1,40);
    B0 = zeros(10,1);
    B1 = ones(10,1);
    Beta = 2 * vertcat(B0,B1,B0,B1);
    Y = X * Beta + 15*randn(100,1);

    fun = @(x0,y0,x1,y1) norm(y1-x1*(x0\y0))^2;  % residual sum of squares
    [in,history] = sequentialfs(fun,X,Y,'cv',5);

    sumin = in + sumin;
end

bar(sumin)

matlabpool close

Starting matlabpool using the 'local4' configuration ... connected to 3 labs.
Sending a stop signal to all the labs ... stopped.

One last comment: As you can see, neither sequentialfs nor all possible subset selection is perfect. Neither of these techniques was able to capture all of the true variables while excluding the distractors. It's important to recognize that these data sets were deliberately designed to be difficult to fit. [We don't want to benchmark performance using "easy" problems where any old technique might succeed.]

Next week's blog post will focus on a radically different way to solve these same types of problems. Stay tuned for a look at regularization techniques such as ridge regression, lasso and the elastic net.

Here are a couple questions that you might want to consider: # All possible subset methods are O(2^n). What is the order of complexity of sequential feature selection? # Other than scaling issues, can you think of any possible problems with sequential feature selection? Are there situations where the algorithm might miss important variables?

Get the MATLAB code

Published with MATLAB® 7.13

Contents

A Brief History of MATLAB to C

In April, 2011, MathWorks introduced MATLAB Coder as a stand-alone product to generate C code from MATLAB code. This tool lets user generate readable, portable, and customizable C code from their MATLAB algorithms.

Many astute readers will notice that C code generation from MATLAB isn't really brand new - and that we've had this capability to generate C code from MATLAB for quite some time now. Yes, that's true - we've been incubating this technology for quite some time till we felt it was ready to debut as a stand alone tool. Here's a timeline of this technology over the past few years:

A Simple Example

Let me introduce the basics of using MATLAB Coder through a very simple example that multiplies two variables. Let's generate C code from the following MATLAB function that multiplies two inputs:

dbtype simpleProduct.m

1     function c = simpleProduct(a,b) %#codegen
2     c = a*b;

To generate C code from this function using MATLAB Coder, I first have to specify the size and data types of my inputs - I can do this through the MATLAB Coder UI (shown in the section below) where I specify my inputs as a [1x5] and [5x2] single precision matrices, and subsequently generate the following C code:

dbtype simpleProduct.c

1     #include "simpleProduct.h"
2     
3     void simpleProduct(const real32_T a[5], const real32_T b[10], real32_T c[2])
4     {
5       int32_T i0;
6       int32_T i1;
7       for (i0 = 0; i0 < 2; i0++) {
8         c[i0] = 0.0F;
9         for (i1 = 0; i1 < 5; i1++) {
10          c[i0] += a[i1] * b[i1 + 5 * i0];
11        }
12      }
13    }

MATLAB is a polymorphic language which means a single function, such as simpleProduct, can accept input arguments of different size and data types and output correct results. This function can behave as a simple scalar product, a dot product, or a matrix multiplier depending on what inputs you pass.

Languages like C are said to be “stongly typed,” which requires you to create a different version of a function for every variation of input data size and types. Hence, when you write C code to implement simpleProduct, you have to know ahead of time the sizes and the data types of your inputs so you can implement the right variant.

While MATLAB Coder helps you move from the highly flexible world of MATLAB to a strongly typed world of C, you still have to specify all of the constraints C expects. You do this in MATLAB Coder by creating a MATLAB Coder project file (simpleProduct.prj in this case) where you can specify the various code generation parameters including the sizes and data types of your inputs as shown below.

You can also edit the default configuration parameters by clicking on the 'More Settings' link which appears on the 'Build' tab of this project UI .

In the example above, I turned off a few default options such as including comments in the generated code just so you can see how compact the code is. You not only get the option to generate comments in the resulting C code, but also inlcude the original MATLAB code as comments in the corresponding sections of the generated code. This helps with traceability of your C code to your original algorithms.

Detailed MATLAB Coder Worfklow

The simple example above quickly illustrates the process of generating code with MATLAB coder and shows how the resulting C code looks.

Naturally, your real-world functions are going to be much more involved and may run into hundreds or even thousands of lines of MATLAB Code. To help you handle that level of complexity, you would need an iterative process, like the 3-step workflow described here, that guides you through the task of code generation incrementally:

I'll now demonstrate this concept using a slightly more involved example. The following MATLAB code implements the Newton-Raphson numerical technique for computing the n-th root of a real valued number.

dbtype newtonSearchAlgorithm.m

1     function [x,h] = newtonSearchAlgorithm(b,n,tol) %#codegen
2     % Given, "a", this function finds the nth root of a
3     % number by finding where: x^n-a=0.
4     
5         notDone = 1;
6         aNew    = 0; %Refined Guess Initialization
7         a       = 1; %Initial Guess
8         count     = 0;
9         h(1)=a;
10        while notDone
11            count = count+1;
12            [curVal,slope] = fcnDerivative(a,b,n); %square
13            yint = curVal-slope*a;
14            aNew = -yint/slope; %The new guess
15            h(count)=aNew;
16            if (abs(aNew-a) < tol) %Break if it's converged
17                notDone = 0;
18            elseif count>49 %after 50 iterations, stop
19                notDone = 0;
20                aNew = 0;
21            else
22                a = aNew;
23            end
24        end
25        x = aNew;
26        
27    function [f,df] = fcnDerivative(a,b,n)
28    % Our function is f=a^n-b and it's derivative is n*a^(n-1).
29    
30        f  = a^n-b;
31        df = n*a^(n-1);

Our first step towards generating C code from this file is to prepare it for code generation. For code generation, each variable inside your MATLAB code must be intialized, which means specifying its size and data type. In this case, I'll choose to initialize the variable h to a static maximum size, h = zeros(1,50), as it doesn't grow beyond a length of 50 inside the for loop.

You can invoke the code generation function using the GUI or the command line (through the codegen command). Without worrying about the details of this approach, let's look at the generated C code:

dbtype newtonSearchAlgorithmMLC.c

1     #include "newtonSearchAlgorithmMLC.h"
2     
3     void newtonSearchAlgorithmMLC(real_T b, real_T n, real_T tol, real_T *x, real_T
4       h[50])
5     {
6       int32_T notDone;
7       real_T a;
8       int32_T count;
9       real_T u1;
10      real_T slope;
11      notDone = 1;
12      *x = 0.0;
13      a = 1.0;
14      count = -1;
15      memset((void *)&h[0], 0, 50U * sizeof(real_T));
16      h[0] = 1.0;
17      while (notDone != 0) {
18        count++;
19        u1 = n - 1.0;
20        u1 = pow(a, u1);
21        slope = n * u1;
22        u1 = pow(a, n);
23        *x = -((u1 - b) - slope * a) / slope;
24        h[count] = *x;
25        if (fabs(*x - a) < tol) {
26          notDone = 0;
27        } else if (count + 1 > 49) {
28          notDone = 0;
29          *x = 0.0;
30        } else {
31          a = *x;
32        }
33      }
34    }

The subfunction fcnDerivative is inlined in the generated C code. You can choose to not inline the code by putting this command coder.inline('never') in the MATLAB file.

Support for Dynamic Sizing

If you have variable in your MATLAB code that needs to vary its size during execution, you can choose to deal with this in three different ways in the generated C code:

The last two options can be enabled by turning on the variable-sizing feature in the configuration parameters. However, do remember that enabling this option also makes the resulting C code bigger in size - so if you can avoid it, you can get much more compact and possibly more efficient code.

Use Cases of MATLAB Coder

The primary use of MATLAB Coder is to enable algorithm developers and system engineers working in MATLAB to quickly generate readable and portable C code. The generated C code can be used in different ways supporting different workflows. Here are a few use cases of MATLAB Coder:

I won't be talking about creating standalone executables or creating libraries from MATLAB Coder in this blog; I'll, however, say a few words on generating MEX functions.

Generating MEX Functions for Simulation Acceleration

The MEX interface enables you to integrate C code into MATLAB. One common use of MEX is to speed up performance bottlenecks in MATLAB code by re-writing them in C and executing them back in MATLAB as MEX functions. MATLAB Coder can save you time and effort by automatically generating MEX functions from your MATLAB code (and allowing you to import legacy C code easily as well).

The speed-up you get through automatic MEX generation can vary quite a bit depending on the application. In some cases, you may not get any speed up at all, or possibly even a slow-down, as MATLAB language has gotten quite smart and efficient in computing many built-in functions by automatically taking advantage of processor specific routines (such as Intel Performance Primitives) and multithreading to utilize multiple cores.

A few guidelines that you can use to determine if your algorithm is a good candidate for speed-up with MEX are:

In these cases, you can expect to see a speed-up. Let's look at a realistic example that illustrates this concept.

The example I chose here is one that my colleague Sarah Zaranek had highlighted in her guest blog on vectorization.

I modified the original MATLAB script from that article into a function as only functions are supported by MATLAB Coder. The other change that I made to the original script was to initialize the output variables subs and vals (as all variables in your code has to be defined/intialized once when used with MATLAB Coder). And by turning on the variable sizing feature (using dynamic memory allocation), I didn't have to preallocate it to a sufficiently large enough size to accomodate its growth, which would use more memory (often a lot) than needed.

The following statements create and use a MATLAB Coder configuration object with support for variable size arrays and dynamic memory allocation to generate a MEX function.

% Create a MEX configuration object
cfg = coder.config('mex');
% Turn on dynamic memory allocation
cfg.DynamicMemoryAllocation = 'AllVariableSizeArrays';
% Generate MEX function
codegen -config cfg gridSpeedUpFcn -args {10}
disp('MEX generation complete!')

Warning: There is no short path form of 'H:\Documents\LOREN\MyJob\Art of
MATLAB' available. Short path names do not include embedded spaces. The
short path name feature may be disabled in your operating system. Short
path names are required for successful code generation, therefore it may
not be possible to successfully complete the code generation process. To
avoid this problem, enable short path name support in your operating
system, or do not attempt code generation in paths containing spaces. 
MEX generation complete!

This generates a MEX function in the current folder. I'll use tic, and toc to estimate the execution times of the original MATLAB function and its MEX version.

I run each function multiple times inside a for loop and use only the last few runs for our computation time calculation so we can minimize the effects of initialization and caching.

MToc = zeros(1,30);
MexToc = zeros(1,30);
for j = 1:30
    tic; gridSpeedUpFcn(10);
    MToc(j) = toc;
end

for j = 1:30
    tic; gridSpeedUpFcn_mex(10);
    MexToc(j) = toc;
end

eff_M_time = mean(MToc(1,21:30));
eff_Mex_time = mean(MexToc(1,21:30));

disp(['MATLAB code execution time (with nx=10): ', num2str(eff_M_time)])
disp(['MEX code executing time (with nx=10): ', num2str(eff_Mex_time)])

disp([' Speed up factor: ', num2str(eff_M_time/eff_Mex_time)]);

MATLAB code execution time (with nx=10): 0.31652
MEX code executing time (with nx=10): 0.0027057
 Speed up factor: 116.9806

The speed up we see in this example is more pronounced for larger values of nx (and would vary between different computers).

Vectorization and Code Generation

Sarah, in her blog article, explains quite nicely on how to vectorize this code using meshgrid and ndgrid to explode the variables in support of vectorization, and get similar levels of speed up.

Sarah paid a small price in code readability for her vectorization efforts, but a much larger price in memory consumption. This points to a common trade-off when vectorizing MATLAB code: execution time vs. memory consumption. The temporary variables which enable us to remove loops through matrix math are often very large, effectively limiting the number of iterations you can take before running out of memory. For instance, if we run her algorithm for a 1,000x1,000 grid instead of a 100x100 grid, we would create 9 double precision variables that are about 7.5 GB each!

Code generation can effectively deliver similar speed-up benefits (and with less effort in some cases), as it's trying to explicitly achieve the same for loop operations in C code that a vectorized code achieves. However, the downside for most users with the code generation approach is that it's available with a separate product (higher cost), and only MATLAB code that supports code generation can automatically be converted to MEX as well.

Customizing and Optimizing the Generated Code

Now with most of the basics of code generation out of the way, a few questions might be on your mind:

Learning More About MATLAB Coder

You can find more information about MATLAB Coder on the product page, which includes demos and recorded webinars. The product documentation is also an excellent resource to find answers to questions you might have when starting with this product.

Please feel free to share your thoughts and comments on topics discussed here in this posting here.

Get the MATLAB code

Published with MATLAB® 7.13

Contents

Two Output Streams

Programs running in interactive MATLAB send error, warning and informational text messages to the MATLAB command window. By default, deployed applications output text to the standard output and standard error streams.

In many cases, particularly for console applications on Windows and UNIX, this works well enough: the output appears in the DOS window or UNIX terminal. However, non-console applications (and any type of shared library or builder component) may be unable to access either of the standard output streams; messages sent to those streams will be lost. To avoid losing messages when the standard streams are unavailable, deployed applications register output and error message handlers. A message handler directs formatted text messages to a persistent or visible location, typically a log file or an output window.

In this post, I'll show you how to send both error and informational messages to a log file.

A Function with Output

To demonstrate output redirection, I wrote the sendmessages function, which displays an informational message and a warning, and then issues an error:

function sendmessages(info, warn, err)
    disp(info);
    warning(warn);
    error(err);

In order to follow the rest of this posting, you'll need the source code from MATLAB Central, particularly the C main program, log.c. See logREADME in the ZIP file for a complete description of distributed files and step-by-step directions for building the application.

Where are the Handlers?

Create a C shared library from the example function with MATLAB Compiler:

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

Open the generated wrapper file, libmsgfcn.c, in a text editor. The wrapper file contains two default handlers, mclDefaultErrorHandler and mclDefaultPrintHandler. The generated code contains two library initialization functions: libmsgfcnInitialize and libmsgfcnInitializeWithHandlers. You redirect program output by passing your custom message handlers to libmsgfcnInitializeWithHandlers. All message handlers conform to the same interface:

 static int MessageHandler(const char *message)

You never need to call the message handler directly. The MATLAB Compiler Runtime (MCR) calls your message handler automatically.

Logging Messages to a File

A common use for a message handler is to record messages and the time they occurred in a log file. This type of message handler requires only a few lines of C code:

static FILE *logFile = NULL;

static int LogFileHandler(const char *s)
{
    if (logFile != NULL)
    {
        time_t now = time(NULL);
        fprintf(logFile, "%s", ctime(&now));
        fprintf(logFile, "%s\n", s);
    }
}

Your main program needs to open and close the log file, of course. Open the log file before calling the library initialization function and close it after calling mclTerminateApplication. See the full implementation in log.c.

Build the application with mbuild:

(UNIX): mbuild -g log.c -L. -lmsgfcn

(Windows): mbuild -g log.c libmsgfcn.lib

Run it from a DOS shell or UNIX terminal window:

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

(Windows): log "information" "look out!" "oops."

The example program creates the log file MessageLog.txt to capture program output.

Fri Aug 19 11:53:32 2011
information

Fri Aug 19 11:53:32 2011
Warning: lookout!
> In sendmessages at 3

Fri Aug 19 11:53:32 2011
Error using sendmessages (line 4)
oops.

Next time, I'll show you how to display log messages in a persistent window, and how to distinguish between errors and informational messages. In the meantime, let me know how your applications manage output. Does your code need to distinguish between program output and error messages?

Get the MATLAB code

Published with MATLAB® 7.13