Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Vectorize silly matrix

Subject: Vectorize silly matrix

From: yoshis88@gmail.com

Date: 20 Apr, 2009 01:17:51

Message: 1 of 11

Dear comp.soft-sys.matlab,

First post here!

Anyway, I have a chunk of code that has been bothering me. I'm
somewhat of a veteran MATLAB/Octave user, and anytime I see a for
loop, the coder in me wants to vectorize it. So, here's something
which I have been picking on for quite a while.

function A = createSillyMatrix(r,c)
for k = 1:r
    for m = 1:c
        A(k,m) = k+m-1
    end
end

Any thoughts?

Subject: Vectorize silly matrix

From: Roger Stafford

Date: 20 Apr, 2009 02:54:01

Message: 2 of 11

yoshis88@gmail.com wrote in message <2b9d378f-1af9-40c1-81dc-1ad16f709c85@c12g2000yqc.googlegroups.com>...
> Dear comp.soft-sys.matlab,
>
> First post here!
>
> Anyway, I have a chunk of code that has been bothering me. I'm
> somewhat of a veteran MATLAB/Octave user, and anytime I see a for
> loop, the coder in me wants to vectorize it. So, here's something
> which I have been picking on for quite a while.
>
> function A = createSillyMatrix(r,c)
> for k = 1:r
> for m = 1:c
> A(k,m) = k+m-1
> end
> end
>
> Any thoughts?

 A = hankel(1:r,r:r+c-1);

Roger Stafford

Subject: Vectorize silly matrix

From: Jos

Date: 20 Apr, 2009 08:44:02

Message: 3 of 11

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gsgo49$ejh$1@fred.mathworks.com>...
> yoshis88@gmail.com wrote in message <2b9d378f-1af9-40c1-81dc-1ad16f709c85@c12g2000yqc.googlegroups.com>...
> > Dear comp.soft-sys.matlab,
> >
> > First post here!
> >
> > Anyway, I have a chunk of code that has been bothering me. I'm
> > somewhat of a veteran MATLAB/Octave user, and anytime I see a for
> > loop, the coder in me wants to vectorize it. So, here's something
> > which I have been picking on for quite a while.
> >
> > function A = createSillyMatrix(r,c)
> > for k = 1:r
> > for m = 1:c
> > A(k,m) = k+m-1
> > end
> > end
> >
> > Any thoughts?
>
> A = hankel(1:r,r:r+c-1);
>
> Roger Stafford

Some other matlab tricks:

cumsum([[1:r].' ones(r,c-1)],2)
bsxfun(@plus,[1:r].',1:c)-1

Jos

Subject: Vectorize silly matrix

From: T

Date: 20 Apr, 2009 10:04:01

Message: 4 of 11

>> r = 500;c=100;
[t(1)] = timeit(@() method1(r,c));
[t(2)] = timeit(@() method2(r,c));
[t(3)] = timeit(@() method3(r,c));
[t(4)] = timeit(@() method4(r,c));
t

t =

    0.1651 0.0023 0.0049 0.0028

where:
Method1
function A = method1(r,c)
for k = 1:r
    for m = 1:c
        A(k,m) = k+m-1;
    end
end

Method2
function A = method2(r,c)
[x,y]=meshgrid(1:c,1:r);
A = x+y-1;

Method3
function A = method3(r,c)
A = hankel(1:r,r:r+c-1);

Method4
function A = method4(r,c)
cumsum([[1:r].' ones(r,c-1)],2);
A = bsxfun(@plus,[1:r].',1:c)-1;

Subject: Vectorize silly matrix

From: Bruno Luong

Date: 20 Apr, 2009 10:27:01

Message: 5 of 11

"T " <tDO.Tdamsma@studentPO.INTtudelftPERI.ODnl> wrote in message <gshhah$jgp$1@fred.mathworks.com>...

>
> Method4
> function A = method4(r,c)
> cumsum([[1:r].' ones(r,c-1)],2);
> A = bsxfun(@plus,[1:r].',1:c)-1;

Jos gave us two different methods, not one.

Bruno

Subject: Vectorize silly matrix

From: T

Date: 20 Apr, 2009 10:48:01

Message: 6 of 11

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gshill$crv$1@fred.mathworks.com>...
> "T " <tDO.Tdamsma@studentPO.INTtudelftPERI.ODnl> wrote in message <gshhah$jgp$1@fred.mathworks.com>...
>
> >
> > Method4
> > function A = method4(r,c)
> > cumsum([[1:r].' ones(r,c-1)],2);
> > A = bsxfun(@plus,[1:r].',1:c)-1;
>
> Jos gave us two different methods, not one.
>
> Bruno

oops, it must be Monday morning...

>> r = 500;c=1000;
[t(1)] = timeit(@() method1(r,c));
[t(2)] = timeit(@() method2(r,c));
[t(3)] = timeit(@() method3(r,c));
[t(4)] = timeit(@() method4(r,c));
[t(5)] = timeit(@() method5(r,c));
t

t =

    3.6523 0.0337 0.0637 0.0494 0.0177
with

function A = method4(r,c)
A = cumsum([[1:r].' ones(r,c-1)],2);

function A = method5(r,c)
A = bsxfun(@plus,[1:r].',1:c)-1;

Subject: Vectorize silly matrix

From: T

Date: 20 Apr, 2009 10:51:01

Message: 7 of 11

and for small matrices:
>> r = 5;c=10;
t =

  1.0e-004 *

    0.1779 0.3368 0.6615 0.1212 0.3166

Subject: Vectorize silly matrix

From: T

Date: 20 Apr, 2009 12:04:01

Message: 8 of 11

I came up with some other options, and I've plotted the performance as function of r
d c for square matrices (r=c), see: http://www.damsma.net/performance.png

the fastest method really depends on the matrix size. Interesting that method 7, with a loop, is pretty fast. Also notice the strange behaviour at i=8 (r,c=512). Anyone have more suggestions for other methods?

1
for k = 1:r
    for m = 1:c
        A(k,m) = k+m-1;
    end
end

2
[x,y]=meshgrid(1:c,1:r);
A = x+y-1;

3
A = hankel(1:r,r:r+c-1);

4
A = cumsum([[1:r].' ones(r,c-1)],2);

5
A = bsxfun(@plus,[1:r].',1:c)-1;

6
A = repmat(1:c,r,1)+repmat((1:r)',1,c)-1;

7
A = ones(r,c);
for i=1:r
    A(i,:)=i:(c+i-1);
end

I tested for i=1:12, c,r = i^3, all using timeit, using matlab 2007b on a 2.41 ghz athlon 64 on WinXP with 1gb ram.

Subject: Vectorize silly matrix

From: Bruno Luong

Date: 20 Apr, 2009 12:23:01

Message: 9 of 11


>
> the fastest method really depends on the matrix size. Interesting that method 7, with a loop, is pretty fast. Also notice the strange behaviour at i=8 (r,c=512). Anyone have more suggestions for other methods?

Silly matrix indeed. I don't believe it's worth to find a sophisticated method. But here we go

5bis
A = bsxfun(@plus,(1:r).',0:c-1);

8
s=ones(1,r*c); s(1+r:r:end)=-r+2;
A = reshape(cumsum(s),r,c);

There is also a fliplr with spdiags, but I shall not bother to code it.

Bruno

Subject: Vectorize silly matrix

From: yoshis88@gmail.com

Date: 20 Apr, 2009 17:50:10

Message: 10 of 11

On Apr 20, 8:04=A0am, "T " <tDO.Tdam...@studentPO.INTtudelftPERI.ODnl>
wrote:
> I came up with some other options, and I've plotted the performance as fu=
nction of r
> d c for square matrices (r=3Dc), see:http://www.damsma.net/performance.pn=
g
>
> the fastest method really depends on the matrix size. Interesting that me=
thod 7, with a loop, is pretty fast. Also notice the strange behaviour at i=
=3D8 (r,c=3D512). Anyone have more suggestions for other methods?
>
> 1
> for k =3D 1:r
> =A0 =A0 for m =3D 1:c
> =A0 =A0 =A0 =A0 A(k,m) =3D k+m-1;
> =A0 =A0 end
> end
>
> 2
> [x,y]=3Dmeshgrid(1:c,1:r);
> A =3D x+y-1;
>
> 3
> A =3D hankel(1:r,r:r+c-1);
>
> 4
> A =3D cumsum([[1:r].' ones(r,c-1)],2);
>
> 5
> A =3D bsxfun(@plus,[1:r].',1:c)-1;
>
> 6
> A =3D repmat(1:c,r,1)+repmat((1:r)',1,c)-1;
>
> 7
> A =3D ones(r,c);
> for i=3D1:r
> =A0 =A0 A(i,:)=3Di:(c+i-1);
> end
>
> I tested for i=3D1:12, c,r =3D i^3, all using timeit, using matlab 2007b =
on a 2.41 ghz athlon 64 on WinXP with 1gb ram.

Wow! And I thought I was good at vectorization XD.
I was actually trying to come up with the cumsum solution on my own,
but had some trouble getting the right behavior.
I figured I might be able to get it with bsxfun, but I always had
trouble visualizing that function working.
I tried getting the repmat solution too, but didn't quite get there.
I didn't even think about meshgrid!

That really is quite curious that method 7 is so efficient. *ponders*
Thank you for the in-depth analyses on method vs time comparisons! It
has certainly opened my eyes.

-- Prashant

Subject: Vectorize silly matrix

From: T

Date: 21 Apr, 2009 07:46:02

Message: 11 of 11

Very true this is a silly problem. But it has a lot of solutions. I count 10 already. I made a new run: http://www.damsma.net/performance.png

I tested now for i=1:17, c,r = i^3, all using timeit, using matlab 2007b on a 2.41 ghz athlog 64 on WinXP with 1gb ram. At i=18, my pc ran out of memory (was already paging a lot). i=18 gives a matrix of 18^6, something like 34M entries. It actually calculated for method 2, (skipped method 1) but then it gave errors.

1
for k = 1:r
    for m = 1:c
        A(k,m) = k+m-1;
    end
end

2
[x,y]=meshgrid(1:c,1:r);
A = x+y-1;

3
A = hankel(1:r,r:r+c-1);

4
A = cumsum([[1:r].' ones(r,c-1)],2);

5
A = bsxfun(@plus,[1:r].',1:c)-1;

5bis
A = bsxfun(@plus,(1:r).',0:c-1);

6
A = repmat(1:c,r,1)+repmat((1:r)',1,c)-1;

7
A = ones(r,c);
for i=1:r
    A(i,:)=i:(c+i-1);
end

7bis
A = ones(r,c);
k=0:c-1;
for i=1:r
    A(i,:)=k+i;
end

8
s=ones(1,r*c); s(1+r:r:end)=-r+2;
A = reshape(cumsum(s),r,c);

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us