Matlab: How to vectorize a nested loop over a 2D set of vectors

Question

I have a function in the following form:

function Out = DecideIfAPixelIsWithinAnEllipsoidalClass(pixel,means,VarianceCovarianceMatrix)  
   ellipsoid = (pixel-means)'*(VarianceCovarianceMatrix^(-1))*(pixel-means);  
   if ellipsoid <= 1
      Out = 1;
   else
      Out = 0;
   end
end

I am doing remote-sensing processes with matlab and I want to classify a LandSatTM images.This picture has 7 bands and is 2048*2048.So I stored them in 3 dimentinal 2048*2048*7 matrix.in this function means is a 7*1 matrix calculated earlier using the sample of the class in a function named ExtractStatisticalParameters and VarianceCovarianceMatrix is a 7*7 matrix in fact you see that:

ellipsoid = (pixel-means)'*(VarianceCovarianceMatrix^(-1))*(pixel-means);

is the equation of an ellipsoid.My problem is that each time you can pass a single pixel(it is a 7*1 vector where each row is the value of the pixel in a seperated band) to this function so I need to write a loop like this:

for k1=1:2048  
   for k2=1:2048  
      pixel(:,1)=image(k1,k2,:); 
      Out = DecideIfAPixelIsWithinAnEllipsoidalClass(pixel,means,VarianceCovarianceMatrix);  
   end  
end

and you know it will take alot of time and energy of the system.Can you suggest me a way to reduce the pressure applied on the systam?

score 6 · Accepted Answer · edited May 23 '17 at 12:01

6

No need for loops!

pMinusMean = bsxfun( @minus, reshape( image, [], 7 ), means' ); %//' subtract means from all pixes
iCv = inv( arianceCovarianceMatrix );
ell = sum( (pMinusMean * iCv ) .* pminusMean, 2 ); % note the .* the second time!
Out = reshape( ell <= 1, size(image(:,:,1)) ); % out is 2048-by-2048 logical image

Update:

After a (somewhat heated) debate in the comments below I add a correction made by Rody Oldenhuis:

pMinusMean = bsxfun( @minus, reshape( image, [], 7 ), means' ); %//' subtract means from all pixes
ell = sum( (pMinusMean / varianceCovarianceMatrix ) .* pminusMean, 2 ); % note the .* the second time!
Out = reshape( ell <= 1, size(image(:,:,1)) );

The key issue in this change is that Matlab's inv() is poorly implemented and it is best to use mldivide and mrdivide (operators / and \) instead.

edited May 23 '17 at 12:01

Community

1
1

answered Jun 17 '13 at 20:48

Shai

111,146
38
238
371

1

...did you just use `inv()` to solve a linear system? Near the top of `doc inv`: "In practice, it is seldom necessary to form the explicit inverse of a matrix. A frequent misuse of `inv` arises when solving the system of linear equations `Ax = b`. One way to solve this is with `x = inv(A)*b`. A better way, from both an execution time and numerical accuracy standpoint, is to use the matrix division operator `x = A\b`. This produces the solution using Gaussian elimination, without forming the inverse. See `mldivide` (`\`) for further information." – Rody Oldenhuis Jun 17 '13 at 21:15
@RodyOldenhuis `inv` is NOT used for solving systems of equations in this question, it is used as the covariance matrix of a gaussian (ellipsoid). – Shai Jun 17 '13 at 21:22
euhmmm...so? What prevents you from writing `sum( pMinusMean/varianceCovarianceMatrix .* pminusMean, 2)`? – Rody Oldenhuis Jun 17 '13 at 21:31
Oh well, in that case, let's just use `pi = 3` from now on shall we, that's far more readable than writing all those digits all the time... Who cares about all those digits anyway, right? Suggesting someone use `inv()` is teaching that person bad habits; you really ***NEVER*** need or should use `inv()`. – Rody Oldenhuis Jun 17 '13 at 21:38
[Never](http://blogs.mathworks.com/loren/2007/05/16/purpose-of-inv/)? Let's not go crazy. – horchler Jun 17 '13 at 22:09
@Shai - ah yes, the `bsxfun` you love so much :) – Dang Khoa Jun 18 '13 at 02:42
@DangKhoa I run `bsxfun` in a loop, at the background - just for fun – Shai Jun 18 '13 at 05:03
@Shai that wouldn't be very efficient then, would it? :) – Eitan T Jun 18 '13 at 06:20
@horchler: No, really, ***NEVER***. Those four reasons Loren mentions are easily nullified: 1. teaching proper use of `mrdivide`, `mldivide`, `lu()`, `qr()` etc. (and *why*) is actually much more valuable and durable. 2. you could just deprecate `inv()`, or 3. you could make `inv() ` nothing more than a wrapper around `mrdivide`/`mldivide`. 4. just use `A\eye(size(A))` if you **really** insist. But for this particular question, `iCv` is used for no other purpose than to solve a post-multiplied system (`X*A = B`), so why on Earth wouldn't you use the best tool available for that (`mrdivide`) ? – Rody Oldenhuis Jun 18 '13 at 14:06
@RodyOldenhuis - I updated my answer to reflect your argumnets. – Shai Jun 18 '13 at 14:15
2

@Shai: Well that wasn't really my intention, but thanks :) It's not that it's "poorly implemented", it's just that no algorithms exist that are as fast and accurate as the ones you can use to compute the products `inv(A)*x` or `x*inv(A)`. Sorry if I come off a bit strong about this issue, but when you see people use `i` or `j` as a variable name, what would you do? :) I'm saying that it's nice to implement a bubble-sort or bogosort for educational reasons once or twice, but really, in the end what you should *use* and *teach* is quicksort (and similar) -- same with `inv()`. – Rody Oldenhuis Jun 18 '13 at 14:41
1

You comment wasn't about `A*X = B`, but rather that one should NEVER use `inv`. I don't don't like blanket statements. As you say, teaching new users why things work one way or another and what the trade-offs are is better than just saying "NEVER." As is mentioned, `inv` is a simple function that has speed advantages in some cases. Precision down to `eps` is not always desired as long as one understands the numerics of a particular problem. "Best" does not always mean most accurate. `trace(A\eye(size(A)))` takes 50% longer than `trace(inv(A))` for large random matrices on my machine. – horchler Jun 18 '13 at 14:49
@horchler: name 1 case, any case, where you think you *need* `inv()`, and mention why, and I'll come up with a better alternative. – Rody Oldenhuis Jun 18 '13 at 19:18
@RodyOldenhuis Suppose you have a 3x3 affine transformation relating image I to image J such that `imtransform(I, H) == J`. How would you find the inverse transform relating `J` to `I`? – Shai Jun 18 '13 at 19:35
@Shai: Do you mean something like [this](http://negativeprobability.blogspot.de/2011/11/affine-transformations-and-their.html)? Sure looks like an ordinary linear system to me there... – Rody Oldenhuis Jun 18 '13 at 19:57

Matlab: How to vectorize a nested loop over a 2D set of vectors

1 Answers1

Update:

Linked