Is it possible to replace ranges of a matrix elements using indices from two vectors?

Hi,

I have a matrix of size 5x5 all the elements of which are zeros. I would like to replace some of these zeros with ones. The elements that I want to replace are referenced by two vectors a & b each of size (5 x 1). These two vectors point to the beginning and ending column indices of the range to be replaced. These ranges vary from one row to another. I tried using X(:,a:b)=1, but it only uses the first value in these vectors. I don't know how to modify my code so that the ranges vary by each row.

> X=zeros(5,5)
X =
     0     0     0     0     0
     0     0     0     0     0
     0     0     0     0     0
     0     0     0     0     0
     0     0     0     0     0
> a=[2,3,1,3,4]'
a =
     2
     3
     1
     3
     4
> b=[4,3,2,4,5]'
b =
     4
     3
     2
     4
     5
> X(:,a:b)=1
X =
     0     1     1     1     0
     0     1     1     1     0
     0     1     1     1     0
     0     1     1     1     0
     0     1     1     1     0
How can I get the result look like this?

X = 
        0    1  1  1  0
  0  0  1  0  0
  1  1  0  0  0
  0  0  1  1  0
  0  0  0  1  1

Thanks! Hawre

Answer 1

[Matt Fig]

Modified FOR loop for speed...

    N = 5;
    X = false(N);
    for ii = 1:N
        X(ii,a(ii):b(ii)) = 1;
    end
    X = double(X);  % Convert to double.

Note that if logicals will work for your application, then the last line is not needed.

%

%

%

EDIT In response to your extra request for a vectorized solution.

X2 = zeros(N,N);
C = b-a+1;
R = zeros(1,sum(C));
R([0;cumsum(C(1:N-1))]+1) = 1;
C = arrayfun(@colon,a,b,'un',0);  % Or try MCOLON
C = [C{:}];  % Not necessary with MCOLON
X2(cumsum(R)+(C-1)*N) = 1;

MCOLON might speed up the call to ARRAYfUN. Using ARRAYFUN as above, this is much slower than the FOR loop for 5000-by-5000 X (N=5000), and occasionally runs out of memory. I tested the above using:

N = 5000;
a = ceil(rand(N,1)*N);
b = min(a+round(rand(N,1)*N),N);

[Hawre Jalal]

Thanks! But is it possible without using loops?

[Matt Fig]

Yes, but I doubt faster. Unless your array is much larger...

[Hawre Jalal]

Would matrix programming be faster for larger matrices for this problem? I need to use this matrix in some simulations that involves millions of observations. I was worried that using a loop will hinder the performance and was wondering if there is a faster way of doing it. Thanks again for your prompt response.

Answer 2

[Andrei Bobrov]

more variant

indcol = ones(5,1)*(1:5);
X = +(a(:,ones(5,1))<=indcol & b(:,ones(5,1))>=indcol)

and

C = arrayfun(@(x,y,z)[z*ones(y-x+1,1), (x:y).'],a,b,(1:length(a))','un',0);
C = cat(1,C{:});
X(sub2ind(size(X),C(:,1),C(:,2))) = 1;

[Oleg Komarov]

+ 1 for the first, also bsxfun version:

indcol = 1:5;

X = bsxfun(@le,a,indcol) .* bsxfun(@ge,b,indcol);

P.S: I think + is obscure as compared to double.

[Oleg Komarov]

Timings of Matt's LOOP and bsxfun LOGICAL (more or less the same):

N = 5000;

a = ceil(rand(N,1)*N);

b = min(a+round(rand(N,1)*N),N);

times = zeros(100,2);

for n = 1:100

tic

indcol = 1:N;

X1 = bsxfun(@le,a,indcol) .* bsxfun(@ge,b,indcol);

times(n,1) = toc;

end

for n = 1:100

tic

X2 = zeros(N);

for ii = 1:N

X2(ii,a(ii):b(ii)) = 1;

end

times(n,2) = toc;

end

isequal(X1,X2)

mean(times(2:end,:))

[Matt Fig]

When I put this into a function, the FOR loop is more than 2.5 times as fast (using a modified loop from my original post...). Did you have your timings in a function or a script, Oleg? I am using 2007b also, so something may have changed....

function [] = address_vects()

N = 4000;

T = [0 0];

for jj = 1:100

a = ceil(rand(N,1)*N);

b = min(a+round(rand(N,1)*N),N);

tic

X = false(N,N);

for ii = 1:size(X,1)

X(ii,a(ii):b(ii)) = 1;

end

X = double(X); % Convert to double.

T(1) = T(1) + toc;

tic

% indcol = ones(N,1)*(1:N);

% X1 = +(a(:,ones(N,1))<=indcol & b(:,ones(N,1))>=indcol);

indcol = 1:N;

X1 = bsxfun(@le,a,indcol) .* bsxfun(@ge,b,indcol);

T(2) = T(2) + toc;

end

reltime = T/min(T)

>> address_vects

reltime =

1 2.6987

[Oleg Komarov]

Vista 32b r2010b. It may be multithreading of bsxfun.

reltime =

1.0497 1.0000

[Matt Fig]

O.k., that is what I thought. Still, 5% is pretty close. I wish BSXFUN had built-in string inputs (think CELLFUN). Then I bet this solution would be MUCh faster.

[Hawre Jalal]

Matt and Oleg, Thank you very much for your input. It helped clarify a bit more of how Matlab handles matrices and vectors.