强制矩阵的列在不同的范围内

Question

我有一个名为 l 的矩阵，大小为 20X3。 我想做的是： 假设我有这个限制：

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

我想强制矩阵 l 的所有元素在限制范围内。 l1_max 和 l1_min 中第一列的值。 l2_max 和 l2_min 中第二列的值。 mu_max & mu_min 内第 3 列的值。

我的做法是这样的：

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

可以做得更好吗？

Answer 1

您不必遍历行，对整列使用矢量化操作：

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

类似地：

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;

类似于 Bas 的想法的另一种方法是应用 min 和 max，如下所示：

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);

看来这两种方法的性能相当。

Answer 2

您甚至不必遍历所有列，对整个矩阵的操作可以在 2 次调用 bsxfun 中完成，与列数无关：

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

这使用了两个技巧：剪辑 min_val 和 max_val 之间的值，您可以使用clipped_x = min(max(x, min_val), max_val)。另一个技巧是使用有点晦涩的bsxfun，它在单例扩展之后应用一个函数。当您在两个矩阵上使用它时，它会在应用函数之前将最小的矩阵“挤压”到与最大矩阵相同的大小，因此上面的示例等效于 M = min(M, repmat(column_max, size(M, 1), 1))，但希望以更有效的方式计算。

Answer 3

以下是测试目前讨论的各种方法的基准。我正在使用 File Exchange 上的 TIMEIT 函数。

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

鉴于以下实现：

1) 遍历所有值并与最小值/最大值进行比较

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) 将每一列与最小值/最大值进行比较

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) 将每列截断到限制

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) (3) 中截断的矢量化版本

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

5) 绝对值比较：-a |x|

（注意：这不适用于您的情况，因为它需要一个对称的限制范围。我只是为了完整性而将其包括在内。此外，它原来是最慢的方法。）

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

---

我的机器上输入矩阵大小为 10000x500 的时间：

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

我想说上述所有方法都足够快，bsxfun 解决方案是最快的:)