Contents

2007 MATLAB World Tour - Australia and Taiwan

Come meet us during our visits to cities in Asia and the Pacific region for the 2007 MATLAB World Tour.

Korea

We will also be attending the MATLAB and Simulink Conference 2007 in Seoul, Korea later this month.

Manhattan Beach, California

On my way home from Asia and Australia, I will be join other colleagues to attend The MathWorks Aerospace and Defense Conference '07.

See You There!

I hope to meet many of you at these events. I plan to take some pictures and post some highlights, as I did from the trips I just completed in Europe.

Get the MATLAB code

Published with MATLAB® 7.4

Contents

Reasons to Keep inv

There are several reasons why we do keep the function inv. These reasons are:

Reasons to Not Use or Remove inv

There are good reasons to not use inv however. The main one is

Here is the second paragraph of the description for inv in the documentation:

In practice, it is seldom necessary to form the explicit inverse of a matrix. A frequent misuse of inv arises when solving the system of linear equations.

One way to solve this is with x = inv(A)*b. A better way, from both an execution time and numerical accuracy standpoint, is to use the matrix division operator x = A\b. This produces the solution using Gaussian elimination, without forming the inverse. See \ and / for further information.

According to Cleve, the best way to compute inv(X'*X) is to avoid computing X'*X, and use the economy-sized QR factorization of X instead.

  [Q,R] = qr(X,0)

  S = inv(R)

  S*S'

To see an extreme example, try

d = 1.e-10
X = [1 1; d 0; 0 d]

d =
  1.0000e-010
X =
    1.0000    1.0000
    0.0000         0
         0    0.0000

inv(X'*X)

Warning: Matrix is singular to working precision.
ans =
   Inf   Inf
   Inf   Inf

[Q,R] = qr(X,0)
S = inv(R)
S*S'

Q =
   -1.0000    0.0000
   -0.0000   -0.7071
         0    0.7071
R =
   -1.0000   -1.0000
         0    0.0000
S =
  1.0e+009 *
   -0.0000   -7.0711
         0    7.0711
ans =
  1.0e+019 *
    5.0000   -5.0000
   -5.0000    5.0000

The computed X' * X is exactly singular and does not have an inverse, but R does have an inverse. We're still using the inv function, but on a nonsingular matrix.

Your Thoughts?

So, we have the function inv that is ripe for misuse. What's the best way we can steer users to better solutions? Please post your thoughts here.

Get the MATLAB code

Published with MATLAB® 7.4

Contents

Italy

Scott and I started off in Milano.

Scott got this beautiful cup of coffee. It reminds me of this submission on the File Exchange.

Here's a shot of my master class. The class was very interactive, with some participants "earning" a MATLAB T-shirt for asking a good question or answering a complex question well.

Spain

Scott and I then went to Madrid.

More T-shirts rewarded for participation. (Hint: this is worth keeping in mind for the future, if you come to visit us on the tour elsewhere.)

France

Next stop, meeting Mike in Toulouse.

We met some customers throughout, including during the outdoors coffee break. Here's Jérôme Briot and Gaetan Guerin. Jérôme is an active participant on the MATLAB newsgroup when he is not doing his research in medical imaging.

We used Mike's clock as our timekeeper in France.

In Paris, we met Charles Lehalle who uses MATLAB heavily for his finance work, such as tools he's posted here.

Ambling around Paris, I saw tables of people playing chess in the park.

And these charming mice showed up at the Musée d'Orsay.

Switzerland

Next stop, Zurich, where Sam

and I met up with another great group of users.

Souvenirs

I didn't buy much on my trip, but here are a few delicacies that I brought home for gifts. Such beautiful packaging!

All in all, great memories, including those from meeting with so many curious customers!

Get the MATLAB code

Published with MATLAB® 7.4

Contents

What error Does

The main thing error is responsible for is issuing a message to the user that he or she has done something incorrectly. The first place you might see error checking in a function is looking at the input arguments. Using nargchk is a common way to see if there's a problem and then calling error afterwards.

Sample Snippet

Let's see one pattern. Here we are checking that the number of inputs, narg, lies between 2 and 7. If not, we'll issue an error.

narg = 4;
msg = nargchk(2,7,narg,'struct');
if ~isempty(msg)
    error(msg)
end

It turns out, this code is more complex than necessary. We can equivalently write

error(nargchk(2,7,narg,'struct'));

because the function error does not produce an error with empty input. This is a gotcha that people occasionally run into. You can't just use a simple statement like error('') to signal an error in your code.

Why Use the 'struct' Option?

I am using the 'struct' argument for nargchk so the output is a structure with both a message identifier and the message itself. Why? Because I, as the application developer, could choose to issue a different message or take a different action, e.g., using a try-catch construct. Using message identifiers for warning users about something allows the users to control that specific warning. Perhaps, after seeing it once, a user might want to turn that warning off, but would still like to be alerted of other possible issues.

Another gotcha

Users are also occasionally tripped by turning a particular warning off, running some other code, and then not expecting to see the turned off warning showing up when they use lastwarn to see if something has happened. Depending on the sequence of events, it's possible that lastwarn will return a message from a warning that was turned off. warning off ID simply turns off the display of the warning message.

Any More Message Wisdom?

If you have more message wisdom, post your thoughts here.

Get the MATLAB code

Published with MATLAB® 7.4

Contents

Recent Sample Question

Recently there was a question on the newsgroup about how to vectorize the access to elements in a structure array. And also got an email from a customer on this topic. Here is one customer's description of his problem:

    In a nutshell, what I am trying to do (in as few lines of code as
    possible) is:

    state = an array of structs, say N items each with (for example)
    components x and y.

    In my case 'state' is reasonably complicated and hence does not
    warrant the use of a simple 2 x N matrix.

    for i=2:N
       state(i)=markovmodel(state(i-1)); % 1. Access individual structs
    end

    plot(state.x) % 2. Access all of the entries of one element as  vector.

The Answer

Let's create an array of structures first.

states(1).x = 3;
states(2).x = 4;
states(3).x = 5;
states(1).y = 42;
states(2).y = 99;
states(3).y = 0;
states

states = 
1x3 struct array with fields:
    x
    y

Let's see what's in states.x

states.x

ans =
     3
ans =
     4
ans =
     5

With an array of structs, you can gather the x values using

allxvals = [states.x]

allxvals =
     3     4     5

This is because states.x produces a comma-separated list. How do I know this? Because, if I leave the semi-colon off, I see a bunch of output with ans =; in other words, I am getting multiple output values. How do I collect multiple values into an array? I can use either square brackets or curly braces.

Using the square brackets, [], just places that list into between the brackets, creating an array of one type, in this case, doubles. You can do likewise with struct arrays with fields that are more suitable for a cell array using curly braces, {} in place of the square brackets.

User Solution

If instead of using

    plot(state.x)
above, the user replaced this with

    plot([state.x])

the informationg gets plotted correctly.

Comments?

My question to you is why this question keeps getting asked. Is the concept unusual enough that people don't know how to even look for the information on-line? What can we do at The MathWorks to make this more visible? Please pass along your thoughts here.

Get the MATLAB code

Published with MATLAB® 7.4

Contents

Recent Example at The MathWorks

Recently at The MathWorks, someone wanted to share some files with another person so they could discuss the content. To do so, the author provided code that looked like this. Here YEAR is a variable of class char.

   load(['T:\1-36U6NN\secprd',YEAR,'.txt'])

You might wonder what the problem with this is. Well, for one thing, the author works on a Windows platform, but the recipient frequently does not. So how could this code have been written differently?

The T:\ drive maps to a file server that is accessible from any of our networked computers. So one possibility is to write code something like this:

      baseLocationWin = 'T:\';
      baseLocationOther = '//home/tester';
      if ispc  % true for windows platforms
          baseLocation = baseLocationWin;
      else
          baseLocation = baseLocationOther;
      end
      path1 = '1-36U6NN';
      YEAR = '2001';
      fname = ['secprd' YEAR '.txt'];
      fullfname = fullfile(baseLocation, path1, fname);

On Windows, where I ran this code, here's the file name.

disp(fullfname)

T:\1-36U6NN\secprd2001.txt

The equivalent UNIX pathname would be

disp([baseLocationOther '/' path1, '/' fname])

//home/tester/1-36U6NN/secprd2001.txt

Reference Material

Here you'll see a list of some MATLAB functions that are useful in writing machine independent applications when it comes to referring to files.

Comments

Have you been in this situation? Do you have other techniques to help users write robust code for handling files? If so, please post here.

Get the MATLAB code

Published with MATLAB® 7.4

Contents

Locations

The tour is a full-day event, planned for 17 (!) cities in Europe, Israel, Australia and Asia.

Format

The format is to have talks by MathWorks senior staff in the mornings, followed by master classes in the afternoons.

Who Should Attend

If there's a session near you, please register. I will be speaking in several of the locations and hope to see many of you there. It will be fun to put a face with a name.

Get the MATLAB code

Published with MATLAB® 7.4

Contents

Pitfalls

I mentioned one pitfall above and in an earlier post.

format short
for ind = 0:.1:1;
    if ind ~= fix(10*ind)/10
        disp(ind - fix(10*ind)/10)
    end
end

  5.5511e-017
  1.1102e-016
  1.1102e-016

If : were implemented naively, there would be another pitfall for users. What would happen if you want to produce some numbers between 1 - eps and 1 + eps, let's say with an increment of eps/4. Just doing repeated addition, we'd never get from the start to the end.

x = 1-eps
x1 = x + eps/4
isequal(x,x1)

x =
    1.0000
x1 =
    1.0000
ans =
     1

It's sometimes helpful to see the results in hex format to see the details. I'll show a few of the relevant numbers here.

format hex
x
x1
real1 = 1
eps1 = eps

x =
   3feffffffffffffe
x1 =
   3feffffffffffffe
real1 =
   3ff0000000000000
eps1 =
   3cb0000000000000

To see why, you might check out the help for eps. eps/4, the increment, is smaller than the distance between x and the next closest floating point number.

Notation and Defaults

Let's assume we have the start value in variable a, the end value in b, and the increment specified in d. If the user doesn't supply a step size d, then d is set to the default value of 1.

Handling Exceptions

The next part of the code weeds out situations in which the user supplies a non-finite value to : such as |NaN or Inf. In this case, MATLAB either returns NaN or an error with message identifier MATLAB:pmaxsize.

Note: There is an exception to this. Does anyone care to chime in with their guess about it? If so, please post here.

Check for Arguments Producing Empty Output

Next there are checks in the code looking for cases where the start value is larger than the end value and the increment is negative, and the complete opposite. In other words, a < b && d < 0 or a > b & d > 0. These cases and the case where d == 0 all produce empty vectors, in other words zeros(0,1).

The Main Event

Finally, here's where the real guts of the calculation take place.

The last two steps ensure that the vector output is symmetric about the mid-point of the vector's interval.

Want to See the Code?

Here it is now.

type colonop

function v = colonop(a,d,b)
% COLONOP  Demonstrate how the built-in a:d:b is constructed.
%
%   v = colonop(a,b) constructs v = a:1:b.
%   v = colonop(a,d,b) constructs v = a:d:b.
%
%   v = a:d:b is not constructed using repeated addition.  If the
%   textual representation of d in the source code cannot be
%   exactly represented in binary floating point, then repeated
%   addition will appear to have accumlated roundoff error.  In
%   some cases, d may be so small that the floating point number
%   nearest a+d is actually a.  Here are two important examples.
%
%   v = 1-eps : eps/4 : 1+eps is the nine floating point numbers
%   closest to v = 1 + (-4:1:4)*eps/4.  Since the spacing of the
%   floating point numbers between 1-eps and 1 is eps/2 and the
%   spacing between 1 and 1+eps is eps,
%   v = [1-eps 1-eps 1-eps/2 1 1 1 1 1+eps 1+eps].
%
%   Even though 0.01 is not exactly represented in binary,
%   v = -1 : 0.01 : 1 consists of 201 floating points numbers
%   centered symmetrically about zero.
%
%   Ideally, in exact arithmetic, for b > a and d > 0,
%   v = a:d:b should be the vector of length n+1 generated by
%   v = a + (0:n)*d where n = floor((b-a)/d).
%   In floating point arithmetic, the delicate computations
%   are the value of n, the value of the right hand end point,
%   c = a+n*d, and symmetry about the mid-point.

if nargin < 3
    b = d;
    d = 1;
end

tol = 2.0*eps*max(abs(a),abs(b));
sig = sign(d);

% Exceptional cases.

if ~isfinite(a) || ~isfinite(d) || ~isfinite(b)
   v = NaN;
   return
elseif d == 0 || a < b && d < 0 || b < a && d > 0
   % Result is empty.
   v = zeros(1,0);
   return
end

% n = number of intervals = length(v) - 1.

if a == floor(a) && d == 1
   % Consecutive integers.
   n = floor(b) - a;
elseif a == floor(a) && d == floor(d)
   % Integers with spacing > 1.
   q = floor(a/d);
   r = a - q*d;
   n = floor((b-r)/d) - q;
else
   % General case.
   n =  round((b-a)/d);
   if sig*(a+n*d - b) > tol
      n = n - 1;
   end
end

% c = right hand end point.

c = a + n*d;
if sig*(c-b) > -tol
   c = b;
end

% v should be symmetric about the mid-point.

v = zeros(1,n+1);
k = 0:floor(n/2);
v(1+k) = a + k*d;
v(n+1-k) = c - k*d;
if mod(n,2) == 0
   v(n/2+1) = (a+c)/2;
end

An Example

Let's try that example from above now.

format short
a = 1-eps;
b = 1+eps;
d = eps/4;
v = colonop(a,d,b)

v =
    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000

Now compare this to MATLAB's answer using :.

vc = a:d:b
isequal(v,vc)

vc =
    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000
ans =
     1

Colon is Numerically Delicate and Sophisticated

There's clearly more to : than meets the eye. If you have thoughts on my challenge from above, please post here.

Get the MATLAB code

Published with MATLAB® 7.4

Contents

Make Data

Create a "large" dataset. This one works for me with 1GB RAM and 1.5GB of swap space running Windows XP and no other applications.

n = 38*2^20;
x = randn(n,1);

Now call another function that will call either a regular function or an in-place function.

inplaceTest(x)

Let's look at the code for inplaceTest.

type inplaceTest