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:
Vectorizing nested for loops with loop dependencies

Subject: Vectorizing nested for loops with loop dependencies

From: Michael

Date: 3 Jun, 2013 11:21:09

Message: 1 of 15

I have a function that takes a nx1 time series "ts" as input and returns a nxn matrix "V".
Now imagine "ts" as a landscape: Whenever point b can be seen from point a (i.e.: a and b can be connected with a straight line), V(a,b) ~= 0. If point b cannot be seen from point a, V(a,b) = 0.
Here's the function:

function V = someFunction(ts, weighted)

len = length(ts);
V = zeros(len,len);

for a = 1 : len,
   for b = 1 : a-1,
       if a-b == 1, %if points a and b are neighbours.
           if weighted == 1,
               V(a,b) = (ts(a) - ts(b)) / (a - b); %Calculate the slope between a and b.
           elseif weighted == 0,
               V(a,b) = 1;
           end
       elseif a-b > 1, %if points a and b are separated by at least one other point.
           for c = b+1 : a-1,
                maxHeight = ts(b) + (ts(a) - ts(b)) * ((b - c) / (b - a));
                if ts(c) > maxHeight,
                    break; % Stop the inner loop the moment an intersection was found.
                elseif c == a-1, %if it's the last iteration of c.
                    if weighted == 1,
                        V(a,b) = (ts(a) - ts(b)) / (a - b); %Calculate slope between a and b.
                    elseif weighted == 0,
                      V(a,b) = 1;
                    end
                end
           end
       end
   end
end

end

I need to vectorize the three nested loops, but I'm having a hard time doing just that (especially the innermost loop gives a lot of trouble).
Any help would be appreciated...

Subject: Vectorizing nested for loops with loop dependencies

From: dpb

Date: 3 Jun, 2013 14:03:08

Message: 2 of 15

On 6/3/2013 6:21 AM, Michael wrote:
> I have a function that takes a nx1 time series "ts" as input and returns
> a nxn matrix "V".
> Now imagine "ts" as a landscape: Whenever point b can be seen from point
> a (i.e.: a and b can be connected with a straight line), V(a,b) ~= 0. If
> point b cannot be seen from point a, V(a,b) = 0.
> Here's the function:
>
> function V = someFunction(ts, weighted)
> L = length(ts);
> V = zeros(L,L);
> for a=1:L
> for b=1:a-1,
> if a-b==1
> if weighted == 1,
> V(a,b) = (ts(a)-ts(b))/(a - b);
> elseif weighted == 0,
> V(a,b) = 1;
> end
> elseif a-b > 1
> for c = b+1:a-1
> maxHeight = ts(b) + (ts(a)-ts(b))*((b-c)/(b-a));
> if ts(c) > maxHeight,
> break; % Stop the inner loop the moment an intersection was found.
> elseif c==a-1, %if it's the last iteration of c.
> if weighted == 1,
> V(a,b) = (ts(a)-ts(b))/(a - b);
> elseif weighted == 0,
> V(a,b) = 1;
> end
             end
> end
> end
> end
> end
> end
>
> I need to vectorize the three nested loops, but I'm having a hard time
> doing just that (especially the innermost loop gives a lot of trouble).
> Any help would be appreciated...

Well, let's just start and see what we're doing...ok, the first portion of

if a-b==1
   if weighted == 1,
     V(a,b) = (ts(a)-ts(b))/(a - b);
   elseif weighted == 0,
     V(a,b) = 1;
   end
...

is simply

if W
   V=diag(diff(ts),-1);
else
   V=diag(ones(1:L-1),-1);
end

and it's outside the looping entirely.

Now

elseif a-b > 1
   for c = b+1:a-1
     maxHeight = ts(b) + (ts(a)-ts(b))*((b-c)/(b-a));
     if ts(c) > maxHeight,
       break; % Stop the inner loop intersection found.
     elseif c==a-1, %if it's the last iteration of c.
       if weighted == 1,
         V(a,b) = (ts(a)-ts(b))/(a - b);
       elseif weighted == 0,
         V(a,b) = 1;
       end
     end
   end
end


I guess I don't in first blush follow what you're actually after...if I
work through the indexing of a small (say 4x4) array I find that for the
lower triangular of lower than the first off-diagonal the following
indices and the associated range for the c loop

a b c
---- ----
3 1 2:2
4 1 2:3
4 2 3:3

But the following elseif of c==a-1 means that only the following indices
in V can be stored into

a b c c==a-1 a b
---- --- ------ ---
3 1 2:2 2 3,2
4 1 2:3 3 4,3
4 2 3:3 3 4,3

So, looks to me like you're doing a whole lot of looping to reset the
same element over and over but it's still only addressing the first
lower diagonal in the end.

In which case it seems like you could compute a logical vector on the
maximum distance condition you want and simply do a single diag() using
that.

I didn't try to work this out on the machine; just tried to read the
coding fairly quickly so I may have made an indexing boo-boo...

--

Subject: Vectorizing nested for loops with loop dependencies

From: Michael

Date: 3 Jun, 2013 15:08:08

Message: 3 of 15

dpb <none@non.net> wrote in message <koi7mg$nrc$1@speranza.aioe.org>...
> On 6/3/2013 6:21 AM, Michael wrote:
> > I have a function that takes a nx1 time series "ts" as input and returns
> > a nxn matrix "V".
> > Now imagine "ts" as a landscape: Whenever point b can be seen from point
> > a (i.e.: a and b can be connected with a straight line), V(a,b) ~= 0. If
> > point b cannot be seen from point a, V(a,b) = 0.
> > Here's the function:
> >
> > function V = someFunction(ts, weighted)
> > L = length(ts);
> > V = zeros(L,L);
> > for a=1:L
> > for b=1:a-1,
> > if a-b==1
> > if weighted == 1,
> > V(a,b) = (ts(a)-ts(b))/(a - b);
> > elseif weighted == 0,
> > V(a,b) = 1;
> > end
> > elseif a-b > 1
> > for c = b+1:a-1
> > maxHeight = ts(b) + (ts(a)-ts(b))*((b-c)/(b-a));
> > if ts(c) > maxHeight,
> > break; % Stop the inner loop the moment an intersection was found.
> > elseif c==a-1, %if it's the last iteration of c.
> > if weighted == 1,
> > V(a,b) = (ts(a)-ts(b))/(a - b);
> > elseif weighted == 0,
> > V(a,b) = 1;
> > end
> end
> > end
> > end
> > end
> > end
> > end
> >
> > I need to vectorize the three nested loops, but I'm having a hard time
> > doing just that (especially the innermost loop gives a lot of trouble).
> > Any help would be appreciated...
>
> Well, let's just start and see what we're doing...ok, the first portion of
>
> if a-b==1
> if weighted == 1,
> V(a,b) = (ts(a)-ts(b))/(a - b);
> elseif weighted == 0,
> V(a,b) = 1;
> end
> ...
>
> is simply
>
> if W
> V=diag(diff(ts),-1);
> else
> V=diag(ones(1:L-1),-1);
> end
>
> and it's outside the looping entirely.
>
> Now
>
> elseif a-b > 1
> for c = b+1:a-1
> maxHeight = ts(b) + (ts(a)-ts(b))*((b-c)/(b-a));
> if ts(c) > maxHeight,
> break; % Stop the inner loop intersection found.
> elseif c==a-1, %if it's the last iteration of c.
> if weighted == 1,
> V(a,b) = (ts(a)-ts(b))/(a - b);
> elseif weighted == 0,
> V(a,b) = 1;
> end
> end
> end
> end
>
>
> I guess I don't in first blush follow what you're actually after...if I
> work through the indexing of a small (say 4x4) array I find that for the
> lower triangular of lower than the first off-diagonal the following
> indices and the associated range for the c loop
>
> a b c
> ---- ----
> 3 1 2:2
> 4 1 2:3
> 4 2 3:3
>
> But the following elseif of c==a-1 means that only the following indices
> in V can be stored into
>
> a b c c==a-1 a b
> ---- --- ------ ---
> 3 1 2:2 2 3,2
> 4 1 2:3 3 4,3
> 4 2 3:3 3 4,3
>
> So, looks to me like you're doing a whole lot of looping to reset the
> same element over and over but it's still only addressing the first
> lower diagonal in the end.
>
> In which case it seems like you could compute a logical vector on the
> maximum distance condition you want and simply do a single diag() using
> that.
>
> I didn't try to work this out on the machine; just tried to read the
> coding fairly quickly so I may have made an indexing boo-boo...
>
> --

dpb,

your solution for the upper part of the code (using diag and diff) is brilliant - I implemented it right away.
What I'm trying to achieve in the lower part is the following:
1) Select two non-neighbouring points (a and b) from ts.
2) Check if there are points between a and b that would intersect a straight line from a to b. If there are no intersecting points, then matrix element V(a,b) is assigned a nonzero value.
If the option "weighted" equals 0, the value should be 1.
If the option "weighted" equals 1, the value should be the slope of the line connecting a and b.
3) Do the previous steps for all possible combinations of points a and b.

I hope this clarifies what the lower part code is supposed to do...

Subject: Vectorizing nested for loops with loop dependencies

From: Michael

Date: 3 Jun, 2013 17:45:11

Message: 4 of 15

dpb <none@non.net> wrote in message <koi7mg$nrc$1@speranza.aioe.org>...
> On 6/3/2013 6:21 AM, Michael wrote:
> > I have a function that takes a nx1 time series "ts" as input and returns
> > a nxn matrix "V".
> > Now imagine "ts" as a landscape: Whenever point b can be seen from point
> > a (i.e.: a and b can be connected with a straight line), V(a,b) ~= 0. If
> > point b cannot be seen from point a, V(a,b) = 0.
> > Here's the function:
> >
> > function V = someFunction(ts, weighted)
> > L = length(ts);
> > V = zeros(L,L);
> > for a=1:L
> > for b=1:a-1,
> > if a-b==1
> > if weighted == 1,
> > V(a,b) = (ts(a)-ts(b))/(a - b);
> > elseif weighted == 0,
> > V(a,b) = 1;
> > end
> > elseif a-b > 1
> > for c = b+1:a-1
> > maxHeight = ts(b) + (ts(a)-ts(b))*((b-c)/(b-a));
> > if ts(c) > maxHeight,
> > break; % Stop the inner loop the moment an intersection was found.
> > elseif c==a-1, %if it's the last iteration of c.
> > if weighted == 1,
> > V(a,b) = (ts(a)-ts(b))/(a - b);
> > elseif weighted == 0,
> > V(a,b) = 1;
> > end
> end
> > end
> > end
> > end
> > end
> > end
> >
> > I need to vectorize the three nested loops, but I'm having a hard time
> > doing just that (especially the innermost loop gives a lot of trouble).
> > Any help would be appreciated...
>
> Well, let's just start and see what we're doing...ok, the first portion of
>
> if a-b==1
> if weighted == 1,
> V(a,b) = (ts(a)-ts(b))/(a - b);
> elseif weighted == 0,
> V(a,b) = 1;
> end
> ...
>
> is simply
>
> if W
> V=diag(diff(ts),-1);
> else
> V=diag(ones(1:L-1),-1);
> end
>
> and it's outside the looping entirely.
>
> Now
>
> elseif a-b > 1
> for c = b+1:a-1
> maxHeight = ts(b) + (ts(a)-ts(b))*((b-c)/(b-a));
> if ts(c) > maxHeight,
> break; % Stop the inner loop intersection found.
> elseif c==a-1, %if it's the last iteration of c.
> if weighted == 1,
> V(a,b) = (ts(a)-ts(b))/(a - b);
> elseif weighted == 0,
> V(a,b) = 1;
> end
> end
> end
> end
>
>
> I guess I don't in first blush follow what you're actually after...if I
> work through the indexing of a small (say 4x4) array I find that for the
> lower triangular of lower than the first off-diagonal the following
> indices and the associated range for the c loop
>
> a b c
> ---- ----
> 3 1 2:2
> 4 1 2:3
> 4 2 3:3
>
> But the following elseif of c==a-1 means that only the following indices
> in V can be stored into
>
> a b c c==a-1 a b
> ---- --- ------ ---
> 3 1 2:2 2 3,2
> 4 1 2:3 3 4,3
> 4 2 3:3 3 4,3
>
> So, looks to me like you're doing a whole lot of looping to reset the
> same element over and over but it's still only addressing the first
> lower diagonal in the end.
>
> In which case it seems like you could compute a logical vector on the
> maximum distance condition you want and simply do a single diag() using
> that.
>
> I didn't try to work this out on the machine; just tried to read the
> coding fairly quickly so I may have made an indexing boo-boo...
>
> --

@dpb:
I've just re-read your reply: If I understand you correctly, you're basically saying that my loops only fill the first lower diagonal of the matrix. Or am I mistaken?
Anyway: If you run the code, you will see that this is clearly not the case - there are also other elements of the matrix (beside the lower diagonal) that are assigned a value.

Subject: Vectorizing nested for loops with loop dependencies

From: dpb

Date: 3 Jun, 2013 18:18:01

Message: 5 of 15

On 6/3/2013 10:08 AM, Michael wrote:
> dpb <none@non.net> wrote in message <koi7mg$nrc$1@speranza.aioe.org>...
...

>> Well, let's just start and see what we're doing...ok, the first
>> portion of
>>
...[snipped for brevity]...
>> ...
>>
>> is simply
>>
>> if W
>> V=diag(diff(ts),-1);
>> else
>> V=diag(ones(1:L-1),-1);
>> end
>>
>> and it's outside the looping entirely.
>>
>> Now
>>
>> elseif a-b > 1
>> for c = b+1:a-1
>> maxHeight = ts(b) + (ts(a)-ts(b))*((b-c)/(b-a));
>> if ts(c) > maxHeight,
>> break; % Stop the inner loop intersection found.
>> elseif c==a-1, %if it's the last iteration of c.
>> if weighted == 1,
>> V(a,b) = (ts(a)-ts(b))/(a - b);
>> elseif weighted == 0,
>> V(a,b) = 1;
>> end
>> end
>> end
>> end
>>
>>
>> I guess I don't in first blush follow what you're actually after...if
>> I work through the indexing of a small (say 4x4) array I find that for
>> the lower triangular of lower than the first off-diagonal the
>> following indices and the associated range for the c loop
>>
>> a b c
>> ---- ----
>> 3 1 2:2
>> 4 1 2:3
>> 4 2 3:3
>>
>> But the following elseif of c==a-1 means that only the following
>> indices in V can be stored into
>>
>> a b c c==a-1 a b
>> ---- --- ------ ---
>> 3 1 2:2 2 3,2
>> 4 1 2:3 3 4,3
>> 4 2 3:3 3 4,3
>>
>> So, looks to me like you're doing a whole lot of looping to reset the
>> same element over and over but it's still only addressing the first
>> lower diagonal in the end.
>>
>> In which case it seems like you could compute a logical vector on the
>> maximum distance condition you want and simply do a single diag()
>> using that.
>>
>> I didn't try to work this out on the machine; just tried to read the
>> coding fairly quickly so I may have made an indexing boo-boo...
>>
...

> dpb,
>
> your solution for the upper part of the code (using diag and diff) is
> brilliant - I implemented it right away.
> What I'm trying to achieve in the lower part is the following:
> 1) Select two non-neighbouring points (a and b) from ts.
> 2) Check if there are points between a and b that would intersect a
> straight line from a to b. If there are no intersecting points, then
> matrix element V(a,b) is assigned a nonzero value. If the option
> "weighted" equals 0, the value should be 1. If the option "weighted"
> equals 1, the value should be the slope of the line connecting a and b.
> 3) Do the previous steps for all possible combinations of points a and b.
>
> I hope this clarifies what the lower part code is supposed to do...

OK, I see on second reading I did goof--the a,b don't change on the c
loop, obviously, so the LT are filled...my bad there.

I don't have time to cogitate over this more at the moment but hopefully
later on will be able to look at it some more, maybe.

I think I'd look at how to precompute the location of the intersection
and then store...mayhaps interp1() properly applied will do something
useful; I'll have to think about that.

Meanwhile maybe somebody else will stumble along and have an "aha!"
moment...

Good luck...

--

Subject: Vectorizing nested for loops with loop dependencies

From: Bruno Luong

Date: 3 Jun, 2013 19:20:12

Message: 6 of 15

What is typical value of len?

Here is an equivalent function for weight > 0.

function V = someFunction(ts)
len = length(ts);
i = 1:len;
V = bsxfun(@minus,ts(:),ts(:)') ./ bsxfun(@minus,i(:),i(:)');
V = tril(V,-1);
for a = 1 : len,
    for b = 1 : a-1,
        i = b+1:a-1;
        if any(ts(i) > interp1([b a], ts([b a]), i))
            V(a,b) = 0;
        end
    end
end
end

% Bruno

Subject: Vectorizing nested for loops with loop dependencies

From: Michael

Date: 3 Jun, 2013 19:40:10

Message: 7 of 15

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <koiq9c$g5u$1@newscl01ah.mathworks.com>...
> What is typical value of len?
>
> Here is an equivalent function for weight > 0.
>
> function V = someFunction(ts)
> len = length(ts);
> i = 1:len;
> V = bsxfun(@minus,ts(:),ts(:)') ./ bsxfun(@minus,i(:),i(:)');
> V = tril(V,-1);
> for a = 1 : len,
> for b = 1 : a-1,
> i = b+1:a-1;
> if any(ts(i) > interp1([b a], ts([b a]), i))
> V(a,b) = 0;
> end
> end
> end
> end
>
> % Bruno

Hi Bruno,
thanks for your effort...The typical value for "len" is between 50 and several hundreds.

Subject: Vectorizing nested for loops with loop dependencies

From: dpb

Date: 3 Jun, 2013 20:25:02

Message: 8 of 15

On 6/3/2013 2:40 PM, Michael wrote:
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message
> <koiq9c$g5u$1@newscl01ah.mathworks.com>...
>> What is typical value of len?
>>
>> Here is an equivalent function for weight > 0.
>>
>> function V = someFunction(ts)
>> len = length(ts);
>> i = 1:len;
>> V = bsxfun(@minus,ts(:),ts(:)') ./ bsxfun(@minus,i(:),i(:)');
>> V = tril(V,-1);
>> for a = 1 : len,
>> for b = 1 : a-1,
>> i = b+1:a-1;
>> if any(ts(i) > interp1([b a], ts([b a]), i))
>> V(a,b) = 0;
>> end
>> end
>> end
>> end
>>
>> % Bruno
>
> Hi Bruno,
> thanks for your effort...The typical value for "len" is between 50 and
> several hundreds.

I was hoping maybe Bruno would be along before I got back...looks like
about what I had in mind as way to attack but didn't have the time to
actually try to write it then--does this suffice/work???

--

Subject: Vectorizing nested for loops with loop dependencies

From: Michael

Date: 3 Jun, 2013 22:28:10

Message: 9 of 15

dpb <none@non.net> wrote in message <koiu2h$r0d$1@speranza.aioe.org>...
> On 6/3/2013 2:40 PM, Michael wrote:
> > "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message
> > <koiq9c$g5u$1@newscl01ah.mathworks.com>...
> >> What is typical value of len?
> >>
> >> Here is an equivalent function for weight > 0.
> >>
> >> function V = someFunction(ts)
> >> len = length(ts);
> >> i = 1:len;
> >> V = bsxfun(@minus,ts(:),ts(:)') ./ bsxfun(@minus,i(:),i(:)');
> >> V = tril(V,-1);
> >> for a = 1 : len,
> >> for b = 1 : a-1,
> >> i = b+1:a-1;
> >> if any(ts(i) > interp1([b a], ts([b a]), i))
> >> V(a,b) = 0;
> >> end
> >> end
> >> end
> >> end
> >>
> >> % Bruno
> >
> > Hi Bruno,
> > thanks for your effort...The typical value for "len" is between 50 and
> > several hundreds.
>
> I was hoping maybe Bruno would be along before I got back...looks like
> about what I had in mind as way to attack but didn't have the time to
> actually try to write it then--does this suffice/work???
>
> --

Thanks Bruno and dpb for your great support.
I just tried out Bruno's code snippet - unfortunately it returned an error for line 9:
Error using >
Matrix dimensions must agree.
Error in someFunction (line 9)
        if any(ts(i) > interp1([b a], ts([b a]), i))

Subject: Vectorizing nested for loops with loop dependencies

From: dpb

Date: 3 Jun, 2013 22:58:35

Message: 10 of 15

On 6/3/2013 5:28 PM, Michael wrote:
...

> Thanks Bruno and dpb for your great support.
> I just tried out Bruno's code snippet - unfortunately it returned an
> error for line 9:
> Error using > Matrix dimensions must agree.
> Error in someFunction (line 9)
> if any(ts(i) > interp1([b a], ts([b a]), i))

Hmmmm....again I'm rushed so not sure the exact fixup since haven't read
Bruno's code carefully, but--

i=b+1:a-1;

so interp1() will return a vector of the length of that number of
elements whereas ts(i) is a single value. Using

if any(ts > interp1([b a], ts([b a]), i))

will fix that problem; problem is I don't know otomh that that is what
you really need at that point.

any() will turn it to a single T/F value, of course.

See if that gets you anywhere--supper's calling; I gotta' run. :)

--

Subject: Vectorizing nested for loops with loop dependencies

From: Michael

Date: 3 Jun, 2013 23:40:09

Message: 11 of 15

dpb <none@non.net> wrote in message <koj72j$jle$1@speranza.aioe.org>...
> On 6/3/2013 5:28 PM, Michael wrote:
> ...
>
> > Thanks Bruno and dpb for your great support.
> > I just tried out Bruno's code snippet - unfortunately it returned an
> > error for line 9:
> > Error using > Matrix dimensions must agree.
> > Error in someFunction (line 9)
> > if any(ts(i) > interp1([b a], ts([b a]), i))
>
> Hmmmm....again I'm rushed so not sure the exact fixup since haven't read
> Bruno's code carefully, but--
>
> i=b+1:a-1;
>
> so interp1() will return a vector of the length of that number of
> elements whereas ts(i) is a single value. Using
>
> if any(ts > interp1([b a], ts([b a]), i))
>
> will fix that problem; problem is I don't know otomh that that is what
> you really need at that point.
>
> any() will turn it to a single T/F value, of course.
>
> See if that gets you anywhere--supper's calling; I gotta' run. :)
>
> --

I'm afraid that didn't fix it - there's still the same error.

Subject: Vectorizing nested for loops with loop dependencies

From: dpb

Date: 4 Jun, 2013 00:31:22

Message: 12 of 15

On 6/3/2013 6:40 PM, Michael wrote:
...


> I'm afraid that didn't fix it - there's still the same error.

Well, need more than that for diagnostics...set a

debug stop on error

and then look at the various terms' sizes and see what's going on.

If that doesn't lead to nirvana, post the code and a sample (short) data
set that can just cut 'n paste and see.

_ALWAYS_ post full error in context instead of trying to describe a symptom.

--

Subject: Vectorizing nested for loops with loop dependencies

From: Bruno Luong

Date: 4 Jun, 2013 04:42:09

Message: 13 of 15

dpb <none@non.net> wrote in message <koj72j$jle$1@speranza.aioe.org>...

> so interp1() will return a vector of the length of that number of
> elements whereas ts(i) is a single value.

ts(i) is a vector, row or column depending on the orientation of ts.

Bruno

Subject: Vectorizing nested for loops with loop dependencies

From: Bruno Luong

Date: 4 Jun, 2013 04:48:09

Message: 14 of 15

Try:

function V = someFunction(ts)

ts = ts(:);
len = length(ts);
i = 1:len;
V = bsxfun(@minus,ts(:),ts(:)') ./ bsxfun(@minus,i(:),i(:)');
V = tril(V,-1);
for a = 1 : len,
    for b = 1 : a-1,
        if any(ts(b+1:a-1) > ts(b) + V(a,b)*(1:a-1-b)')
            V(a,b) = 0;
        end
    end
end

end

Subject: Vectorizing nested for loops with loop dependencies

From: Michael

Date: 4 Jun, 2013 10:54:09

Message: 15 of 15

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <kojri9$g08$1@newscl01ah.mathworks.com>...
> Try:
>
> function V = someFunction(ts)
>
> ts = ts(:);
> len = length(ts);
> i = 1:len;
> V = bsxfun(@minus,ts(:),ts(:)') ./ bsxfun(@minus,i(:),i(:)');
> V = tril(V,-1);
> for a = 1 : len,
> for b = 1 : a-1,
> if any(ts(b+1:a-1) > ts(b) + V(a,b)*(1:a-1-b)')
> V(a,b) = 0;
> end
> end
> end
>
> end

Sorry guys, I admit that I should have done "dbstop if error" the first time I encountered the error.
Anyway: Bruno, your last code snippet works fine. I also did some visual inspection of the matrix that is returned by the function - everything works as expected.
Comparing the old code to the new one, there's about a 30% reduction in calculation time (with n=50).

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