如何编码依赖于其他矢量化表达式的矢量化表达式？

Question

举例来说,如果我有三个表达式:A,B和C如下:

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

其中k和h是一些常数,m并且n是所需的大小C.i是先前获得的值,i+1是下一个值.现在,如果我使用for循环,那么我可以将其编码为:

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

它工作得很好.但我想用矢量形式编码,就像在不必使用循环一样.但问题是我不知道如何绕过依赖:

• A在其先前的值和以前的C值
• B它的前一个值和之前的C值A
• C在以前的值A和B

Answer 1

这是一种基于矩阵的方法来获得向量的n第th个值[A;B;C].我不会完全称它为矢量化,但这可以为你加速:

[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)])

实际上,您可以使用M适当的功率(正或负)从任何_k获得任何[A,B,C]_n ...[A,B,C]

Answer 2

首先,请原谅我滥用Matlab语法来表达数学内容.

请考虑以下代码,我们的示例与您的示例完全相同.注意A,B,C是行的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

这只是上述代码的矩阵向量符号.我认为它甚至更慢.请注意,我们还可以计算:X(:,i+1) = M^i * X(:,1)哪个更慢.

请注意,我们可以使用特征值分解:

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

因此

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

所以,V*D^i*inv(V)对于一个明确的公式i+1中日长期X.我建议用分析方法计算这些,并再次插入你的代码中的公式.

编辑:我写了一些应该接近分析解决系统的代码,你可以比较运行时.这似乎与你的第一个方法结束预分配仍然是最快如果你需要的所有条款.如果你只需要其中一个,我建议的方法肯定更快.

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