matlab 重复循环,在MATLAB中优化重复估计(当前是一个循环)

结果

我加快了您的原始版本,因为您的编辑3实际上不起作用(并且还做了不同的事情).

所以,在我的电脑上：

你的(原版)版本：8.428473秒.

我在下面给出了混淆的单行：0.964589秒.

首先,除了给人留下深刻印象之外,我会按照我的写作给出：

%%// Some example data

xdim = 500;

ydim = 500;

n_timepoints = 10; % for example

estimate = zeros(xdim,ydim); %// initialization with explicit size

picture = randn(xdim,ydim,n_timepoints);

%%// Your original solution

%// (slightly altered to make my version's results agree with yours)

tic

Y = randn(n_timepoints,xdim*ydim);

ii = 1;

for x = 1:xdim

for y = 1:ydim

X = squeeze(picture(x,y,:)); %// or similar creation of X matrix

B = (X'*X)^(-1)*X' * Y(:,ii);

ii = ii+1;

%// sometimes you keep everything and do

%// estimate(x,y,:) = B(:);

%// sometimes just the first element is important and you do

estimate(x,y) = B(1);

end

end

toc

%%// My version

tic

%// UNLEASH THE FURY!!

estimate2 = reshape(sparse(1:xdim*ydim*n_timepoints, ...

builtin('_paren', ones(n_timepoints,1)*(1:xdim*ydim),:), ...

builtin('_paren', permute(picture, [3 2 1]),:))\Y(:), ydim,xdim).'; %'

toc

%%// Check for equality

max(abs(estimate(:)-estimate2(:))) % (always less than ~1e-14)

分解

首先,这是您应该实际使用的版本：

%// Construct sparse block-diagonal matrix

%// (Type "help sparse" for more information)

N = xdim*ydim; %// number of columns

M = N*n_timepoints; %// number of rows

ii = 1:N;

jj = ones(n_timepoints,1)*(1:N);

s = permute(picture, [3 2 1]);

X = sparse(ii,jj(:), s(:));

%// Compute ALL the estimates at once

estimates = X\Y(:);

%// You loop through the *second* dimension first, so to make everything

%// agree, we have to extract elements in the "wrong" order, and transpose:

estimate2 = reshape(estimates, ydim,xdim).'; %'

下面是xdim = ydim = n_timepoints = 2的图片和相应矩阵X的示例：

>> clc, picture, full(X)

picture(:,:,1) =

-0.5643 -2.0504

-0.1656 0.4497

picture(:,:,2) =

0.6397 0.7782

0.5830 -0.3138

ans =

-0.5643 0 0 0

0.6397 0 0 0

0 -2.0504 0 0

0 0.7782 0 0

0 0 -0.1656 0

0 0 0.5830 0

0 0 0 0.4497

0 0 0 -0.3138

你可以看到为什么稀疏是必要的 – 它主要是零,但会很快变大.完整的矩阵将快速消耗所有RAM,而稀疏的矩阵不会消耗比原始图像矩阵更多的RAM.

有了这个矩阵X,就出现了新问题

X·b = Y

现在包含所有问题

X1 · b1 = Y1

X2 · b2 = Y2

...

哪里

b = [b1; b2; b3; ...]

Y = [Y1; Y2; Y3; ...]

所以,单一命令

X\Y

将立即解决您的所有系统.

这将所有艰苦的工作卸载到一组高度专业化,编译为机器特定的代码,优化的每路算法,而不是解释的,通用的,总是两步 – 远离MATLAB中的硬件循环.

将它转换为X是矩阵的版本应该是直截了当的;你最终会得到像blkdiag那样的东西,mldivide也可以使用与上面完全相同的方式.