加速矩阵计算

Question

我正在研究线性模型预测控制，我只需要计算控制器的一些矩阵。计算其中一个需要花费大量的时间，我想问一下，是否有更好的方法来编码这个计算。我正在使用MATLAB，但我也理解FORTRAN。

我想要计算一个矩阵(，Φ，)，但是计算它需要花费很长的时间。矩阵。矩阵的形式为(右)：Φ矩阵

这里是一本书，我在书中找到了这张图片，以防你需要参考(尤其是第8页)。

现在，我用MATLAB编写的代码如下所示：(在编辑之后移动它)

考虑到我将有相当大的NS，Np和Nc变量，它将花费大量的时间来做这个计算。有没有一种最佳的方法(或者至少比我的更好)来加速这个计算？

编辑

在考虑了@Daniel's & After 2682877的评论之后，我测试了这个

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)将越来越大。

丹尼尔的第二条评论，我知道我做的计算比我应该做的更多，但是我缺乏经验限制了我提升这个想法的能力。

@杜拉斯的评论，我不太熟悉parfor，但我会测试它。

参考答案:一旦我理解了你的建议，我就会对它们进行检验(.)然后回到你身边。

结果很明显，我最初的想法只是比这里建议的要差一点。谢谢你们！你真是帮了大忙！

Answer 1

计算C * A^(i-1-(j-1)) * B的结果是有限的，它只依赖于i和j之间的差异。为了不重复计算它，我的解决方案迭代这个差值，然后根据这两个变量计算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

业绩比较：

​

​

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