我发现自己需要对图像中的每个像素进行最小二乘（或类似的基于矩阵的运算）。每个像素都有一组与之关联的数字，因此可以将其排列为 3D 矩阵。
（下一位可以跳过）
快速解释我所说的最小二乘估计的含义：
假设我们有一些由 Y = Ax^2 + Bx + C 建模的二次系统，我们正在寻找那些 A、B、C 系数。对于 X 和相应 Y 的几个样本（至少 3 个），我们可以通过以下方式估计它们：
- 将（比如说 10 个）X 样本排列成一个矩阵，例如
X = [x(:).^2 x(:) ones(10,1)]; - 将 Y 个样本排列成一个相似的矩阵：
Y = y(:); - 通过求解来估计系数 A、B、C：
coeffs = (X'*X)^(-1)*X'*Y;
如果您愿意，请自行尝试：
A = 5; B = 2; C = 1;
x = 1:10;
y = A*x(:).^2 + B*x(:) + C + .25*randn(10,1); % added some noise here
X = [x(:).^2 x(:) ones(10,1)];
Y = y(:);
coeffs = (X'*X)^-1*X'*Y
coeffs =
5.0040
1.9818
0.9241
如果我在那里失去了你，请重新开始关注
*主要改写*我已对其进行了修改，使其尽可能接近我遇到的实际问题，并且仍然使其成为最小的工作示例。
问题设置
%// Setup
xdim = 500;
ydim = 500;
ncoils = 8;
nshots = 4;
%// matrix size for each pixel is ncoils x nshots (an overdetermined system)
%// each pixel has a matrix stored in the 3rd and 4rth dimensions
regressor = randn(xdim,ydim, ncoils,nshots);
regressand = randn(xdim, ydim,ncoils);
所以我的问题是我必须对图像中的每个像素进行 (X'*X)^-1*X'*Y （最小二乘或类似）操作。虽然它本身是矢量化/矩阵化的，但我必须为每个像素执行此操作的唯一方法是在 for 循环中，例如：
原始代码风格
%// Actual work
tic
estimate = zeros(xdim,ydim);
for col=1:size(regressor,2)
for row=1:size(regressor,1)
X = squeeze(regressor(row,col,:,:));
Y = squeeze(regressand(row,col,:));
B = X\Y;
% B = (X'*X)^(-1)*X'*Y; %// equivalently
estimate(row,col) = B(1);
end
end
toc
Elapsed time = 27.6 seconds
针对评论和其他想法进行编辑
我尝试了一些事情：
1. 重新塑造成一个长向量并删除了双for循环。这节省了一些时间。
2. 预先通过 -ing 图片删除squeeze（和在线转置）permute：这节省了更多时间。
当前示例：
%// Actual work
tic
estimate2 = zeros(xdim*ydim,1);
regressor_mod = permute(regressor,[3 4 1 2]);
regressor_mod = reshape(regressor_mod,[ncoils,nshots,xdim*ydim]);
regressand_mod = permute(regressand,[3 1 2]);
regressand_mod = reshape(regressand_mod,[ncoils,xdim*ydim]);
for ind=1:size(regressor_mod,3) % for every pixel
X = regressor_mod(:,:,ind);
Y = regressand_mod(:,ind);
B = X\Y;
estimate2(ind) = B(1);
end
estimate2 = reshape(estimate2,[xdim,ydim]);
toc
Elapsed time = 2.30 seconds (avg of 10)
isequal(estimate2,estimate) == 1;
罗迪奥尔登休斯的方式
N = xdim*ydim*ncoils; %// number of columns
M = xdim*ydim*nshots; %// number of rows
ii = repmat(reshape(1:N,[ncoils,xdim*ydim]),[nshots 1]); %//column indicies
jj = repmat(1:M,[ncoils 1]); %//row indicies
X = sparse(ii(:),jj(:),regressor_mod(:));
Y = regressand_mod(:);
B = X\Y;
B = reshape(B(1:nshots:end),[xdim ydim]);
Elapsed time = 2.26 seconds (avg of 10)
or 2.18 seconds (if you don't include the definition of N,M,ii,jj)
所以问题是：
有没有（甚至）更快的方法？
（我不这么认为。）