Thread Subject: create matrix as a function of previous columns

I would like to create a matrix in which each column is a function of the
previous column. It is probably easiest to describe what I would like to
do with the following code :

M is 20 x 100 matrix
k1 is 20 x 1 vector of constants
k2 is 20 x 1 vector of constants
k3 is 20 x 1 vector of constants

The first column of M is 'primed' with values.

for c = 2:100; % step through each column of M starting at column 2
M(:, c) = (k1 .* M(:, c-1)) + (k2 .* M(:, c-1)) + k3;
end;

This portion of code is called several times (with larger dimensions) and
is causing quite a slowdown and I would like to vectorize it to improve
performance.

I've been trying to do it with combinations of cumsum() and cumprod() and
I've gotten close but can't get it quite right.

Any suggestions would be greatly appreciated.

Thank you,
Curt

Curt <97wideglide@gmail.com> writes:

> I would like to create a matrix in which each column is a function of the
> previous column. It is probably easiest to describe what I would like to
> do with the following code :
>
> M is 20 x 100 matrix
> k1 is 20 x 1 vector of constants
> k2 is 20 x 1 vector of constants
> k3 is 20 x 1 vector of constants
>
> The first column of M is 'primed' with values.
>
> for c = 2:100; % step through each column of M starting at column 2
> M(:, c) = (k1 .* M(:, c-1)) + (k2 .* M(:, c-1)) + k3;
> end;
>
> This portion of code is called several times (with larger dimensions) and
> is causing quite a slowdown and I would like to vectorize it to improve
> performance.
>
> I've been trying to do it with combinations of cumsum() and cumprod() and
> I've gotten close but can't get it quite right.

You could start with eliminating the excess operations. k1 and k2
associate, so:

M(:,c) = (k1+k2) .* M(:,c-1) + k3;

Maybe you could write down the closed-form expansion of this expression
and compute it directly?

kM + c
(kkM + kc) + c
(kkkM + kkc + kc) + c

I'm not doing this carefully, but I think it comes out to (for the jth column)

Mk.^j + sum_i=0^{j-1} c k^i

Mk.^j + c (k^j - 1) / (k^(j-1) - 1)

or, back to your problem

M(:,j) = M(:,1).*(k1+k2).^j + (k3.^j - 1) ./ (k3.^(j-1) - 1)


So now you can repmat your k1,k2,k3 columns, generate a j matrix for the
whole thing, and do this all in one expression.

Looking at this again, I think I'm probably off by one on the j parameter.


-Peter

Peter Boettcher <boettcher@ll.mit.edu> wrote in message
<muyod4xbuk2.fsf@G99-Boettcher.llan.ll.mit.edu>...
> Curt <97wideglide@gmail.com> writes:
>
> > I would like to create a matrix in which each column is a function of the
> > previous column. It is probably easiest to describe what I would like to
> > do with the following code :
> >
> > M is 20 x 100 matrix
> > k1 is 20 x 1 vector of constants
> > k2 is 20 x 1 vector of constants
> > k3 is 20 x 1 vector of constants
> >
> > The first column of M is 'primed' with values.
> >
> > for c = 2:100; % step through each column of M starting at column 2
> > M(:, c) = (k1 .* M(:, c-1)) + (k2 .* M(:, c-1)) + k3;
> > end;
> >
> > This portion of code is called several times (with larger dimensions) and
> > is causing quite a slowdown and I would like to vectorize it to improve
> > performance.
> >
> > I've been trying to do it with combinations of cumsum() and cumprod()
and
> > I've gotten close but can't get it quite right.
>
> You could start with eliminating the excess operations. k1 and k2
> associate, so:
>
> M(:,c) = (k1+k2) .* M(:,c-1) + k3;
>
> Maybe you could write down the closed-form expansion of this expression
> and compute it directly?
>
> kM + c
> (kkM + kc) + c
> (kkkM + kkc + kc) + c
>
> I'm not doing this carefully, but I think it comes out to (for the jth column)
>
> Mk.^j + sum_i=0^{j-1} c k^i
>
> Mk.^j + c (k^j - 1) / (k^(j-1) - 1)
>
> or, back to your problem
>
> M(:,j) = M(:,1).*(k1+k2).^j + (k3.^j - 1) ./ (k3.^(j-1) - 1)
>
>
> So now you can repmat your k1,k2,k3 columns, generate a j matrix for the
> whole thing, and do this all in one expression.
>
> Looking at this again, I think I'm probably off by one on the j parameter.
>
>
> -Peter

It looks like a call to bsxfun to compute the powers,
then a cumsum.

John

On Wed, 16 Jul 2008 13:36:13 -0400, Peter Boettcher wrote:

On Wed, 16 Jul 2008 13:36:13 -0400, Peter Boettcher wrote:

> Curt <97wideglide@gmail.com> writes:
>
>> I would like to create a matrix in which each column is a function of
>> the previous column. It is probably easiest to describe what I would
>> like to do with the following code :
>>
>> M is 20 x 100 matrix
>> k1 is 20 x 1 vector of constants
>> k2 is 20 x 1 vector of constants
>> k3 is 20 x 1 vector of constants
>>
>> The first column of M is 'primed' with values.
>>
>> for c = 2:100; % step through each column of M starting at column 2
>> M(:, c) = (k1 .* M(:, c-1)) + (k2 .* M(:, c-1)) + k3;
>> end;
>>
>> This portion of code is called several times (with larger dimensions)
>> and is causing quite a slowdown and I would like to vectorize it to
>> improve performance.
>>
>> I've been trying to do it with combinations of cumsum() and cumprod()
>> and I've gotten close but can't get it quite right.
>
> You could start with eliminating the excess operations. k1 and k2
> associate, so:
>
> M(:,c) = (k1+k2) .* M(:,c-1) + k3;
>
> Maybe you could write down the closed-form expansion of this expression
> and compute it directly?
>
> kM + c
> (kkM + kc) + c
> (kkkM + kkc + kc) + c
>
> I'm not doing this carefully, but I think it comes out to (for the jth
> column)
>
> Mk.^j + sum_i=0^{j-1} c k^i
>
> Mk.^j + c (k^j - 1) / (k^(j-1) - 1)
>
> or, back to your problem
>
> M(:,j) = M(:,1).*(k1+k2).^j + (k3.^j - 1) ./ (k3.^(j-1) - 1)
>
>
> So now you can repmat your k1,k2,k3 columns, generate a j matrix for the
> whole thing, and do this all in one expression.
>
> Looking at this again, I think I'm probably off by one on the j
> parameter.
>
>
> -Peter

Thank you very much for your response. I made a small mistake in the
description of my problem which changed it slightly. Sadly, I've been
working on this all day (and night) and now to solve it I think I need a
solution to another related problem.

The mistake I made is that K1 is actually a 20 x 100 matrix of constants
instead of a vector so I have to consider each individual Kn instead of
just repmating K1 and raising to the jth power.

The corrected code is :

for c = 2:100; % step through each column of M starting at column 2
  M(:, c) = (k1(:, c-1) .* M(:, c-1)) + (k2 .* M(:, c-1)) + k3;
end;

and the closed form expansion then looks like :

M
K1M + c
K2K1M + K2c + c
K3K2K1M + K3K2c + K3c + c
K4K3K2K1M + K4K3K2c + K4K3c + K4c + c

and for the jth column

(Kn...K1M) + (Kn...K2c) + (Kn...K3c) + ... + Knc + c

The first operand (Kn...K1M) is straightforward but the others have me
puzzled.

So, the problem I think I need to solve is given a matrix of values how
can I produce another matrix who's elements are the product of the
current element and the element in the previous column plus 1 ?

the following matrix :

v1 v4 v7
v2 v5 v8
v3 v6 v9

would produce :

v1 v4(v1+1) v7(v4(v1+1) +1)
v2 v5(v2+1) v8(v5(v2+1) +1)
v3 v6(v3+1) v9(v6(v3+1) +1)

I would really appreciate any advice anyone has for me.

Thanks,
Curt