计算 O(n) 中的最大归一化平均子数组

Question

这是最大和子数组问题的一种变体，但有一点变化（也类似于大小为 k 的最大平均子数组）。我不想找到最大总和，而是想找到归一化平均值，即sum/sqrt(length)。

例如，如果这是我的数组 ，c(1, -3, -4, 2, 6, 8, -4, -3, -5, 1, 7, -9, 3, 2)则最大和子数组问题将输出 max_sum = 16，并且子数组将为c(2, 6, 8)，您可以使用 Kadane 算法在线性时间内实现此目的。就我而言，输出为9.899495，子数组为c(6, 8)。

这显然可以很容易地在 中完成O(n^2)，但我正在尝试看看是否可以在 中完成O(n)。我尝试提出 Kadane 算法的变体，但它并不完全准确（与O(n^2)算法相比）。我使用过r，但欢迎您使用任何语言或伪代码，我的模拟如下：

set.seed(1) #I tried from 1 to 20 for comparison
n=30000
z <- rnorm(n)
normalized_mean_subarray(z)

理想情况下，输出应该具有最大值、子数组开始和子数组结束，但单独的最大值也一样大！提前致谢！

编辑

我将提供一个 O(n^2) 方法作为参考。对于n=10000，它会运行大约1.5分钟，所以检查结果是可以的。

你们中的一些人要求提供一个带有负值的示例。当n=2000和 时 set.seed(1)，max_interval 为 4.171375，从 1270 到 1295。如果检查该间隔 ( ) 中的值print(z[1270:1295])，您将能够看到许多负值。

set.seed(1)
n <- 2000
z <- rnorm(n)

Z <- function(I) {
  (length(I)**-0.5) * sum(z[I])
}

max_interval <- 0

start_time <- Sys.time()

for (window in 1:(n-1)) {
  normalizer <- (window**-0.5)
  base <- Z(1:window)
  max_interval_window <- 0

  for (i in 2:(n-window)){
    new_interval <- base + normalizer*(z[i+window-1] - z[i-1])
    max_interval_window <- ifelse(new_interval > max_interval_window, new_interval, max_interval_window)
    base <- new_interval
  }

  max_interval <- max(max_interval, max_interval_window)
}

end_time <- Sys.time()
print(end_time - start_time)

max_interval

这只是实现 O(n^2) 的一种方法，而且可能速度较慢，目标始终是在线性时间内完成此操作（希望是可行的）。

Answer 1

这是一个精确的解决方案，根据经验，在随机高斯输入上速度很快（O(n log n)），如果我付出更多努力，可能会证明如此。在我的笔记本电脑上，处理一百万个元素通常需要不到五秒的时间。（另一方面，最坏的情况是 O(n\xc2\xb2 log n)，我隐隐怀疑 I\xe2\x80\x99m 描述的算法并不是实例化基本思想的最佳方法。）

\n

这是一种分而治之的算法。每个征服步骤都会寻找穿过阵列中间的最佳解决方案。如果我们有一个快速的 (max, +) 卷积，我们可以应用左侧累积和和右侧累积和来确定每个长度的最佳解决方案，但没有人知道如何比 (max, +) 卷积更快地进行 (max, +) 卷积最坏情况下的二次时间。然而，我们处理的不是最坏情况的数据。由于随机高斯 n 元素向量的累积和的累积最大向量通常具有大约 \xe2\x88\x9 个不同元素，因此我们可以丢弃其他元素并用强力进行卷积。

\n

在 R 中（不是我选择的武器，但也许我\xe2\x80\x99m 会变得更好）：

\n

# From https://stackoverflow.com/a/76143942/2144669 by ThomasIsCoding\nf1 <- function(x) {\n    v <- c(0, cumsum(x))\n    msk <- lower.tri(matrix(NA, length(v), length(v)))\n    p <- outer(v, v, `-`)[msk]\n    q <- sqrt(abs(row(msk) - col(msk)))[msk]\n    max(p/q)\n}\n\n\nf2 <- function(x) {\n    if (length(x) > 50) {\n        mid <- length(x)%/%2\n        x_left <- x[1:mid]\n        x_right <- x[(mid + 1):length(x)]\n        max(f2(x_left), f2(x_right), max_crossover(x_left, x_right))\n    } else {\n        f1(x)\n    }\n}\n\nmax_crossover <- function(x_left, x_right) {\n    right <- hull(x_right)\n    best <- -Inf\n    for (left in hull(rev(x_left))) {\n        combos <- left + right\n        best <- max(best, Re(combos)/sqrt(Im(combos)))\n    }\n    best\n}\n\nhull <- function(x) {\n    cumsums <- cumsum(x)\n    is_cummax <- cumsums == cummax(cumsums)\n    complex(real = cumsums[is_cummax], imaginary = (1:length(x))[is_cummax])\n}\n\n\nf2(rnorm(1e+06))\n

\n