r - 查找最大和最小矩阵边际总变异性

Question

是否有一种更优雅的方法可以根据二进制矩阵的填充和大小来计算其最大或最小变化率（CV）？考虑到所有行和列的总数必须不为零。例如

foo(n_col, n_row, fill){ get maximum possible CV }

假设我们有一个名为 m的矩阵，其中所有列和行的总数均为 > 0，但矩阵的填充量最少。

m <- matrix(rep(0,25), nrow = 5)
diag(m) <- 1
#     [,1] [,2] [,3] [,4] [,5]
#[1,]    1    0    0    0    0
#[2,]    0    1    0    0    0
#[3,]    0    0    1    0    0
#[4,]    0    0    0    1    0
#[5,]    0    0    0    0    1

variability1 <- sd(colSums(m))/mean(colSums(m))
variability1
# [1] 0
# the maximum and minimum for this fill is zero 
# considering that  all column and row totals must be > 0

也许我们可以像增加填充量一样检查最大值：

# find out which matrix elements are zeros
empty <- which(m < 1)
# vector for results
variability <- rep(NA, length(empty))
#
for(i in 1:length(variability)){
m[empty[[i]] ] <- 1
variability[[i]] <- sd(colSums(m))/mean(colSums(m))
}
# we get what should the maximum CV for each given level of matrix fill...
c(variability1, variability)

我认为像这样以矩阵方式填充矩阵会在边际列总计中保持最大的可变性吗？是否有一种更简单的方法来处理不同大小，填充和形状的矩阵的最大和最小可变性？

Answer 1

下面提供了问题的另一种表示形式，它是对二进制矩阵的列总和向量的选择进行优化，从而使给定fill的可变性最大化。提供了有关该公式有效性的非正式论据以及解决该问题的算法。生成的算法与OP的断言一致


像这样逐列填充矩阵，以保持边际列总计的最大可变性


首先，通过fill二进制矩阵1将n_row定义为n_col中m的数目。根据问题陈述的约束，m是二进制矩阵，所有行和列之和都大于零，所以fill是范围[max(n_row, n_col),n_row*n_col]的整数。

然后问题是对于给定值fill在范围[max(n_row, n_col),n_row*n_col]中，找到最大值

 sd(colSums(m))/mean(colSums(m))

在所有 m上，使得 m是具有 fill个数 1且所有行列总和都大于零的二进制矩阵。

我们注意到最好用 m的列总和的向量而不是 m本身来指定此优化问题的范围。这是因为存在不同的 m，它们具有相同的列和向量，因此具有相同的目标值。将列总和的向量表示为 x，可以将上述优化问题重申为最大化之一：

sd(x)/mean(x)

这样， x的每个元素都是 [1, n_row]范围内的整数，而 sum(x)是 fill。

此外，由于 sum(x)被约束为等于 fill，因此对于给定的 mean(x)，分母项 x在所有 fill上都是恒定的。因此，最大化的等效目标函数就是 sd(x)或等效地 x的方差。

为了最大化 x的方差，我们需要选择 x，以使其值之间的差异最大，同时仍然满足 x的约束。在这里，我们可以相对于 fill归纳地考虑这个问题。假设对于给定的 fill，我们有 x的解决方案，可以在满足其约束的同时最大化 x的方差。问题就变成了：当我们将 fill增加到 fill + 1时，新的 x使它的方差最大化的是什么？因为我们有 sum(x)=fill并且 x中的每个元素都是整数的约束，所以递增 fill意味着我们必须递增 x的一个且仅一个元素。现在放宽对 x中每个元素的上限约束（即 x[i] <= n_row中所有 i的 [1,n_col]），然后问题就变成了： x中哪个元素要增加，从而最大程度地增加了。 x的方差。对于这个问题的答案，我们可以看一下 var(x)的泰勒级数展开式：

var(x + dx) = var(x) + gradient(var(x)) %*% dx + 1/2 * t(dx) %*% Hessian(var(x)) %*% dx

其中 dx是长度为 n_col的向量，其中一个元素等于 1，所有其他元素为 0（即指示符向量）。由于 var(x)在 x中是平方的，因此二阶扩展就足够了。此外，由于 dx是指示符向量，因此仅Hessian矩阵的对角线元素很重要。这些是由：

gradient(var(x))[i] = 2*(x[i]-mean(x))/(n_col-1),      for all i in [1,n_col]
Hessian(var(x))[i,i] = 2/n_col                  ,      for all i in [1,n_col]

由于Hessian的所有对角项均相同，因此对于 dx的任何选择，泰勒级数的二阶项均相同。因此，只有一阶项在确定 x中要递增哪个元素以使 x的方差增加最大化时很重要。从梯度项来看，很明显，我们应该选择增加 i中具有最大当前值 x的 x[i]个元素，以最大程度地增加 x的方差。现在，我们在 x的每个元素上重新引入上限约束。然后，最佳选择是增加 i中具有最大当前值 x的第 x[i] < n_row个元素。请注意，如果 x中有多个具有相同最大值 x[i] < n_row的此类元素，则选择这些元素中的任何一个都会导致 x的方差具有相同的最大增加。

到目前为止，我们已经显示的是，给定一个 fill和 x的解决方案，它在满足其约束的同时最大化 x的方差，我们有一条规则 dx将最大化 x的增量增量cc>表示 fill + 1。仍然需要说明的是，该规则会产生一个新的 x，它是最优的 x，它使新 x的 fill + 1的方差最大。现在，我们通过矛盾来证明这一点。具体来说，如果这个新的 x不能使 x的 fill + 1方差最大化，则必须存在另一个针对 x_1的列和 fill的向量，以及另一个规则 dx_1

var(x_1 + dx_1) > var(x + dx)

但是，由于 x使 var(x)最大化 fill，并且梯度和Hessian方程对于任何 x成立，因此我们具有：

var(x_1 + dx_1) = var(x_1) + gradient(var(x_1)) %*% dx_1 + 1/2 * t(dx_1) %*% Hessian(var(x_1)) %*% dx_1
                <= var(x_1) + 2*(max(x_1)-mean(x_1))/(n_col-1) + constant
                <= var(x) + 2*(max(x)-mean(x))/(n_col-1) + constant
                = var(x + dx)

因此矛盾。为了更清楚地解释步骤：


从步骤1到步骤2，在 dx_1处 x_1的最佳选择是使 x_1中的最大元素递增，从而使 gradient(var(x_1)) %*% dx_1 <= 2*(max(x_1)-mean(x_1))/(n_col-1)递增的选择。同样，对于所有 x和 dx对于给定的 fill，二阶项都是常数，因此我们简单地将其声明为 constant。
从第2步到第3步，我们假设（cc） var(x_1) <= var(x)使 x的方差最大，（ii）对于 fill的最优规则 gradient(var(x)) %*% dx = 2*(max(x)-mean(x))/(n_col-1) dx，并且（ iii）假定 fill使 max(x_1) <= max(x)处的方差最大。要查看后者，请考虑 x的列总和 fill的公共前导向量。 x_-1递增到 fill-1和 x_-1之间的规则差异仅仅是 x中要递增的元素的不同选择。从 x_1方差的泰勒级数展开中可以明显看出， x_-1中要递增到 x_-1的元素的选择必须大于或等于要去 x_-1的选择，因为< cc>。因此， x。现在，将此推理扩展到任何先前的 x_1处列总和的任何公共向量，包括初始 var(x_1) <= var(x)，其中列总和的初始向量为所有 max(x_1) <= max(x)。然后，为一条路径选择上面定义的最佳增量规则以达到 fill-k >= max(n_row, n_col)；而对于其他路径，请选择一条任意的增量规则路径以达到 fill。由于最佳规则总是在每个步骤中递增状态的最大元素（以上限为准），因此很显然再次 1。


最后，为了完成数学归纳，我们从初始填充开始，其中 x是所有 x_1。由于给定了初始填充，对于 max(x_1) <= max(x)没有其他选择，因此这对 x进行了优化。现在，最佳增量规则 1是选择 var(x)的第一个元素进行增量，因为所有元素都相等。由于增加 x的任何其他元素将导致相同的方差，因此所得的 dx会使初始填充的方差最小，并且加1。

上面的论点立即建议使用以下算法在列和的向量上分配 x值：


循环浏览列总和 x + dx的向量中的每个元素。
对于每个 x元素 fill。请注意，我们从填充中减去 x，以便保留它们以至少用 i填充列求和向量的其余元素，并将数量限制为 x[i] <- min(n_row, fill - (ncol_-i))以满足问题的约束。
更新 (n_col-i)


该算法和相关的参数验证了OP的断言：


像这样逐列填充矩阵，以保持边际列总计的最大可变性


在R中，代码如下所示：

foo <- function(n_col, n_row, fill) {
  ## preallocate the vector of column sums x and initialize to NA
  x <- rep(NA, n_col)
  for (i in seq_len(n_col)) {
    x[i] <- pmin.int(n_row, fill-(n_col-i))
    fill <- fill - x[i]
  }
  ## compute the variability given the vector of column sums x
  sd(x)/mean(x)
}

认识到循环中 1的重复减量可以用 n_row代替，以上简化为：

foo <- function(n_col, n_row, fill) {
  x <- pmin.int(pmax.int(cumsum(c(fill-n_col+1,rep(-n_row+1,n_col-1))),1),n_row)
  ## compute the variability given the vector of column sums x
  sd(x)/mean(x)
}

使用此功能，我们可以恢复OP的结果：

n_col=5
n_row=5
variability <- sapply(max(n_col,n_row):(n_col*n_row), function(fill) foo(n_col, n_row, fill))
print(variability)
## [1] 0.0000000 0.3726780 0.6388766 0.8385255 0.9938080 0.8660254 0.8131156 0.8122329 0.8426501
##[10] 0.7319251 0.6666667 0.6404344 0.6443795 0.5414886 0.4707512 0.4330127 0.4259177 0.3049184
##[19] 0.1944407 0.0931695 0.0000000