如何编写依赖于其他向量化表达式的向量化表达式？

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

这是一种基于矩阵的方法来获取[A;B;C] 向量的n-th 值。我不会完全称它为矢量化，但这可以大大加快你的速度：

[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 的适当幂（正或负）从任何[A,B,C]_k 获取任何[A,B,C]_n ...

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) 是X 的i+1th 项的显式公式。我建议对这些进行分析计算，然后将你得到的公式再次插入你的代码中。

编辑：我写了一些应该接近分析解决系统的代码，您可以比较运行时。看来最终用你的第一种方法预分配仍然是最快的 IF 你需要 ALL 的条件。如果你只需要其中一个，我建议的方法肯定更快。

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