基于边际均值求解满足覆盖率约束的最优行列子集方法

• 针对纵向研究数据集构建的二进制可用矩阵df：行对应观测个案，列对应测量时间点，取值为1代表该位置数据完整可用。需求为筛选最优行列子集，要求子集内每一行、每一列的1值占比至少达到2/3。
• 核心难点：行筛选与列筛选存在双向依赖——剔除不符合占比要求的行后，满足阈值的列集合会变动，反之亦然。该问题属于整数规划范畴，但之前未找到可落地的实现方案；最初设想的「占比≥0.66约束下同时最大化rowMeans和colMeans」思路未经过验证，不确定是否能得到全局最优解。
• 已尝试的暴力枚举方案基于combn()遍历所有行、列组合，在5行10列的小样本上可以得到3行6列的正确结果，但存在两个致命缺陷：
  • 计算量随矩阵规模指数增长，71行×155列的实际数据集完全无法运行
  • 将行、列筛选拆分为两个独立步骤，无法实现联合优化，得到的结果是次优的

该问题是典型的0-1整数规划问题，不需要暴力枚举，用开源整数规划求解器即可在秒级处理当前规模的数据，核心逻辑如下：

1. 定义两个二进制决策变量：
   • r[i]：取1代表保留第i行，取0代表剔除
   • c[j]：取1代表保留第j列，取0代表剔除
2. 优化目标：优先最大化保留的行、列总数量（可根据研究偏好调整权重，比如需要更多观测就给行变量加更高权重），可附加极小的稠密性权重，优先选择同等规模下完整率更高的子集。
3. 约束条件用大M法实现：
   • 对所有行i：如果行i被保留，则该行在保留列中的1值占比≥2/3；如果行被剔除，该约束自动失效
   • 对所有列j：如果列j被保留，则该列在保留行中的1值占比≥2/3；如果列被剔除，该约束自动失效
     该建模方式是行列联合优化，求解得到的是全局最优解，不存在分步筛选的次优问题，且计算复杂度随矩阵规模线性增长，完全适配71×155的数据集。

依赖lpSolve开源求解器，无需安装商业软件：

library(lpSolve)

select_optimal_subset <- function(df, thresh = 2/3) {
  n_row <- nrow(df)
  n_col <- ncol(df)
  M <- max(n_row, n_col) + 1 # 大M法常数，取比最大维度大1即可避免浮点误差
  n_vars <- n_row + n_col # 决策变量：前n_row个为行选择变量，后n_col个为列选择变量

  # 目标函数：最大化保留的行列总数，可按需调整行/列权重，比如c(rep(2,n_row), rep(1,n_col))即优先保留更多行
  obj <- c(rep(1, n_row), rep(1, n_col))

  const_list <- list()
  const_rhs <- c()
  const_dir <- c()

  # 构建行约束：保留行的1值占比≥阈值
  for (i in 1:n_row) {
    coef <- rep(0, n_vars)
    coef[i] <- -M
    for (j in 1:n_col) {
      coef[n_row + j] <- df[i, j] - thresh
    }
    const_list <- c(const_list, list(coef))
    const_rhs <- c(const_rhs, -M)
    const_dir <- c(const_dir, ">=")
  }

  # 构建列约束：保留列的1值占比≥阈值
  for (j in 1:n_col) {
    coef <- rep(0, n_vars)
    coef[n_row + j] <- -M
    for (i in 1:n_row) {
      coef[i] <- df[i, j] - thresh
    }
    const_list <- c(const_list, list(coef))
    const_rhs <- c(const_rhs, -M)
    const_dir <- c(const_dir, ">=")
  }

  const_mat <- do.call(rbind, const_list)

  # 求解0-1整数规划
  res <- lp(
    direction = "max",
    objective.in = obj,
    const.mat = const_mat,
    const.dir = const_dir,
    const.rhs = const_rhs,
    all.bin = TRUE
  )

  if (res$status != 0) stop("未找到可行解，请检查阈值设置是否合理")

  # 解析结果
  r_selected <- which(res$solution[1:n_row] >= 0.99)
  c_selected <- which(res$solution[(n_row+1):n_vars] >= 0.99)
  sub_mat <- df[r_selected, c_selected, drop = FALSE]

  # 返回结果与校验指标
  list(
    selected_row_index = r_selected,
    selected_col_index = c_selected,
    subset_dim = dim(sub_mat),
    row_complete_ratio = round(rowMeans(sub_mat), 3),
    col_complete_ratio = round(colMeans(sub_mat), 3)
  )
}

# 测试（5行10列示例，和之前暴力枚举的测试场景一致）
set.seed(123)
test_mat <- matrix(sample(0:1, 50, replace = T, prob = c(0.3, 0.7)), nrow = 5)
test_res <- select_optimal_subset(test_mat)
print(test_res$subset_dim) # 输出保留的行列数

• 71行×155列的数据集总决策变量仅226个，约束共226条，普通消费级笔记本上数秒即可得到结果，无指数爆炸问题
• 如果需要严格按整数计数（比如3列至少2个1、4列至少3个1，避免浮点比例的计数误差），只需要把约束中的比例形式改写为整数形式即可
• 如果对保留方向有偏好，比如希望尽可能多保留测量波次、尽可能多保留有效样本，直接调整目标函数obj中对应变量的权重即可，权重越高求解时越优先保留对应维度

内容的提问来源于stack exchange，提问作者JannisOverflows