The engineer who attempted teaching me MatLab used a short script to demonstrate its matrix manipulation prowess. MatLab (shortened from Matrix Laboratory) is built from the ground up to do matrix mathematics. I wanted to translate his script to VBA, to see, 1. If it’s possible, and 2. How Excel compares. This is the original MatLab script:

%

% This is an exercise with MatLab Graphics and Flow Control. The results

% should be a graphic drawn by an Iterated Function System.

%

% It takes about 2-3 minutes to complete the drawing

%

figure % Establish a graphic window

pause(5.); % An opportunity to re-size the Figure window

hold on % Points accumulate in the window as an animation

x=[-5 0 5]; % 3 Initial points size the output axes

y=[0 1 10];

plot(x,y,’.g’,’MarkerSize’,4)

z = [x(2);y(2)];

%

pause(3.); % Pause before drawing

%

% 4 Matrices are key to the IFS Fractal description

%

fernM1 = [0 0

0 0.16];

fernM2 = [0.2 -0.26

0.23 0.22];

fernM3 = [-0.15 0.28

0.26 0.24];

fernM4 = [0.85 0.04

-0.04 0.85];

%

% Here is the fern

% We use tic/toc timing to measure execution time

%

tic

for kk = 1:20000

dd = rand();

if dd LTE 0.025

z = fernM1 * z;

plot(z(1),z(2),’.g’,’MarkerSize’,4)

pause(0.01); % Pause values less than 10ms makes no difference

elseif (dd GT 0.025) & (dd LTE 0.125)

z = fernM2 * z + [0 ; 1.6];

plot(z(1),z(2),’.g’,’MarkerSize’,4)

elseif (dd GT 0.125) & (dd LTE 0.225)

z = fernM3 * z + [0 ; 0.44];

plot(z(1),z(2),’.g’,’MarkerSize’,4)

else

z = fernM4 * z + [0 ; 1.6];

plot(z(1),z(2),’.g’,’MarkerSize’,4)

end

end

display(toc)

tic

for kk = 1:50000

dd = rand();

if dd LTE 0.025

z = fernM1 * z;

plot(z(1),z(2),’.g’,’MarkerSize’,4)

elseif (dd GT 0.025) & (dd LTE 0.125)

z = fernM2 * z + [0 ; 1.6];

plot(z(1),z(2),’.g’,’MarkerSize’,4)

elseif (dd GT 0.125) & (dd LTE 0.225)

z = fernM3 * z + [0 ; 0.44];

plot(z(1),z(2),’.g’,’MarkerSize’,4)

else

z = fernM4 * z + [0 ; 1.6];

plot(z(1),z(2),’.g’,’MarkerSize’,4)

end

end

display(toc)

display(‘The fern is done’)

Above, GT is “Greater Than”, LTE is “Less Than or Equals” and & is an ampersand. Here is my translation.

Dim x(1 To 1, 1 To 3) As Double

Dim y(1 To 1, 1 To 3) As Double

Dim z() As Variant

Dim FernM1(1 To 2, 1 To 2) As Double

Dim FernM2(1 To 2, 1 To 2) As Double

Dim FernM3(1 To 2, 1 To 2) As Double

Dim FernM4(1 To 2, 1 To 2) As Double

Dim kk As Long, dd As Double

Dim ZXArray(1 To 32000, 1 To 1) As Double

Dim ZYArray(1 To 32000, 1 To 1) As Double

Dim Fern As New Chart

ReDim z(1 To 2, 1 To 1) As Variant

x(1, 1) = -5: x(1, 2) = 0: x(1, 3) = 5

y(1, 1) = 0: y(1, 2) = 1: y(1, 3) = 10

z(1, 1) = x(1, 2)

z(2, 1) = y(1, 2)

FernM1(1, 1) = 0: FernM1(1, 2) = 0

FernM1(2, 1) = 0: FernM1(2, 2) = 0.16

FernM2(1, 1) = 0.2: FernM2(1, 2) = -0.26

FernM2(2, 1) = 0.23: FernM2(2, 2) = 0.22

FernM3(1, 1) = -0.15: FernM3(1, 2) = 0.28

FernM3(2, 1) = 0.26: FernM3(2, 2) = 0.24

FernM4(1, 1) = 0.85: FernM4(1, 2) = 0.04

FernM4(2, 1) = -0.04: FernM4(2, 2) = 0.85

For kk = 1 To 32000

dd = Rnd()

If dd LTE 0.025 Then

z = Application.WorksheetFunction.MMult(FernM1, z)

ElseIf dd GT 0.025 And dd LTE 0.125 Then

z = Application.WorksheetFunction.MMult(FernM2, z)

z(2, 1) = z(2, 1) + 1.6

Elseif dd GT 0.125 And dd LTE 0.225 Then

z = Application.WorksheetFunction.MMult(FernM3, z)

z(2, 1) = z(2, 1) + 0.44

Else

z = Application.WorksheetFunction.MMult(FernM4, z)

z(2, 1) = z(2, 1) + 1.6

End If

ZXArray(kk, 1) = z(1, 1)

ZYArray(kk, 1) = z(2, 1)

Next kk

Worksheets(“Sheet1”).Range(“A1:A32000”) = ZXArray

Worksheets(“Sheet1”).Range(“B1:B32000”) = ZYArray

Application.ScreenUpdating = False

Set Fern = Charts.Add

With Fern

.ChartType = xlXYScatter

.SetSourceData Source:=Sheets(“Sheet1”).Range(“A1:B32000”), PlotBy:=xlColumns

.Location Where:=xlLocationAsNewSheet, Name:=”Fern”

.Legend.Delete

With .SeriesCollection(1)

.MarkerBackgroundColorIndex = 10

.MarkerForegroundColorIndex = 10

.MarkerStyle = xlDiamond

.Smooth = False

.MarkerSize = 2

.Shadow = False

End With

With .PlotArea

.Border.LineStyle = xlNone

With .Interior

.ColorIndex = 2

.PatternColorIndex = 1

.Pattern = xlSolid

End With

End With

With .Axes(xlCategory)

.HasMajorGridlines = False

.HasMinorGridlines = False

.MajorTickMark = xlNone

.MinorTickMark = xlNone

.TickLabelPosition = xlNone

End With

With .Axes(xlValue)

.HasMajorGridlines = False

.HasMinorGridlines = False

.MajorTickMark = xlNone

.MinorTickMark = xlNone

.TickLabelPosition = xlNone

.MinimumScale = 0

.MaximumScale = 10

.MinorUnitIsAuto = True

.MajorUnit = 2

.Crosses = xlAutomatic

.ReversePlotOrder = False

.ScaleType = xlLinear

.DisplayUnit = xlNone

.Border.LineStyle = xlNone

End With

End With

Application.ScreenUpdating = True

End Sub

Again, GT is “Greater Than” and LTE is “Less Than or Equals”. Make the substitutions after you paste in the code to a new workbook. We have to do this or the WordPress parser treats them as HTML tags and doesn’t present the complete information, trying as it does to helpfully shorten everyone’s code.

Some observations: You don’t dimension variables in MatLab; and case matters, with x being a different variable from X. MatLab can plot at least 70,000 points, while Excel is limited to 32,000. The Excel limit is supposedly a series limit, but I could never plot two series totaling over 32K points. MatLab doesn’t need a spreadsheet to plot, so if there’s a way in Excel to plot values without going through a spreadsheet, would someone please so describe. I don’t think you can. As expected, the Sheet1 ranges are the same dimensions as the ZX and ZY arrays, which are at the series limits, and allows copying the arrays to the spreadsheet.

Matrix Multiplication in VBA requires dimensioning the VB arrays as if they were cell arrays so that the worksheet function MMULT() can be used. Z() is a variant re-dimed to the dimensions of the matrix multiplication, meeting the specifications MMULT() requires in a spreadsheet–two rows by one column. I can’t make the code run on a Mac, where we’re stuck at VBA5. The Mac won’t accept the assignment of the array.

And finally, while it’s less than half the points, it’s **much less** than half the time to run. It’s almost instantaneous on my machine. That may mean the MatLab script is not optimum. I’ll pass ~~on~~ along any observations MatLab veterans may make.

Give it a try. You’ll see why it was chosen as a demo. It’s worth the trouble. A picture is worth one kilo-word, and you’ll have it in a blink.

*…mrt*

So where is the extra time used in MATLAB’s implementation? Is it slower in general or, for example, is it faster at generating the points and slower at plotting them?

Anyone doing matrix maths and manipulation should really look at APL. Ken Iverson’s language has been around since the 60s; it is clear, concise and compact.

Cheers

Andy

For loops in Matlab are notoriously slow. You are supposed to take advantage of the various matrix and vector structures and avoid for loops as much as possible. So you would do something like rand([50000,1]) outside of a for loop to get dd to be a vector with 50000 numbers in it. Then you could use dd(i) to access each element in it within the for loop. I haven’t taken the time to figure out a way to get rid of everything else in the for loops, but taking the rand out of the loop may improve the speed a little.

Impressive

Very cool, Michael. BTW, I took the number of data points up to 128,000 in Excel 2010 32-bit and it worked great.

Tim –

Wow. That’s great. Did you know what you were going to get? ;-)

…mrtFirst off, that’s really cool. I used to try and do things like this comparing excel and vba, but since my vba abilities weren’t anywhere near as good as my matlab abilities there was no comparison.

Couple observations from my experience w/ Matlab:

1) Matlab handles arrays better than Excel in every possible regard – from creating them to doing any sort of analysis or calculation, it is easier and faster in Matlab. Case in point – scalar to multidimensional calculations can be done directly in Matlab whereas you need to loop through each element individually in Excel

2) You can dimension arrays in Matlab, which speeds things up quite a bit w/ large data sets

I don’t have Matlab on my computer any more so I can’t check it, but I’m guessing the reason there’s such a discrepancy between the two in terms of time is that the Matlab version is plotting each point as it’s calculated (matlab equivalent of leaving screenupdating on) vs. VBA does all the calculations then plots the result.

I agree with Bjacobowski. I took out the plotting routines and just saved the data as the iteration proceeded. Then I did the plot at the end. The second loop (50,000 iterations) went from about 28 seconds on my laptop to about 0.08 seconds. I didn’t try the VBA version.

If I do this with both loops and take out the pauses, the whole thing takes about 0.1 seconds.

Michael,

The reason that your MATLAB code is ‘slow’ is that every single iteration of the loop that computes the fern you are updating the plot. It is generally not a good idea to do this. I modified your code so that it simply computes the fern, and then plots it (see below). Running on even a fairly old tablet from IBM (X61) in R2010b the times I get for the 20000 iterations and then 50000 are

>> fern

ans =

0.1405

ans =

0.3483

The fern is done

code – apologies if the HTML gets mangled …

% fern.m

%

% This is an exercise with MatLab Graphics and Flow Control. The results

% should be a graphic drawn by an Iterated Function System.

%

% It takes about 2-3 minutes to complete the drawing

%

figure(1) % Establish a graphic window – use the same on each time

clf % Clear it and set the axis

axis([-5 5 0 10]);

hold on % Points accumulate in the window as an animation

% Initial values

z = [0;1];

%

% 4 Matrices are key to the IFS Fractal description

%

fernM1 = [0 0 ; 0 0.16];

fernM2 = [0.2 -0.26 ; 0.23 0.22];

fernM3 = [-0.15 0.28 ; 0.26 0.24];

fernM4 = [0.85 0.04 ; -0.04 0.85];

%

% Here is the fern

% We use tic/toc timing to measure execution time

%

% How many points are we making this time? Allocate the array

N = 20000;

allZ = zeros(2, N);

tic

for kk = 1:N

dd = rand();

if dd 0.025) && (dd 0.125) && (dd <= 0.225)

z = fernM3 * z + [0 ; 0.44];

else

z = fernM4 * z + [0 ; 1.6];

end

allZ(:, kk) = z;

end

display(toc)

% Plot what we have so far …

plot(allZ(1,:), allZ(2,:), ‘.g’,’MarkerSize’, 4)

pause(2); % See the sparse fern for 2 seconds

% Now do more points

N = 50000;

allZ = zeros(2, N);

tic

for kk = 1:N

dd = rand();

if dd 0.025) && (dd 0.125) && (dd <= 0.225)

z = fernM3 * z + [0 ; 0.44];

else

z = fernM4 * z + [0 ; 1.6];

end

allZ(:, kk) = z;

end

display(toc)

% plot the rest (hold on so

plot(allZ(1,:), allZ(2,:), ‘.g’,’MarkerSize’, 4)

display(‘The fern is done’)

The following comments have been received from Doug Jenkins over at Newton Excel Bach (see http://newtonexcelbach.wordpress.com/about/). Doug is a structural engineer and uses Excel, as he says, for functions without a dollar sign. He is a frequent commentator here, but at the moment is DDoE-challenged.

——-

Expect some flack from the Matlab users (but on the other hand no “real” Matlab user would be seen dead at a site called Daily Dose of Excel, so maybe not :))

There’s no doubt that for real matrix work Matlab will handle far bigger matrices, and deal with them far quicker, than the built in Excel functions, but if you are wanting to do that sort of thing in Excel there are tools available which will give you Matlab like functionality. For instance see:

http://newtonexcelbach.wordpress.com/2010/05/27/compiled-alglib-matrix-functions-for-excel/

I took the liberty of playing with your code to make the following changes:

Read the number of graph points, and the basic “fern” parameters from the spreadsheet.

Check if the chart sheet exists

Alternative code using variants and doubles (with timer for each!)

VBA matrix multiply instead of using worksheetfunction.

Interesting what small changes to the parameters do to the fern plot. Excel 2010 will handle up to 1

million data points by the way (although the fern looks nicer with about 30000!).

I’d appreciate it if you would post the comments at DDoE for me.

——-

I’m pleased to do so. Doug will be putting up his code over at his place. The fern is traveling across the Pacific.

…mrtJos –

Thank you, but to be clear, it’s not my Matlab code, I’m just the translator. I don’t have the imagination to conceive either the elegance of the math or the beauty of the result. The author is a PhD mathematician who shared his time with such as me. I’ll pass along the comments and perhaps he’ll join us.

…mrtThat’s some pretty cool stuff – good translating!

Mat

What about minitab?

I’ve heard from the author. He’s not the sort to blog, but he did give me this wiki pointer:

http://en.wikipedia.org/wiki/Barnsley%27s_fern

Michael Barnsley has several patents on fractal compression. At the above link are a good explanation of the math behind the fern, and two examples of “mutant” ferns from different coefficients.

…mrt> so if there’s a way in Excel to plot values without going through a spreadsheet, would someone please so describe.

Define names:

Names.Add “x”, ZXArray

Names.Add “y”, ZYArray

Create the chart based on a dummy series eg: [a1:b2], then update the series values

.XValues = “Sheet1!x”

.Values = “Sheet1!y”

Now you should be able to remove all worksheets and just leave the chart sheet.

Btw, nice work! Lori

Interesting post,

The biggest performance hog in the Matlab code is, as Jos points out, the updating of the plot.

With some help from the gurus over in the comp.soft-sys.matlab newsgroup, here you have an Matlab optimized version courtesy of Bruno Luong.

Using this code the average of 50 trials gives me 0.0893 seconds to calc 50000 points (this is on a Intel Q9550, 4Gb ram, Vista, Matlab R2009a).

Now that should surely knock the socks of any VBA routine :-P

(still haven’t got completely rid of the for loop though)

% Initial values

z = [0;1];

%

% 4 Matrices are key to the IFS Fractal description

%

fernM1 = [0 0 ; 0 0.16];

fernM2 = [0.2 -0.26 ; 0.23 0.22];

fernM3 = [-0.15 0.28 ; 0.26 0.24];

fernM4 = [0.85 0.04 ; -0.04 0.85];

N = 50000;

allZ = zeros(2, N);

dd = rand(1,N);

[na, where] = histc(dd, [0 0.025 0.125 0.225 1]);

fern = { fernM1,fernM2,fernM3,fernM4 };

a = cat(2, [0; 0], [0 ; 1.6], [0 ; 0.44], [0 ; 1.6]);

tic

for kk = 1:N

w = where(kk);

z = fern{w}*z + a(:,w);

allZ(:,kk) = z;

end

display(toc)

plot(allZ(1,:), allZ(2,:), ‘.g’,’MarkerSize’, 4)

disp(‘The fern is done’);

Peder –

Thank you, and thanks to Bruno also. I forwarded your script to the author. There are things in that script I know weren’t in the course.

Lori –

Thank you for the compliment, and thank you for the tip. That’s a nice trick.

…mrtThose are curly apostrophes in Peder’s script, compliments of WordPress. I can’t speak to Matlab, but I know Octave, the Gnu Matlab, doesn’t like them.

Hi,

VBA plots matrices directly seems tô MATLAB.