如何解决使用矢量化求解多个线性系统
| 抱歉,如果这很明显,但我搜索了一会儿却没有找到任何东西(或错过了它)。 我正在尝试用A 4x4矩阵和B 4x1向量解形式为Ax = B的线性系统。 我知道对于单个系统,我可以使用mldivide
来获得x:x=A\\B
。
但是,我正在尝试解决大量系统(可能> 10000),并且我不愿意使用for循环,因为在许多MATLAB问题中,我被告知它明显比矩阵公式慢。
然后我的问题是:是否有一种方法可以使用A 4x4x N和B矩阵4x N的矢量化来解决Ax = B?
PS:我不知道这是否重要,但是B向量对于所有系统都是相同的。
解决方法
您应该使用for循环。如果
A
保持不变而B
发生变化,则预先计算因式分解并重新使用它可能会有好处。但是对于您的问题,其中A
改变而B
保持不变,则除了解决N个线性系统外别无选择。
您也不必太担心循环的性能成本:MATLAB JIT编译器意味着循环在最新版本的MATLAB上通常可以同样快。
, 我认为您无法进一步优化。正如@Tom所解释的那样,由于A
是一个变化,因此事先将各种A
\分解为无益。
给定您提到的尺寸,除了循环解决方案之外,它还非常快:
A = rand(4,4,10000);
B = rand(4,1); %# same for all linear systems
tic
X = zeros(4,size(A,3));
for i=1:size(A,3)
X(:,i) = A(:,:,i)\\B;
end
toc
经过的时间是0.168101秒。
, 这是问题所在:
您正在尝试在3D矩阵上执行2D操作(mldivide)。无论您如何看待它,都需要按索引引用矩阵,这是时间惩罚的开始……这不是问题所在的for循环,而是人们如何使用它们。
如果您可以以不同的方式构造问题,那么也许可以找到一个更好的选择,但是现在您有几个选择:
1-混合
2-并行处理(编写parfor循环)
3-CUDA
, 这是一个相当神秘的解决方案,它利用了MATLAB独特的优化功能。使用对角线的4x4块构造一个巨大的4k x 4k稀疏矩阵。然后同时解决所有问题。
在我的机器上,它获得与@ Amro / Tom的for循环解决方案相同的解决方案,单精度精度最高,但是速度更快。
n = size(A,1);
k = size(A,3);
AS = reshape(permute(A,[1 3 2]),n*k,n);
S = sparse( ...
repmat(1:n*k,n,1)\',...
bsxfun(@plus,reshape(repmat(1:n:n*k,1),[],0:n-1),...
AS,...
n*k,n*k);
X = reshape(S\\repmat(B,k,k);
对于一个随机的例子:
For k = 10000
For loop: 0.122570 seconds.
Giant sparse system: 0.032287 seconds.
如果您知道4x4矩阵是正定的,则可以在S上使用use11ѭ来提高精度。
真傻但是即使在使用JIT的情况下,matlab的for循环在2015年仍然有多慢。当k不太大时,此解决方案似乎找到了一个最佳位置,因此所有内容仍适合内存。
, 我知道这个职位已有几岁了,但是我还是会捐出两分钱。您可以将所有A矩阵放入一个较大的块对角矩阵中,其中大矩阵的对角线上将有4x4块。右边将是您的所有b向量一遍又一遍地堆叠在一起。设置好之后,它将被表示为稀疏系统,并且可以使用mldivide选择的算法来有效地解决。这些块在数值上是解耦的,因此,即使其中存在单个块,使用mldivide时,非单个块的答案也应该是正确的。在MATLAB Central上有采用这种方法的代码:
http://www.mathworks.com/matlabcentral/fileexchange/24260-multiple-same-size-linear-solver
我建议尝试看看这种方法是否比循环更快。我怀疑它会更有效,尤其是对于大量的小型系统。特别是,如果在N个矩阵上都有A系数的公式,则可以使用MATLAB向量运算(无循环)构建左手边,这可以进一步节省成本。正如其他人指出的那样,矢量化操作并不总是更快,但根据我的经验,它们通常是。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。