How to code a vectorized expression that depends o

问题:

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:

• A on its previous value and previous C value
• B on it previous values and previous C value of A
• C on the previous values of A and B

回答1:

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 ...

回答2:

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