CODEWITHSUNDEEP
matlabvectorization
nashynash
nashynash

Reputation: 375

How to code a vectorized expression that depends on other vectorized expression?

If for example I have three expressions: A, B and C as follows:

A(i+1) = A(i) + C(i).k
B(i+1) = B(i) + A(i).h
C(i+1) = A(i) + B(i)

where k and h are some constants and m and n is the desired size of C. i is the previous obtained value, i+1 is the next value. Now, if I use for loop, then I can code it as:

A(1)= 2;
B(1)= 5;
C(1)= 3;
for i=1:10
    A(i+1) = A(i) + C(i)*2;
    B(i+1) = B(i) + A(i)*3;
    C(i+1) = A(i) + B(i);
end

And it works just fine. But I want to code it in a vector form, as in without having to use a loop. But the problem is I do not know how to get around the dependency of:

Upvotes: 5

Views: 206

Answers (2)

flawr
flawr

Reputation: 11628

First, forgive me for abusing Matlab syntax for expressing mathematical stuff.

Consider following code, where we do exactly the same as in your example. Note that A,B,C are the rows of X. 

X = zeros(3,N+1);
X(:,1) = [2,5,3];
M= [1,0,2;3,1,0;1,1,0];
for i=1:N
X(:,i+1) = M*X(:,i);
end

This is just a matrix vector notation of the above code. I think it is even slower. Note that we could also compute: X(:,i+1) = M^i * X(:,1) which is even slower.

Notice that we can use the eigenvalue decomposition:

[V,D] = eigs(M);
X(:,i+1) = [V*D*inv(V)]^i * X;

Therefore 

X(:,i+1) = V*D^i*inv(V) * X;

So V*D^i*inv(V) is an explicit formula for the i+1th term of X. I suggest computing those analytically, and plug the formula you get into your code again.

EDIT: I wrote some code that should be close to analyitcally solving the system, you can compare the runtimes. It seems in the end preallocation with your first method is still the fastest IF you need ALL the terms. If you only need one of them, my suggested method is certainly quicker.

clear;clc
N = 10000000;
tic
    A(1)= 2;
    B(1)= 5;
    C(1)= 3;
    A = zeros(1,N+1);
    B=A;C=A;
    for i=1:N
    A(i+1) = A(i) + C(i)*2;
    B(i+1) = B(i) + A(i)*3;
    C(i+1) = A(i) + B(i);
    end
toc

tic
    X = zeros(3,N+1);
    X(:,1) = [2,5,3];
    M= [1,0,2;3,1,0;1,1,0];
    for i=1:N
    X(:,i+1) = M*X(:,i);
    end
toc



tic
    M= [1,0,2;3,1,0;1,1,0];
    [V,D]=eig(M); 
    v=0:N;
    d=diag(D);
    B=bsxfun(@power,repmat(d,1,N+1),v);
    Y=bsxfun(@times,V * B, V \[2;5;3]);
toc


tic
    M= [1,0,2;3,1,0;1,1,0];
    [V,D]=eig(M); 
    v=0:N;
    d=diag(D);
    Y = ones(3,N+1);
    for i=1:N
    Y(:,i+1) = d.*Y(:,i);
    end
    Y=bsxfun(@times,V * B, V \[2;5;3]);
toc

Upvotes: 5

Dev-iL
Dev-iL

Reputation: 24169

Here's a matrix-based way to obtain the n-th value of the [A;B;C] vector. I wouldn't exactly call it vectorization, but this could speed things up considerably for you:

[A,B,C] = deal(zeros(11,1));
A(1)= 2;
B(1)= 5;
C(1)= 3;

%% // Original method
for k=1:10
  A(k+1) = A(k) + C(k)*2;
  B(k+1) = B(k) + A(k)*3;
  C(k+1) = A(k) + B(k);
end

%% // Matrix method:
%// [ A ]     [1  0  2][ A ]
%// | B |  =  |3  1  0|| B |
%// [ C ]     [1  1  0][ C ]
%//      i+1                i
%// 
%// [ A ]     [1  0  2][ A ]        [1  0  2]   ( [1  0  2][ A ] )
%// | B |  =  |3  1  0|| B |    =   |3  1  0| * ( |3  1  0|| B | )
%// [ C ]     [1  1  0][ C ]        [1  1  0]   ( [1  1  0][ C ] )
%//      i+2                i+1                                 i

%// Thus, this coefficient matrix taken to the n-th power, multiplied by the input 
%// vector will yield the values of A(n+1), B(n+1), and C(n+1):
M = [1 0 2
     3 1 0
     1 1 0];
isequal(M^10*[A(1);B(1);C(1)],[A(11);B(11);C(11)])

In reality you can use M to the appropriate power (positive or negative) to obtain any [A,B,C]n from any [A,B,C]k ...

Upvotes: 7

Related Questions