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

简介:该问题涉及用于线性回归计算的算法的改进.

我有一个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的主要优势,即MATrix操作.出于这个原因,它不能很好地扩展,并且执行时间比我想象的要长得多.

存在基于矩阵划分执行此计算的替代方法,如此处和此处所示,其涉及如下内容:

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

也就是说,1在观测值上附加s 矢量,然后mldivide求解方程组.

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

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

问题#2(奖金):基于矩阵划分的解决方案背后的原理是什么？

包含矩阵划分的解决方案背后的秘密成分是Vandermonde矩阵.这个问题讨论了线性问题(线性回归),并且这些问题总是可以表示为矩阵问题,\(mldivide)可以在均方误差意义上解决^‡.解决类似问题的这种算法在这个答案中得到了证明和解释.

下面是比较两种选择原来的解决方案聊天建议基准测试代码^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)给出X(:,n)任何n.移动A到右手侧mldivide,这意味着该部门A\X(:,n)可以在一个随时随地都可以做n有A\X.对于超定系统(线性回归问题)也是如此,其中一般没有精确解,并mldivide找到最小化均方误差的矩阵.在这种情况下,同样的操作A\X(:,n)(方法2)可以一次就完成了所有n与A\X(方法3).

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

对于500*500(或2.5E5)元素的情况,从加速Method 1到Method 3大约是x3500!

这也是有趣的观察输出的profile(在这里,为500*500的情况下):

方法1

方法2

方法3

从上面可以看出,通过重新排列元素squeeze并flipud占据运行时的大约一半(!)Method 2.还可以看出,溶液从细胞转化为基质的过程中损失了一些时间.

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

笔记:

"压缩"和"显式"版本之间的差异很小Method 3,有利于"显式"版本.因此,它不包括在比较中.

尝试了一种解决方案,其中输入Method 3为gpuArray-ed.这并没有提供改进的性能(甚至有些降低它们),可能是由于错误的实现,或者与在RAM和VRAM之间来回复制矩阵相关的开销.