StackOverflow2013

Note that there are some explanatory texts on larger screens.

plurals

POOptimizing repetitive estimation (currently a loop) in MATLAB
text
Body
copied! <p>I've found myself needing to do a least-squares (or similar matrix-based operation) for every pixel in an image. Every pixel has a set of numbers associated with it, and so it can be arranged as a 3D matrix. </p> <p>(<strong>This next bit can be skipped</strong>)</p> <p>Quick explanation of what I mean by least-squares estimation : </p> <p>Let's say we have some quadratic system that is modeled by Y = Ax^2 + Bx + C and we're looking for those A,B,C coefficients. With a few samples (at least 3) of X and the corresponding Y, we can estimate them by:</p> <ol> <li>Arrange the (lets say 10) X samples into a matrix like <code>X = [x(:).^2 x(:) ones(10,1)];</code></li> <li>Arrange the Y samples into a similar matrix: <code>Y = y(:);</code></li> <li>Estimate the coefficients A,B,C by solving: <code>coeffs = (X'*X)^(-1)*X'*Y;</code></li> </ol> <p>Try this on your own if you want:</p> <pre><code>A = 5; B = 2; C = 1; x = 1:10; y = A*x(:).^2 + B*x(:) + C + .25*randn(10,1); % added some noise here X = [x(:).^2 x(:) ones(10,1)]; Y = y(:); coeffs = (X'*X)^-1*X'*Y coeffs = 5.0040 1.9818 0.9241 </code></pre> <p><strong>START PAYING ATTENTION AGAIN IF I LOST YOU THERE</strong></p> <p>*<em>MAJOR REWRITE</em>*I've modified to bring it as close to the real problem that I have and still make it a minimum working example.</p> <p><strong>Problem Setup</strong> </p> <pre><code>%// Setup xdim = 500; ydim = 500; ncoils = 8; nshots = 4; %// matrix size for each pixel is ncoils x nshots (an overdetermined system) %// each pixel has a matrix stored in the 3rd and 4rth dimensions regressor = randn(xdim,ydim, ncoils,nshots); regressand = randn(xdim, ydim,ncoils); </code></pre> <p>So my problem is that I have to do a (X'*X)^-1*X'*Y (least-squares or similar) operation for every pixel in an image. While that itself is vectorized/matrixized the only way that I have to do it for every pixel is in a for loop, like:</p> <p><strong>Original code style</strong> </p> <pre><code>%// Actual work tic estimate = zeros(xdim,ydim); for col=1:size(regressor,2) for row=1:size(regressor,1) X = squeeze(regressor(row,col,:,:)); Y = squeeze(regressand(row,col,:)); B = X\Y; % B = (X'*X)^(-1)*X'*Y; %// equivalently estimate(row,col) = B(1); end end toc Elapsed time = 27.6 seconds </code></pre> <p><strong>EDITS in reponse to comments and other ideas</strong><br> I tried some things:<br> 1. Reshaped into a long vector and removed the double <code>for</code> loop. This saved some time.<br> 2. Removed the <code>squeeze</code> (and in-line transposing) by <code>permute</code>-ing the picture before hand: This save alot more time.<br></p> <p><strong>Current example:</strong> <br></p> <pre><code>%// Actual work tic estimate2 = zeros(xdim*ydim,1); regressor_mod = permute(regressor,[3 4 1 2]); regressor_mod = reshape(regressor_mod,[ncoils,nshots,xdim*ydim]); regressand_mod = permute(regressand,[3 1 2]); regressand_mod = reshape(regressand_mod,[ncoils,xdim*ydim]); for ind=1:size(regressor_mod,3) % for every pixel X = regressor_mod(:,:,ind); Y = regressand_mod(:,ind); B = X\Y; estimate2(ind) = B(1); end estimate2 = reshape(estimate2,[xdim,ydim]); toc Elapsed time = 2.30 seconds (avg of 10) isequal(estimate2,estimate) == 1; </code></pre> <p><strong>Rody Oldenhuis's way</strong></p> <pre><code>N = xdim*ydim*ncoils; %// number of columns M = xdim*ydim*nshots; %// number of rows ii = repmat(reshape(1:N,[ncoils,xdim*ydim]),[nshots 1]); %//column indicies jj = repmat(1:M,[ncoils 1]); %//row indicies X = sparse(ii(:),jj(:),regressor_mod(:)); Y = regressand_mod(:); B = X\Y; B = reshape(B(1:nshots:end),[xdim ydim]); Elapsed time = 2.26 seconds (avg of 10) or 2.18 seconds (if you don't include the definition of N,M,ii,jj) </code></pre> <p><strong>SO THE QUESTION IS:</strong> <br> <strong>Is there an (even) faster way?</strong></p> <p>(I don't think so.)</p>

Querying!

Guidance

An individual column

Larger individual text columns get their own page to allow for proper reading.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload