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

Question

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

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

DL = R*IT^P, where DL and IT are obtained experimentally and R and P must be estimated.

上式可以在应用对数后转化为简单的线性模型：

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

下面介绍的是求解该方程组的最 "naive" 方法，它主要涉及迭代所有“第三维向量”并拟合阶数 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 的主要优势，即 MATrix 运算。出于这个原因，它不能很好地扩展并且执行时间比我认为的要长得多。

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

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

即在观测值上附加一个1s的向量，然后mldivide求解方程组。

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

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

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

Answer 1

包含矩阵除法的解决方案背后的秘密成分是 Vandermonde matrix. The question discusses a linear problem (linear regression), and those can always be formulated as a matrix problem, which \ (mldivide) can solve in a mean-square error sense^‡. Such an algorithm, solving a similar problem, is demonstrated and explained in 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) 可以一次性完成所有 A\X 的 n。 overdetermined system (linear regression problem), in which there is no exact solution in general, and mldivide finds the matrix that minimizes the mean-square error 也是如此。同样在这种情况下，操作 A\X(:,n)（方法 2）可以一次性完成所有 n 和 A\X（方法 3）。

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

对于500*500（或2.5E5）个元素的情况，从Method 1到Method 3的加速比约为x3500!

观察output of profile也很有意思（这里以500*500为例）：

方法一

方法二

方法三

从上面可以看出，通过 squeeze and flipud 重新排列元素大约占 Method 2 运行时间的一半 (!)。还可以看出，将解决方案从单元格转换为矩阵会损失一些时间。

由于第 3 种解决方案避免了所有这些陷阱以及整个循环（这主要意味着在每次迭代时重新评估脚本）- 它毫不奇怪地导致了相当大的加速。

备注：

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

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

Vectorizing the solution of a linear equation system in MATLAB

statistics

optimization

matlab

regression

vectorization

方法一

方法二

方法三

备注：