Vectorization of solving a system of linear equations in MATLAB

Question

Summary: This question is about improving the linear regression calculation algorithm.

---

I have an array of 3D ( dlMAT) representing monochrome photographs of the same scene, performed with different exposure times (vector IT). Mathematically, each vector along the 3rd dimension dlMATrepresents a separate linear regression problem, which must be solved. The equation, the coefficients of which must be estimated, has the form:

DL = R*IT^Pwhere DLand ITobtained experimentally Rand Pshould be evaluated.

The above equation can be converted to a simple linear model after applying the logarithm:

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

Below is the most “naive” way to solve this system of equations, which essentially involves iterating over all the “third dimension vectors” and setting the order polynomial 1to (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

The above approach seems to be a good candidate for vectorization, because it does not use the core power of MATLAB, which is a MATrix operation. For this reason, it does not scale very well and takes much longer than I think it should.

---

There are alternative ways to do this calculation based on matrix division, as shown here and here , which include something like this:

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

That is, adding a vector to the observations 1, and then mldivide, to solve the system of equations.

, , - ( ).

№ 1:. , , (, , , )?

№2 (): ?

Answer 1

, , - . ( ), , \ (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- , ( script ) - , .

:

• "" "" Method 3 "" . .
• , Method 3 gpuArray -ed. ( ), , - , , VRAM.