Vectorize matrix multiplication

Question

I have a school project about running a SOR algorithm on Octave, but mine is very inefficient. So I have this snippet of code:

for ii=1:n
    r = 1/A(ii,ii);
    for jj=1:n
        if (ii!=jj)
            A(ii,jj) = A(ii,jj)*r;
        end;
    end;
    b(ii,1) = b(ii,1)*r;
    x(ii,1) = b(ii,1);
end;

How can I vectorize this? My first attempt was this:

for ii=1:n
    r = 1/A(ii,ii);
    A(find(eye(length(A))!=1)) = A(find(eye(length(A))!=1))*r;
    b(ii,1) = b(ii,1)*r;
    x(ii,1) = b(ii,1);
end;

But I'm not sure it helped much. Is there a better and/or more efficient way to do it?

Thanks!

Answer 1

You can totally avoid loops I believe. You have to see that you are taking as the reciprocal of diagonal elements of A, then why do you want to use a loop. Do it directly. That's the first step. Remove Inf and now you want to multiply r to corresponding non-diagonal elements of corresponding rows, right?

Therefore, use repmat to construct such a matrix which will have replicated elements along columns, since you are multiplying same r to (1,2), (1,3), ... , (1,n). But R has non-zero diagonal elements. Therefore make them zero. Now you will get your A except that the diagonal elements will be zero. Therefore, you just have to add them back from original A. That can be done by A=A.*R+A.*eye(size(A,1)).

Vectorization comes from experience and most importantly analyzing your code. Think at each step whether you want to use the loop, if not replace that step with the equivalent command, other code will follow (for example, I constructed a matrix R, whereas you were constructing individual elements r. So I only thought about converting r -> R and then the rest of the code will just fall into place).

The code is as follows:

R=1./(A.*eye(size(A,1))); %assuming matrix A is square and it does not contain 0 on the main diagonal
R=R(~isinf(R));
R=R(:);
R1=R;
R=repmat(R,[1,size(A,2)]);
R=R.*(true(size(A,1))-eye(size(A,1)));
A=A.*R+A.*eye(size(A,1)); %code verified till here since A comes out to be the same

b = b.*R1;
x=b;

Answer 2

I suppose there are matrices:

A (NxN)
b (Nx1)

The code:

d = diag(A);
A = diag(1 ./ d) * A + diag(d - 1);
b = b ./ d;
x = b;

Answer 3

Randomly bumped onto this problem and at first glance looked interesting given the fact that the problem was tagged as a vectorization problem.

I was able to come up with a bsxfun based vectorized solution that also uses diagonal indexing. This solution seems to give me a 3-4x speedup over the loop code with decent sized inputs.

Assuming that you are still somewhat interested in seeing speedup improvement on this problem, I would be eager to know the kind of speedups you would be getting with it. Here's the code -

diag_ind = 1:size(A,1)+1:numel(A);
diag_A = A(diag_ind(:));
A = bsxfun(@rdivide,A,diag_A);
A(diag_ind) = diag_A;
b(:,1) = b(:,1)./diag_A;
x(:,1) = b(:,1);

Let me know!