获取带有约束的线段上的点的位置

Question

我正在为 ZenUML 设计一个布局引擎。一项要求（简化后）是这样的：

1. 有一条线段；
2. 这条线段上有n个（n < 100，如果影响性能的话可以减少到n < 30）个点，顺序是固定的；(例如P1 ~ Pn)
3. 某些点之间存在已知的最小距离；（例如m_d（P2，P4）= 500）
4. 线段的长度应尽可能小；
5. （很好）相邻点之间的间隙应尽可能均匀（通过标准差^来衡量，并且不得违反1~4）。
6. （新增）在最坏的情况下（或多或少的限制），它必须在不超过 1 秒的时间内给出 30 分的结果。

^standard deviation被用作我可以用我有限的数学知识精确描述的最佳替代方案。然而，正如 Roman 和 David 指出的，在某些情况下，有更有吸引力的结果不能满足最低要求standard deviation。看到这个评论

有姐妹提问https://math.stackexchange.com/q/4377433。这与用数学语言（矩阵）描述的问题相同。

例如，对于有 4 个点（P1 ~ P4）的线段。

m_d（P1，P3）= 200，m_d（P1，P4）= 900。

• A。一种解决方案是 P1 = 0、P2 = 0、P3 = 200、P4 = 900。
• b. 更好的解决方案是 P1 = 0、P2 = 100、P3 = 200、P4 = 900。（将 P3 放在中间）
• C。更好的解决方案是 P1 = 0、P2 = 300、P3 = 600、P4 = 900。（均匀分配 P2 和 P3）。

这个问题的一些背景。看看下图。

• 我们需要为所有生命线（垂直线）创建一个布局。#1、#2
• 这些线之间的间隙由消息的长度决定。#3
• 我们需要确保图表的总宽度较小（以在参与者较多时节省空间）。#4
• 在此图中，B 应移至更靠近 C 的位置。#5

我已经实现了以下满足约束#1~#4 的内容。有人可以帮忙解决约束#5吗？它可以部分满足（如示例中的 b）或完全满足（如示例中的 c）。我们可以用它standard deviation来衡量解决方案的质量。 更新：感谢 Roman Svistunov，获得非最佳位置的基本实现被大大简化。

// Known minimum distance is stored in a matrix
// e.g.
//    m = [
//      [0, 100, 0, 900],
//      [0, 0,   0, 0],
//      [0, 0,   0, 0],
//      [0, 0,   0, 0],
//    ]

let f = function (n: number, m: Array<Array<number>>): number {
  let array = Array(n).fill(0);
  for (let i = 0; i < n; i++) {
    array[i] = f(i, m) + m[i][n];
  }
  return Math.max(0, ...array);
}

Answer 1

I\xe2\x80\x99ll 会研究一些理论，然后给出一些满足 1\xe2\x80\x934 和 6 的 JavaScript 代码，5 的修改版本，另外

\n

\n
7. 对数值问题具有鲁棒性
\n
8. 没有外部依赖。
\n

\n

我们可以将 1\xe2\x80\x934 表示为线性规划。例如，您的示例实例将翻译为

\n

minimize x4 - x1\nsubject to\n\n(known minimum distances)\na: x3 - x1 >= 200\nb: x4 - x1 >= 900\n\n(fixed order)\nc: x2 - x1 >= 0\nd: x3 - x2 >= 0\ne: x4 - x3 >= 0\n\nx1, x2, x3, x4 unbounded\n

\n

每个线性程序都有一个具有相同最佳目标值的对偶程序（在满足此处的条件下）。上述程序的对偶如下：

\n

maximize 200 a + 900 b\nsubject to\n\nx1: -a - b - c = -1\nx2: c - d = 0\nx3: a + d - e = 0\nx4: b + e = 1\n\na, b, c, d, e >= 0\n

\n

这是一个没有容量限制的最大成本流问题。因此，某些路径是最优流，因此我们可以将标准 1\xe2\x80\x934 理解为最长路径问题。

\n

一旦我们有了最佳目标（此处为 900），我们就可以将标准 4 转变为约束

\n

x4 - x1 <= 900\n

\n

并思考如何做出最令人赏心悦目的布局。正如所写，我们可以通过求解二次规划来实现标准 5：

\n

minimize (x2 - x1)**2 + (x3 - x2)**2 + (x4 - x3)**2\nsubject to\n\n(most efficient layout)\nx4 - x1 <= 900\n\n(known minimum distances)\nx3 - x1 >= 200\nx4 - x1 >= 900\n\n(fixed order)\nx2 - x1 >= 0\nx3 - x2 >= 0\nx4 - x3 >= 0\n\nx1, x2, x3, x4 unbounded\n

\n

我写的目标不是标准差，但这并不重要。最小化标准差相当于最小化方差。反过来，随机变量的方差X是E[X**2] - E[X]**2，并且E[X]**2在这里是一个常数（E[X] = 900/4 = 150在示例中）。最后，最大化E[X**2]相当于最大化同等可能结果的平方和。

\n

不过，我会提出一个不同的标准。

\n

5\xe2\x80\xb2。最大化（按字典顺序）从小到大排序的间隙长度列表。

\n

这意味着我们要尽可能避免出现小差距。在 Roman Svistunov 的一个例子中，假设我们有约束

\n

x3 - x1 >= 2\nx4 - x2 >= 5.5\nx5 - x3 >= 8\n

\n

具有最小标准差集的解决方案

\n

x1 = 0\nx2 = 0.75\nx3 = 2\nx4 = 6.25\nx5 = 10\n

\n

而优化替换目标集的解决方案

\n

x1 = 0\nx2 = 1\nx3 = 2\nx4 = 6.5\nx5 = 10\n

\n

我们\xe2\x80\x99已经变得x2更加居中，但代价是变得x4不那么居中。由于x2比 更接近它的邻居x4，这对我来说似乎是一个合理的交易，但你显然拥有最终决定权。

\n

为了优化 5\xe2\x80\xb2，我们进行如下操作。首先制定一个线性程序，例如

\n

maximize z\nsubject to\n\n(most efficient layout)\nx5 - x1 <= 10\n\n(known minimum distances)\nx3 - x1 >= 2\nx4 - x2 >= 5.5\nx5 - x3 >= 8\n\n(fixed order)\nx2 - x1 >= z\nx3 - x2 >= z\nx4 - x3 >= z\nx5 - x4 >= z\n\nx1, x2, x3, x4, x5 unbounded\n

\n

我们在这里求解最佳值，z = 1并记录哪些间隙阻止其变大。我们替换这些z变量并再次求解，重复直到不再出现z。

\n

maximize z\nsubject to\n\n(most efficient layout)\nx5 - x1 <= 10\n\n(known minimum distances)\nx3 - x1 >= 2\nx4 - x2 >= 5.5\nx5 - x3 >= 8\n\n(fixed order)\nx2 - x1 = 1\nx3 - x2 = 1\nx4 - x3 >= z\nx5 - x4 >= z\n\nx1, x2, x3, x4, x5 unbounded\n

\n

最优值是z = 3.5.

\n

maximize z\nsubject to\n\n(most efficient layout)\nx5 - x1 <= 10\n\n(known minimum distances)\nx3 - x1 >= 2\nx4 - x2 >= 5.5\nx5 - x3 >= 8\n\n(fixed order)\nx2 - x1 = 1\nx3 - x2 = 1\nx4 - x3 >= z\nx5 - x4 = 3.5\n\nx1, x2, x3, x4, x5 unbounded\n

\n

最优值是z = 4.5，\xe2\x80\x99 是最后一次迭代。

\n

正如前面所观察到的，我们可以使用扩展的最长路径算法来求解这些线性程序。在 JavaScript 中：

\n

// Dual numbers represent linear functions of time.\nfunction Dual(p, v) {\n  return { position: p, velocity: v };\n}\n\n// Adds dual numbers x and y.\nfunction dualPlus(x, y) {\n  return Dual(x.position + y.position, x.velocity + y.velocity);\n}\n\nconst epsilon = Math.sqrt(Number.EPSILON);\n\n// Compares dual numbers x and y by their position at a time infinitesimally\n// after now.\nfunction dualLessThan(x, y) {\n  let d = x.position - y.position;\n  return d < -epsilon || (Math.abs(d) <= epsilon && x.velocity < y.velocity);\n}\n\n// Tracks delta, the length of time for which every pair of arguments to\n// .dualLessThan() maintains the present relation.\nfunction DeltaTracker() {\n  return {\n    delta: Infinity,\n    dualLessThan: function (x, y) {\n      let lessThan = dualLessThan(x, y);\n      if (lessThan) {\n        [x, y] = [y, x];\n      }\n      if (x.velocity < y.velocity) {\n        this.delta = Math.min(\n          this.delta,\n          (x.position - y.position) / (y.velocity - x.velocity)\n        );\n      }\n      return lessThan;\n    },\n  };\n}\n\n// Converts the adjacency matrix to an adjacency list representation.\nfunction graphFromMatrix(n, matrix) {\n  let graph = [];\n  for (let j = 0; j < n; j++) {\n    graph.push([]);\n    for (let i = 0; i < j; i++) {\n      if (matrix[i][j] > 0) {\n        graph[j].push({ i: i, length: Dual(matrix[i][j], 0) });\n      }\n    }\n  }\n  return graph;\n}\n\n// Essentially the usual algorithm for longest path, but tracks delta.\nfunction longestPathTable(graph, gaps) {\n  let tracker = DeltaTracker();\n  let maximum = Dual(0, 0);\n  let table = [];\n  for (let j = 0; j < graph.length; j++) {\n    let argument = null;\n    if (j > 0) {\n      maximum = dualPlus(maximum, gaps[j - 1]);\n    }\n    for (let edge of graph[j]) {\n      let x = dualPlus(table[edge.i].maximum, edge.length);\n      if (tracker.dualLessThan(maximum, x)) {\n        argument = edge.i;\n        maximum = x;\n      }\n    }\n    table.push({ argument: argument, maximum: maximum });\n  }\n  return [tracker.delta, table];\n}\n\n// Essentially the usual algorithm for decoding the longest path from the\n// dynamic program table.\nfunction freezeCriticalGaps(table, graph, gaps) {\n  let j = table.length - 1;\n  while (j > 0) {\n    let critical = true;\n    let argument = table[j].argument;\n    if (argument !== null) {\n      j = argument;\n    } else {\n      j--;\n      gaps[j].velocity = 0;\n    }\n  }\n}\n\n// Changes the time from now to now + delta.\nfunction advanceGaps(gaps, delta) {\n  for (let i = 0; i < gaps.length; i++) {\n    gaps[i].position += gaps[i].velocity * delta;\n  }\n}\n\n// Extracts the final result from the dynamic program table.\nfunction positionsFromTable(table) {\n  let positions = [];\n  for (let entry of table) {\n    positions.push(entry.maximum.position);\n  }\n  return positions;\n}\n\n// Entry point for the layout algorithm.\nfunction layOut(n, matrix) {\n  let graph = graphFromMatrix(n, matrix);\n  let gaps = [];\n  for (let j = 1; j < n; j++) {\n    gaps.push(Dual(0, 1));\n  }\n  while (true) {\n    let [delta, table] = longestPathTable(graph, gaps);\n    if (delta == Infinity) {\n      return positionsFromTable(table);\n    }\n    if (table[n - 1].maximum.velocity > 0) {\n      freezeCriticalGaps(table, graph, gaps);\n    } else {\n      advanceGaps(gaps, delta);\n    }\n  }\n}\n\n// Various test cases.\nconsole.log(\n  layOut(4, [\n    [0, 200, 0, 900],\n    [0, 0, 0, 0],\n    [0, 0, 0, 0],\n    [0, 0, 0, 0],\n  ])\n);\n\nconsole.log(\n  layOut(5, [\n    [0, 0, 2, 0, 0],\n    [0, 0, 0, 5.5, 0],\n    [0, 0, 0, 0, 8],\n    [0, 0, 0, 0, 0],\n    [0, 0, 0, 0, 0],\n  ])\n);\n\nconsole.log(\n  layOut(4, [\n    [0, 0, 200, 0],\n    [0, 0, 0, 150],\n    [0, 0, 0, 0],\n    [0, 0, 0, 0],\n  ])\n);\n\nconsole.log(\n  layOut(4, [\n    [0, 0, 200, 0],\n    [0, 0, 0, 200],\n    [0, 0, 0, 0],\n    [0, 0, 0, 0],\n  ])\n);\n\nconsole.log(\n  layOut(4, [\n    [0, 0, 200, 0],\n    [0, 0, 0, 300],\n    [0, 0, 0, 0],\n    [0, 0, 0, 0],\n  ])\n);\n

\n

Answer 2

首先，我们应该找到该段的最小可能长度。这可以通过贪心算法来完成：将第一个点放置在 0（不失一般性），并用尽可能小的坐标将后面的点逐个相加，而不会违反到所有先前点的距离。我们会将最后一个点固定在其位置上，并且不会进一步移动它。我认为正确性是显而易见的。

现在棘手的部分是优化相邻距离的标准差。我们将引入 n-1 个变量——距离，我将它们命名为Di。尝试最小化标准差与最小化方差相同。

因为E(D)当我们固定第一个点和最后一个点时 是恒定的，所以最小化标准差与最小化差值之和相同。

对距离的约束都是线性的，因为无论何时m_d(Pi,Pj)你都会添加约束(sum Dk for k=i to j) >= m_d(Pi,Pj)。为了保证距离总和是我们需要的，我们还需要添加sum Dk小于或等于以及大于或等于第一步中计算的最小距离的合同。因此，这个问题是一个二次优化问题，可以用库来解决。

请注意，这些点的坐标可以很容易地从距离作为前缀和来恢复。

即使对于 30 个点，如果您想投入时间进行进一步优化（如果您打算增加点数），任何二次规划算法都可能足够快，因为这个问题实际上可以在多项式时间内精确解决。

第一个也是最重要的优化。由于目标函数是平方和，目标二次形式的矩阵是正定的（因为它是单位矩阵），因此存在多项式算法，如本文所示。

我想提到的最后一个优化可能不相关，但它对于大多数可能的 m_d 矩阵都非常好。为了减少二次规划问题的规模，我们可以在二次规划之前找到一些点的最佳位置。

首先，我们将减少矩阵，以确保对于任何一个矩阵，更新右侧i<j<k都是正确的。m_d(Pi,Pj)+m_d(Pj,Pk)>=m_d(Pi,Pk)请注意，此计算m_d(P1,Pn)也是该段的最小可能长度。此后，我们可以说，对于任何 i，如果m_d(P1,Pi)+m_d(Pi,Pn)=m_d(P1,Pn)为真，则点 Pi 只有一个可能的位置，并且，这允许我们删除其中一个变量，因为现在从 P(i-1) 到 Pi 的距离以及从 P(i-1) 到 Pi 的距离Pi 到 P(i+1) 可以用一个变量替换 - 从 P(i-1) 到 P(i+1) 的距离：第一个距离就是 Pi 位置 - P(i-1) 位置 可以导出通过归纳法，P(i+1)位置就是P(i-1)位置+这个新距离。

编辑2

这个伪代码有两个版本，都假设您使用某种二次规划库，第一个是直观的（并在一些库中使用），第二个是使用矩阵作为输入的库（但不太直观） ）：

function findOptimal(Matrix md, int n) {
  int minLength = findMinLength(md, n);
  Constraints constraints; //from library
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < m; ++j) {
      if (j < i)
        //Here sum of variables is not their current sum
        //but some object for library
        constraints.AddConstranit(sum of variables i to j >= md[j][i])
      else
        constraints.AddContraint(sum of variables j to i >= md[j][i])
    }
  }
  constraints.AddConstraint(sum of variables 0 to n-1 <= minLength)
  variables = QuadraticOptimizationSolve(sum of squares of variables, contracts)
  positions = Array<number>(n);
  positions[0] = 0;
  for (int i = 0; i < n; ++i) {
    positions[i + 1] = positions[i] + variables[i];
  }
  return positions;
}

由于矩阵在库中作为输入更为常见，因此以下是使用它们的伪代码：

function findOptimal(Matrix md, int n) {
  int minLength = findMinLength(md, n);
  Matrix optimization = Matrix(n - 1, n - 1);
  for (int i = 0; i < n - 1; ++i) {
    for (int j = 0; j < n - 1; ++j) {
      if (i == j)
        optimization[i][j] = 1
      else
        optimization[i][j] = 0
    }
  }
  Matrix constraints = Matrix(n * (n - 1) / 2 + 1, n)
  for (int i = 1; i < n; ++i) {
    for (int j = 0; j < i; ++j) {
      constraintRow = i * (i - 1) / 2 + j;
      //constriant is sum of distances from j to i is 
      //at least md[j][i]
      for (int k = j; k < i; ++k) {
        //these are coefficients in linear combination
        //negative because we need >= while optimizers usually use <= contraints
        if (k >= k and k < i)
          constraints[constraintRow][k] = -1
        else
          constraints[constraintRow][k] = 0
      }
      constraints[constraintsRow][n - 1] = -md[i][j];
    }
  }
  variables = QuadraticOptimizationSolve(optimization, constraints);
  positions = Array<number>(n);
  positions[0] = 0;
  for (int i = 0; i < n; ++i) {
    positions[i + 1] = positions[i] + variables[i];
  }
  return positions;
}
Also this solution assumes in your library the constants vector in constraints is the last column, but sometimes it might be separate vector

编辑3

这是针对您的问题的 javascript 解决方案。每秒轻松处理 30 个点（实际上，在我的设置中，它每秒最多可以处理 100 个点）。它使用与您不同的 findMinLength 函数实现（您的函数不保存递归调用的结果，因此是指数函数）。另外，我无法找到一个好的二次规划库，所以我最终使用了这个解决方案（在编辑 4 中进行了编辑）。这是针对 Node js 的，但合并了所有文件（并删除了 require 和导出，对于 vsmall，您需要在其文件中将 vsmall 定义为 epsilon），并且该库似乎可以工作。看起来上一个库只是这个库的浏览器改编版，而这个库有一个问题，在改编后几年就得到了解决，但手动完成它很容易（因为它没有依赖项）。

function findOptimal(Matrix md, int n) {
  int minLength = findMinLength(md, n);
  Constraints constraints; //from library
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < m; ++j) {
      if (j < i)
        //Here sum of variables is not their current sum
        //but some object for library
        constraints.AddConstranit(sum of variables i to j >= md[j][i])
      else
        constraints.AddContraint(sum of variables j to i >= md[j][i])
    }
  }
  constraints.AddConstraint(sum of variables 0 to n-1 <= minLength)
  variables = QuadraticOptimizationSolve(sum of squares of variables, contracts)
  positions = Array<number>(n);
  positions[0] = 0;
  for (int i = 0; i < n; ++i) {
    positions[i + 1] = positions[i] + variables[i];
  }
  return positions;
}

例如，调用findOptimal(4, m)where m 如您的示例中所示（有 4 个点），该函数返回[0,300,600,900]

编辑

正如 @david-eisenstat 中所指出的，第五点存在一种解释，即最大化最小间隙，但有一个更简单且渐近更快的解决方案，即答案二分搜索。

我们知道最小间隙至少为 0，最多为段长度 / n。我们可以通过二分搜索最大化它，并使用检查器确保每个间隙至少是我们的值。如果findMinLength是求解 1-4 的函数，则以下是优化该度量的解（伪代码）：

function findOptimal(Matrix md, int n) {
  minLength = findMinLength(md, n) //The answer to minimize segment length
  minGapL = 0 //Lower bound in maximized minimum gap length
  minGapR = ceiling(minLength / n) //Upper bound in maximized minimum gap
  while minGapL < minGapR {
    minGapTest = (minGapL + minGapR) / 2
    mdCopy = md.copy()
    //This loop ensures that every gap is at least minimum gap tested
    for (int i = 1; i < n; ++i) {
      mdCopy[i - 1][i] = max(mdCopy[i - 1][i], miGapTest)
    }
    int testLength = findMinLength(mdCopy, n)
    if (testLength == minLength) {
      //if it is possible to get tested minimum gap update lower bound
      minGapL = minGapTest
    } else {
      //otherwise update lower bound
      minGapR = minGapTest - 1;
    }
  }
  return (minLength, minGapL);
}

通过此函数返回长度和最大化 minGap，我们可以像循环中那样通过矩阵修改来重建答案，然后解决 1-4 的问题。该解决方案是多项式的（如上所示，1-4 问题可以在线性时间内解决，二分搜索是答案的对数，并且答案是输入大小的指数，它也是输入大小的线性）。