加速矩阵计算

Question

我正在从事线性模型预测控制，需要为控制器计算一些矩阵，但其中一个矩阵的计算时间很长，我想问一下是否有更好的编码方式来计算它。我正在使用MATLAB，但我也了解FORTRAN。
我想计算一个矩阵（Φ），但我现在的方法计算它需要太长时间。Φ矩阵的形式如下（右边那个）：MPC_matrices.
这里是我找到这张图片的书籍，如果需要参考可以看看第8页。
现在我在MATLAB中编写的代码如下：（移动到EDIT之后）
考虑到我的NS、Np和Nc变量会非常大，这个计算过程会花费很长时间。有没有更优化的方法（或者至少比我的方法更好）来加速这个计算过程呢？

EDIT
在考虑了@Daniel和user2682877的评论后，我测试了以下内容：

clear;clc
Np = 80;
Nc = Np / 2;
m  = 3;
q  = 1;
Niter = 30;
MAT = zeros(Niter,5);
for I=1:Niter
    NS = 10 * I;
    A = rand(NS,NS);
    B = rand(NS,m);
    C = rand(1,NS);
    tic
    Phi1 = zeros(Np*q,Nc*m);
    CB = C * B;
    for i=1:Np
        for j=1:Nc
            if j<i
                Phi1( (q*i-(q-1)):(q*i) , (m*j-(m-1)):(m*j) ) = C * A^(i-1-(j-1)) * B;
            elseif j==i
                Phi1( (q*i-(q-1)):(q*i) , (m*j-(m-1)):(m*j) ) = CB;
            end
        end
    end
    t1 = toc;

% 丹尼尔的建议

    tic
    Phi2=zeros(Np*q,Nc*m);
    CB = C * B;
    for diffij=0:Np-1
        if diffij>0
            F=C * A^diffij * B;
        else
            F=CB;
        end
        for i=max(1,diffij+1):min(Np,Nc+diffij)
            j=i-diffij;
            Phi2( (q*i-(q-1)):(q*i) , (m*j-(m-1)):(m*j) ) = F;
        end
    end
    t2 = toc;

% user2682877 的建议

    tic
    Phi3=zeros(Np*q,Nc*m);
    temp = B;
    % 1st column
    Phi3( (q*1-(q-1)):(q*1) , (m*1-(m-1)):(m*1) ) = C * B;
    for i=2:Np
        % reusing temp
        temp = A * temp;
        Phi3( (q*i-(q-1)):(q*i) , (m*1-(m-1)):(m*1) ) = C * temp;
    end
    % remaining columns
    for j=2:Nc
        for i=j:Nc
            Phi3( (q*i-(q-1)):(q*i) , (m*j-(m-1)):(m*j) ) =...
                Phi3( (q*(i-j+1)-(q-1)):(q*(i-j+1)) , (m*1-(m-1)):(m*1) );
        end
    end
    t3 = toc;

    MAT(I,:) = [I, NS, t1, t2 ,t3];
    fprintf('I iteration = %g \n', I);
end
figure(1)
clf(1)
hold on
plot(MAT(:,2),MAT(:,3),'b')
plot(MAT(:,2),MAT(:,4),'r')
plot(MAT(:,2),MAT(:,5),'g')
hold off
legend('My <Unfortunate> Idea','Daniel`s suggestion','user2682877 suggestion')
xlabel('NS variable')
ylabel('Time, s')

这是生成的图形： 请记住，现在NS = 300，但当我升级我的模型（我打算在状态空间模型中包含越来越多的方程和变量）时，这些变量（主要是NS和Np）将越来越大。 关于@Daniel的第二个评论，我知道我执行的计算比我应该执行的计算多，但我的经验不足限制了我对升级的想法。 关于@durasm的评论，我不太熟悉parfor，但我会测试它。 关于答案：只要我理解了你们的建议，我会尽快测试并回复你们。

结果 很明显，我的最初想法只比这里建议的稍差一点..谢谢大家！你们非常有帮助！

- ASK22

---

你能否用一些示例数据来完善你的代码呢？可能随机数据就可以了。根据你的描述，我尝试了 NS=100; m=99; A = rand(NS,NS); B = rand(NS,m); C = rand(1,NS); 但这并没有初始化所有内容。 - Daniel

你肯定可以进行优化。对于 j<i，C * A^(i-1-(j-1)) * B 只有 Np 种不同的结果，但你却大约计算了 Np*Nc/2 次。 - Daniel

我认为parfor不是正确的选择。这段代码涉及大量矩阵乘法，需要使用多核代码。 - Daniel

@ASK22: 比较返回的Phi，还有一点问题是用户2682877的解决方案没有填满完整的矩阵。 - Daniel

@Daniel，某些元素的结果似乎有误，但我会在我（或user2682877）修复后尽快更新结果。谢谢。 - ASK22

Answer 1

你的计算结果只能在有限的一组 C * A^(i-1-(j-1)) * B 中获得，这些结果仅取决于 i 和 j 之间的差异。为了避免重复计算，我的解决方案迭代这个差异和 i，然后根据这两个变量计算 j。

Phi=zeros(Np*q,Nc*m);
CB = C * B;
for diffij=0:Np-1
    if diffij>0
        F=C * A^diffij * B;
    else
        F=CB;
    end
    for i=max(1,diffij+1):min(Np,Nc+diffij)
        j=i-diffij;

        Phi( (q*i-(q-1)):(q*i) , (m*j-(m-1)):(m*j) ) = F;

    end
end

性能比较：

（注意：此处为原文，无需翻译）

- Daniel

---

评论去哪了？这是一个好建议，可以尝试在循环中使用X=B和X=A*X，而不是重复使用A^diffij*B进行乘法。不确定是否更快，但至少值得一试。 - Daniel

Answer 2

也许你可以尝试这个：

Phi=zeros(Np*q,Nc*m);
temp = B;
% 1st column
Phi( (q*1-(q-1)):(q*1) , (m*1-(m-1)):(m*1) ) = C * B;
for i=2:Np
    % reusing temp
    temp = A * temp;
    Phi( (q*i-(q-1)):(q*i) , (m*1-(m-1)):(m*1) ) = C * temp;
end
% remaining columns
for j=2:Nc
    for i=j:Np
        Phi( (q*i-(q-1)):(q*i) , (m*j-(m-1)):(m*j) ) = Phi( (q*(i-j+1)-(q-1)):(q*(i-j+1)) , (m*1-(m-1)):(m*1) );
    end
end

- user2682877

---

我认为内部循环存在问题。在原始代码中，它是 i=1:Np。 - Daniel

@daniel 迭代 i=1 是在上一行完成的。 - user2682877

我不确定这个答案有什么问题，但运行问题中的基准代码会产生不同的结果。 - Daniel

确实你是对的！我在内部循环中写成了Nc而不是Np。我已经编辑了我的帖子。 - user2682877

我看到你更新了代码，但它仍然生成错误的值。 - Daniel