求助 程序的简化

zhailiangjun

各位大侠，这是我的写的一个程序，主要的目的是把一个(N+1)^2*(N+1)^2方阵约化至(N+1)*(N+1)的方阵。但是里面有几层的循环，做起来很浪费时间，请大家帮忙看看有没有办法调整一下。

这个程序的思路其实是很简单的，rho是一个(N+1)^2*(N+1)^2的方阵，是由(N+1)^2的列矩阵和(N+1)^2行矩阵相乘得到的。rho2是一个(N+1)*(N+1)的方阵，是由rho约化得来的。a（程序中用a1和a2表示）和b分别是两组长度为(N+1)的基矢，具有相同的形式，每一组中的第n个矢量可以写成一个单位方矩阵的第n列元素。rho2中的元素，
             rho2(k1,k2)=\sum [ kron(a(k1),b(i)) ]' * rho *kron(a(k2),b(i))。
kron(a(k1),b(i)) 表示 a 这组基矢中的第 k1 个矢量和 b 中第 i 个矢量的直积，\sum是对b的所有基矢的求和。
就是这样了，请大家帮忙看看吧。多谢了。
function rho2=reduce(psi1,psi11,N)
format long
% global N
rho=psi1*psi11;
% psi1是一(N+1)^2维列矩阵，psi11为一(N+1)^2维行矩阵
rho2=zeros(N+1,N+1);
I=eye(N+1,N+1);
for k1=1:N+1
    a1=I(:,k1);
    for k2=1:N+1
        a2=I(:,k2);
        for k3=1:N+1
            b=I(:,k3);
            phi1=kron(a1,b);
            phi2=kron(a2,b);
            rho2(k1,k2)=phi1'*rho*phi2+rho2(k1,k2);
        end
    end
end