forcing the columns of a matrix within different limits

Question

I have a matrix named l having size 20X3. What I wanted to do was this : Suppose I have this limits:

l1_max=20; l1_min=0.5;
l2_max=20; l2_min=0.5;
mu_max=20; mu_min=0.5;

I wanted to force all the elements of the matrix l within the limits. The values of 1st column within l1_max & l1_min. The values of 2nd column within l2_max & l2_min. The values of 3rd column within mu_max & mu_min.

What I did was like this:

for k=1:20
    if l(k,1)>l1_max 
        l(k,1) = l1_max;
    elseif l(k,1)<l1_min
        l(k,1) = l1_min;
    end

    if l(k,2)>l2_max 
        l(k,2) = l2_max;
    elseif l(k,2)<l2_min
        l(k,2) = l2_min;
    end

    if l(k,3)>mu_max 
        l(k,3) = mu_max;
    elseif l(k,3)<mu_min
        l(k,3) = mu_min;
    end
end

Can it be done in a better way ?

Answer 1

You don't have to loop over rows, use vectorized operations on entire columns:

l(l(:, 1) > l1_max, 1) = l1_max;
l(l(:, 1) < l1_min, 1) = l1_min;

Similarily:

l(l(:, 2) > l2_max, 2) = l2_max;
l(l(:, 2) < l2_min, 2) = l2_min;
l(l(:, 3) > l2_max, 3) = mu_max;
l(l(:, 3) < l2_min, 3) = mu_min;

An alternative method, which resembles to Bas' idea, is to apply min and max as follows:

l(:, 1) = max(min(l(:, 1), l1_max), l1_min);
l(:, 2) = max(min(l(:, 2), l2_max), l2_min);
l(:, 3) = max(min(l(:, 3), mu_max), mu_min);

It appears that both approaches have comparable performance.

Answer 2

You don't even have to loop over all columns, the operation on the whole matrix can be done in 2 calls to bsxfun, independent of the number of columns:

column_max = [l1_max, l2_max, mu_max];
column_min = [l1_min, l2_min, mu_min];

M = bsxfun(@min, M, column_max); %clip to maximum
M = bsxfun(@max, M, column_min); %clip to minimum

This uses two tricks: to clip a value between min_val and max_val, you can do clipped_x = min(max(x, min_val), max_val). The other trick is to use the somewhat obscure bsxfun, which applies a function after doing singleton expansion. When you use it on two matrices, it 'extrudes' the smallest one to the same size as the largest one before applying the function, so the example above is equivalent to M = min(M, repmat(column_max, size(M, 1), 1)), but hopefully calculated in a more efficient way.

Answer 3

Below is a benchmark to test the various methods discussed so far. I'm using the TIMEIT function found on the File Exchange.

function [t,v] = testClampColumns()
    % data and limits ranges for each column
    r = 10000; c = 500;
    M = randn(r,c);
    mn = -1.1 * ones(1,c);
    mx = +1.1 * ones(1,c);

    % functions
    f = { ...
        @() clamp1(M,mn,mx) ;
        @() clamp2(M,mn,mx) ;
        @() clamp3(M,mn,mx) ;
        @() clamp4(M,mn,mx) ;
        @() clamp5(M,mn,mx) ;
    };

    % timeit and check results
    t = cellfun(@timeit, f, 'UniformOutput',true);
    v = cellfun(@feval, f, 'UniformOutput',false);
    assert(isequal(v{:}))
end

Given the following implementations:

1) loop over all values and compare against min/max

function M = clamp1(M, mn, mx)
    for j=1:size(M,2)
        for i=1:size(M,1)
            if M(i,j) > mx(j)
                M(i,j) = mx(j);
            elseif M(i,j) < mn(j)
                M(i,j) = mn(j);
            end
        end
    end
end

2) compare each column against min/max

function M = clamp2(M, mn, mx)
    for j=1:size(M,2)
        M(M(:,j) < mn(j), j) = mn(j);
        M(M(:,j) > mx(j), j) = mx(j);
    end
end

3) truncate each columns to limits

function M = clamp3(M, mn, mx)
    for j=1:size(M,2)
        M(:,j) = min(max(M(:,j), mn(j)), mx(j));
    end
end

4) vectorized version of truncation in (3)

function M = clamp4(M, mn, mx)
    M = bsxfun(@min, bsxfun(@max, M, mn), mx);
end

5) absolute value comparison: -a < x < a <==> |x| < a

(Note: this is not applicable to your case, since it requires a symmetric limits range. I only included this for completeness. Besides it turns out to be the slowest method.)

function M = clamp5(M, mn, mx)
    assert(isequal(-mn,mx), 'Only works when -mn==mx')
    idx = bsxfun(@gt, abs(M), mx);
    v = bsxfun(@times, sign(M), mx);
    M(idx) = v(idx);
end

---

The timing I get on my machine with an input matrix of size 10000x500:

>> t = testClampColumns
t =
    0.2424
    0.1267
    0.0569
    0.0409
    0.2868

I would say that all the above methods are acceptably fast enough, with the bsxfun solution being the fastest :)