在MATLAB中向量化线性方程组的解 [英] Vectorizing the solution of a linear equation system in MATLAB

问题描述

总结:该问题涉及线性回归计算算法的改进.

Summary: This question deals with the improvement of an algorithm for the computation of linear regression.

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

I have a 3D (dlMAT) array representing monochrome photographs of the same scene taken at different exposure times (the vector IT) . Mathematically, every vector along the 3rd dimension of dlMAT represents a separate linear regression problem that needs to be solved. The equation whose coefficients need to be estimated is of the form:

DL = R*IT^P，其中DL和IT是通过实验获得的，并且必须估算R和P.

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

在应用对数后，上述方程式可以转换为简单的线性模型:

The above equation can be transformed into a simple linear model after applying a logarithm:

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

下面介绍的是解决此方程组的最幼稚"的方法，该方法本质上涉及对所有三维向量"进行迭代，并将阶次为1的多项式拟合为(IT,DL(ind1,ind2,:):

Presented below is the most "naive" way to solve this system of equations, which essentially involves iterating over all "3rd dimension vectors" and fitting a polynomial of order 1 to (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运算.因此，它的伸缩性不是很好，执行时间比我想象的要长得多.

The above approach seems like a good candidate for vectorization, since it does not utilize MATLAB's main strength that is MATrix operations. For this reason, it does not scale very well and takes much longer to execute than I think it should.

存在基于矩阵除法执行此计算的替代方法，如这里，其中包含以下内容:

There exist alternative ways to perform this computation based on matrix division, as demonstrated here and here, which involve something like this:

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

也就是说，将一个1的向量附加到观测值，然后跟随 mldivide 求解方程组.

That is, appending a vector of 1s to the observations, followed by mldivide to solve the equation system.

我面临的主要挑战是如何使我的数据适应算法(反之亦然).

The main challenge I'm facing is how to adapt my data to the algorithm (or vice versa).

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

Question #1: How can the matrix-division-based solution be extended to solve the problem presented above (and potentially replace the loops I am using)?

问题2(奖励):这种基于矩阵除法的解决方案的原理是什么?

Question #2 (bonus): What is the principle behind this matrix-division-based solution?

推荐答案

包括矩阵除法在内的解决方案的秘密要素是• . 此答案中演示并解释了解决此类问题的算法.

The secret ingredient behind the solution that includes matrix division is the 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} :

Below is benchmarking code that compares the original solution with two alternatives suggested in chat^{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).用mldivide将A移到右侧，这意味着使用A\X一次可以对所有n进行除法A\X(:,n). 超定系统(线性回归问题)也是如此，其中没有确切的解决方案通常，mldivide会找到最小化均方误差的矩阵..同样在这种情况下，对于n和A\X(方法3)，所有n都可以一次性完成A\X(:,n)(方法2)操作.

The reason why method 2 can be vectorized into method 3 is essentially that matrix multiplication can be separated by the columns of the second matrix. If A*B produces matrix X, then by definition A*B(:,n) gives X(:,n) for any n. Moving A to the right-hand side with mldivide, this means that the divisions A\X(:,n) can be done in one go for all n with A\X. The same holds for an 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. In this case too, the operations A\X(:,n) (method 2) can be done in one go for all n with A\X (method 3).

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

The implications of improving the algorithm when increasing the size of dlMAT can be seen below:

对于500*500(或2.5E5)元素，从Method 1到Method 3的加速约为 x3500 ！

For the case of 500*500 (or 2.5E5) elements, the speedup from Method 1 to Method 3 is about x3500!

观察 profile 的输出(在这里，对于500 * 500):

It is also interesting to observe the output of profile (here, for the case of 500*500):

从上面可以看到，通过 squeeze 和 flipud 约占一半(！) Method 2的运行时.还可以看出，溶液从细胞到基质的转化损失了一些时间.

From the above it is seen that rearranging the elements via squeeze and flipud takes up about half (!) of the runtime of Method 2. It is also seen that some time is lost on the conversion of the solution from cells to matrices.

由于第三个解决方案避免了所有这些陷阱，并且完全避免了循环(这主要意味着每次迭代都需要重新评估脚本)，因此毫无疑问会大大提高速度.

Since the 3rd solution avoids all of these pitfalls, as well as the loops altogether (which mostly means re-evaluation of the script on every iteration) - it unsurprisingly results in a considerable speedup.

1. Method 3的压缩"版本和显式"版本之间几乎没有什么区别，而支持显式"版本.因此，它没有包含在比较中.
2. 尝试了对Method 3的输入进行-c编辑的解决方案.这样做可能不会由于错误的实现或与在RAM和VRAM之间来回复制矩阵而产生的开销而提高了性能(甚至使其性能有所下降).

1. There was very little difference between the "compressed" and the "explicit" versions of Method 3 in favor of the "explicit" version. For this reason it was not included in the comparison.
2. A solution was attempted where the inputs to Method 3 were gpuArray-ed. This did not provide improved performance (and even somewhat degradaed them), possibly due to wrong implementation, or the overhead associated with copying matrices back and forth between RAM and VRAM.

这篇关于在MATLAB中向量化线性方程组的解的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！