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:
Can my code be made more efficient using MATLAB's vectorization?

Subject: Can my code be made more efficient using MATLAB's vectorization?

From: Jeff

Date: 14 Feb, 2013 17:49:15

Message: 1 of 5

Is there a better, more MATLAB-y way, to write this code (that is, using vectors)? Eventually, this will be used for large amounts of data (from a program running on a cluster), so I need it as efficient as possible.

U and W are T rows x N columns matrices. Variable evals is a 1 by N matrix of the eigenvalues (of some other matrix which doesn't matter for this post). For now, the largest value of T is about 1000, but eventually it will be many thousands (maybe even millions).

So the first thing I need to do (variable P_p) is subtract the each entry on the row from the one to its right (for now, I only have periodic conditions coded). I think I have that variable coded fairly efficiently (yes? no?).

if strcmp(boundary,'periodic')
    nMinusOne=[N 1:N-1];
    nPlusOne=[2:N 1];
elseif strcmp(boundary,'fixed')
    nMinusOne=[1:N-1];
    nPlusOne=[2:N];
end

P_p = ((U(:,nPlusOne)-U(:,n)).^2);
Q_p = Udot(:,n).^2;
H = P_p + Q_p;


But for variable P_f I need to square each element of W, then multiply everything in column 1 by eval(1), everything in column 2 by evals(2), ... everything in column n by evals(n) (where, obviously, n=1..N). Finally, I need to divide each element of P_f by N. In math terms, I'm trying to do something like

(W_{i,j}^2)*eval_j
------------------
N

And this code is the most efficient I could think of. It's not terribly fast (on my quad core PC). Can it be improved?

P_f = zeros(t_end,N);
for m=1:t_end
    P_f = W(m,:).^2.*evals;
end
P_f = -P_f./N;

Q_f = (Wdot.^2)./n;
H_f = P_f + Q_f;

Subject: Can my code be made more efficient using MATLAB's vectorization?

From: Bruno Luong

Date: 14 Feb, 2013 20:26:14

Message: 2 of 5

"Jeff" wrote in message <kfj82r$k50$1@newscl01ah.mathworks.com>...

>
> P_f = zeros(t_end,N);
> for m=1:t_end
> P_f = W(m,:).^2.*evals;
> end
> P_f = -P_f./N;
>

Bad: you post an untested code (anyone can spot an obvious mistake), so I do the same:

P_f = bsxfun(@times, W.^2, evals/(-N));

% Bruno

Subject: Can my code be made more efficient using MATLAB's vectorization?

From: dpb

Date: 14 Feb, 2013 20:32:21

Message: 3 of 5

On 2/14/2013 11:49 AM, Jeff wrote:
> Is there a better, more MATLAB-y way, to write this code (that is,
> using vectors)? Eventually, this will be used for large amounts of
> data (from a program running on a cluster), so I need it as efficient
> as possible.
>
> U and W are T rows x N columns matrices. Variable evals is a 1 by N
> matrix of the eigenvalues (of some other matrix which doesn't matter
> for this post). For now, the largest value of T is about 1000, but
> eventually it will be many thousands (maybe even millions).
>
> So the first thing I need to do (variable P_p) is subtract the each
> entry on the row from the one to its right (for now, I only have
> periodic conditions coded). I think I have that variable coded fairly
> efficiently (yes? no?).
>
> if strcmp(boundary,'periodic')
> nMinusOne=[N 1:N-1];
> nPlusOne=[2:N 1];
> elseif strcmp(boundary,'fixed')
> nMinusOne=[1:N-1];
> nPlusOne=[2:N];
> end

...

First, a couple of questions...

a) The periodic solution as written will be N columns whereas the fixed
will be N-1. Is this intended/desired?

b) As written the difference between n+ - n- of the array for periodic
will be that of column n minus n-2 rather than the "one on its right"
which implies a first difference. To illustrate if N were 5, say, the
two column vectors you generate are

    [5 1 2 3 4]
    [2 3 4 5 1]

which results in deltas of 2-5, 3-1, 4-2, etc., ... Is this intended or
did you intend

    [5 1 2 3 4]
    [1 2 3 4 5]

instead?

W/ those answers, then can progress to next step (altho I'll note you
may find

doc diff
doc circshift

of interest)...

--

Subject: Can my code be made more efficient using MATLAB's vectorization?

From: Jeff

Date: 16 Feb, 2013 22:51:05

Message: 4 of 5

Hi. Slow response because I forgot to add the thread to my watchlist. But thanks for helping.

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message
>
> P_f = bsxfun(@times, W.^2, evals/(-N));
>

Thanks Bruno. Despite my sloppy code, I think you gave me what I need: bsxfun.


dpb <none@non.net> wrote in message <kfjhkn$8ea$1@speranza.aioe.org>...
> First, a couple of questions...
>
> a) The periodic solution as written will be N columns whereas the fixed
> will be N-1. Is this intended/desired?

dpb:
a) That is not what was intended. The fixed boundary code is not complete (hopefully this is the "untested code" Bruno was referring too) and should not have been included. Eventually, the fixed code will look something like

P_p = ([U(:,nPlusOne); 0]-U(:,n)).^2;

except that 0 will really have to be something like zeros(T,1). I'm working on the periodic code first and will come back to that.

dpb <none@non.net> wrote in message <kfjhkn$8ea$1@speranza.aioe.org>...
> b) As written the difference between n+ - n- of the array for periodic
> will be that of column n minus n-2 rather than the "one on its right"
> which implies a first difference. To illustrate if N were 5, say, the
> two column vectors you generate are
>

b) Are you talking about the code for P_p? I am subtracting the n columns from the n+1 columns (that is, u(n+1) minus u(n)). I think you read my code as u(n+1) minus u(n). I hope so, anyway. Because if you read my code right and you think it is doing u(n+1) minus u(n-1), then in trouble!

I think my code already does what circshift does, at least for the periodic case. diff seems to be more for taking derivatives, so I don't think it will help.


Sorry to all who read this for the sloppy code. I'm learning the math and matlab coding at the same time.

Subject: Can my code be made more efficient using MATLAB's vectorization?

From: dpb

Date: 16 Feb, 2013 23:55:53

Message: 5 of 5

On 2/16/2013 4:51 PM, Jeff wrote:
...

> dpb <none@non.net> wrote in message <kfjhkn$8ea$1@speranza.aioe.org>...
>> First, a couple of questions...
>>
>> a) The periodic solution as written will be N columns whereas the
>> fixed will be N-1. Is this intended/desired?
>
> dpb:
> a) That is not what was intended. The fixed boundary code is not
> complete (hopefully this is the "untested code" Bruno was referring too)
> and should not have been included. Eventually, the fixed code will look
> something like
>
> P_p = ([U(:,nPlusOne); 0]-U(:,n)).^2;
>
> except that 0 will really have to be something like zeros(T,1). I'm
> working on the periodic code first and will come back to that.
>
> dpb <none@non.net> wrote in message <kfjhkn$8ea$1@speranza.aioe.org>...
>> b) As written the difference between n+ - n- of the array for periodic
>> will be that of column n minus n-2 rather than the "one on its right"
>> which implies a first difference. To illustrate if N were 5, say, the
>> two column vectors you generate are
>>
>
> b) Are you talking about the code for P_p? I am subtracting the n
> columns from the n+1 columns (that is, u(n+1) minus u(n)). I think you
> read my code as u(n+1) minus u(n). I hope so, anyway. Because if you
> read my code right and you think it is doing u(n+1) minus u(n-1), then
> in trouble!
>
> I think my code already does what circshift does, at least for the
> periodic case. diff seems to be more for taking derivatives, so I don't
> think it will help.
...

Yes, but you didn't define n so I presumed there was a typo and since
you computed both a nPlusOne and a nMinusOne vector that the other
vector was intended to be nMinusOne. If it's just 1:n then you don't
need a subscript and :,n would reference one column (the nth one) unless
n=[1:size(U,2)]

I haven't compared for relative speed but U(:,nPlusOne)-U is identical
to circshift(U,[0 1])-U

Well, looking I see that circshift is an m-file so probably the
special-case precomputed indices will outperform the packaged more
general function.

diff() will take care of the subtraction as well unless you're sure you
do want the difference array to be the same size--I was just checking
that you think the wrapped position really does have meaning...if that's
what you want, go for it.

--

Tags for 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