在 MATLAB 中向量化线性方程组的解

Question

总结：这个问题涉及线性回归计算算法的改进。

---

我有一个 3D (dlMAT) 数组，表示在不同曝光时间拍摄的同一场景的单色照片（矢量 IT）。在数学上，dlMAT 的第三维上的每个向量都代表一个需要解决的单独的线性回归问题。需要估计系数的方程为：

DL = R*IT^P，其中DL和IT是通过实验获得的，R和P必须是估计的。

上面的等式经过一个对数可以转化为一个简单的线性模型：

log(DL) = log(R) + P*log(IT)  =>  y = a + b*x

下面介绍的是求解这个方程组的最“天真”的方法，它本质上涉及迭代所有“第三维向量”并将1 到(IT,DL(ind1,ind2,:) 阶的多项式拟合：

%// Define some nominal values:
R = 0.3;
IT = 600:600:3000;
P = 0.97;
%// Impose some believable spatial variations:
pMAT = 0.01*randn(3)+P;
rMAT = 0.1*randn(3)+R;
%// Generate "fake" observation data:
dlMAT = bsxfun(@times,rMAT,bsxfun(@power,permute(IT,[3,1,2]),pMAT));
%// Regression:
sol = cell(size(rMAT)); %// preallocation
for ind1 = 1:size(dlMAT,1)
  for ind2 = 1:size(dlMAT,2)
    sol{ind1,ind2} = polyfit(log(IT(:)),log(squeeze(dlMAT(ind1,ind2,:))),1);
  end
end
fittedP = cellfun(@(x)x(1),sol);      %// Estimate of pMAT
fittedR = cellfun(@(x)exp(x(2)),sol); %// Estimate of rMAT

上述方法似乎是矢量化的一个很好的候选，因为它没有利用 MATLAB 的主要优势，即矩阵运算。出于这个原因，它不能很好地扩展，并且执行时间比我认为的要长。

---

存在基于矩阵除法执行此计算的替代方法，如 here 和 here 所示，它们涉及以下内容：

sol = [ones(size(x)),log(x)]\log(y);

也就是说，将1s 的向量附加到观测值后，再附加mldivide 来求解方程组。

我面临的主要挑战是如何使我的数据适应算法（反之亦然）。

问题 #1：如何扩展基于矩阵除法的解决方案以解决上述问题（并可能替换我正在使用的循环）？

问题 #2（奖励）：这种基于矩阵除法的解决方案背后的原理是什么？

Answer 1

包含矩阵除法的解决方案背后的秘密成分是Vandermonde matrix。该问题讨论了一个线性问题（线性回归），这些问题总是可以表述为矩阵问题，\ (mldivide) 可以在均方误差意义上解决^‡ .这种算法解决了类似的问题，在this answer 中进行了演示和解释。

以下是基准代码，将原始解决方案与聊天中建议的两个替代方案进行比较^1、2：

function regressionBenchmark(numEl)
clc
if nargin<1, numEl=10; end
%// Define some nominal values:
R = 5;
IT = 600:600:3000;
P = 0.97;
%// Impose some believable spatial variations:
pMAT = 0.01*randn(numEl)+P;
rMAT = 0.1*randn(numEl)+R;
%// Generate "fake" measurement data using the relation "DL = R*IT.^P"
dlMAT = bsxfun(@times,rMAT,bsxfun(@power,permute(IT,[3,1,2]),pMAT));
%% // Method1: loops + polyval
disp('-------------------------------Method 1: loops + polyval')
tic; [fR,fP] = method1(IT,dlMAT); toc;
fprintf(1,'Regression performance:\nR: %d\nP: %d\n',norm(fR-rMAT,1),norm(fP-pMAT,1));
%% // Method2: loops + Vandermonde
disp('-------------------------------Method 2: loops + Vandermonde')
tic; [fR,fP] = method2(IT,dlMAT); toc;
fprintf(1,'Regression performance:\nR: %d\nP: %d\n',norm(fR-rMAT,1),norm(fP-pMAT,1));
%% // Method3: vectorized Vandermonde
disp('-------------------------------Method 3: vectorized Vandermonde')
tic; [fR,fP] = method3(IT,dlMAT); toc;
fprintf(1,'Regression performance:\nR: %d\nP: %d\n',norm(fR-rMAT,1),norm(fP-pMAT,1));
function [fittedR,fittedP] = method1(IT,dlMAT)
sol = cell(size(dlMAT,1),size(dlMAT,2));
for ind1 = 1:size(dlMAT,1)
  for ind2 = 1:size(dlMAT,2)
    sol{ind1,ind2} = polyfit(log(IT(:)),log(squeeze(dlMAT(ind1,ind2,:))),1);
  end
end

fittedR = cellfun(@(x)exp(x(2)),sol);
fittedP = cellfun(@(x)x(1),sol);

function [fittedR,fittedP] = method2(IT,dlMAT)
sol = cell(size(dlMAT,1),size(dlMAT,2));
for ind1 = 1:size(dlMAT,1)
  for ind2 = 1:size(dlMAT,2)
    sol{ind1,ind2} = flipud([ones(numel(IT),1) log(IT(:))]\log(squeeze(dlMAT(ind1,ind2,:)))).'; %'
  end
end

fittedR = cellfun(@(x)exp(x(2)),sol);
fittedP = cellfun(@(x)x(1),sol);

function [fittedR,fittedP] = method3(IT,dlMAT)
N = 1; %// Degree of polynomial
VM = bsxfun(@power, log(IT(:)), 0:N); %// Vandermonde matrix
result = fliplr((VM\log(reshape(dlMAT,[],size(dlMAT,3)).')).');
%// Compressed version:
%// result = fliplr(([ones(numel(IT),1) log(IT(:))]\log(reshape(dlMAT,[],size(dlMAT,3)).')).');

fittedR = exp(real(reshape(result(:,2),size(dlMAT,1),size(dlMAT,2))));
fittedP = real(reshape(result(:,1),size(dlMAT,1),size(dlMAT,2)));

之所以可以将方法2向量化成方法3，本质上是矩阵乘法可以通过第二个矩阵的列进行分隔。如果A*B 产生矩阵X，则根据定义A*B(:,n) 为任何n 提供X(:,n)。将A 与mldivide 一起移动到右侧，这意味着A\X(:,n) 的划分可以一次性完成所有n 和A\X。 overdetermined system（线性回归问题）也是如此，其中一般没有精确解，mldivide 找到 minimizes the mean-square error 的矩阵。在这种情况下，A\X(:,n)（方法 2）的操作也可以一次性完成所有 n 和 A\X（方法 3）。

在增加dlMAT 的大小时改进算法的含义如下：

对于500*500（或2.5E5）元素，从Method 1到Method 3的加速大约是x3500！

观察profile的output也很有趣（这里是针对500*500的情况）：

方法一

方法二

方法3

从上面可以看出，通过squeeze 和flipud 重新排列元素占用了Method 2 运行时间的大约一半（！）。还可以看出，将解从单元格转换为矩阵会浪费一些时间。

由于第 3 种解决方案完全避免了所有这些陷阱以及循环（这主要意味着在每次迭代时重新评估脚本） - 不出所料，它会显着加快速度。

注意事项：

1. Method 3 的“压缩”版本和“显式”版本之间几乎没有区别，有利于“显式”版本。因此，它未包含在比较中。
2. 尝试了一个解决方案，其中Method 3 的输入为gpuArray-ed。这并没有提供改进的性能（甚至在某种程度上降低了性能），可能是由于错误的实现，或者与在 RAM 和 VRAM 之间来回复制矩阵相关的开销。