线性N元等式最小二乘问题

Question

假设我想找到两条任意高维线的“交点”。这两条线实际上不会相交，但我仍然希望找到最接近所有线的交点（即尽可能靠近所有线的点）。

假设这些线具有方向向量A、B和初始点C、D，我可以通过简单地设置一个线性最小二乘问题来找到最接近的交点：将线交点方程转换为最小二乘问题。

Ax + C = By + D

转换为最小二乘形式

[A, -B] @ [[x, y]] = D - C 

其中@代表矩阵与向量的乘法，然后我可以使用例如np.linalg.lstsq来解决它。

但是如何找到3条或更多任意直线的“最交点”？如果我按照同样的规则，现在有：

Ax + D = By + E = Cz + F

我能想到的唯一方法是将其分解为三个方程式：

Ax + D = By + E
Ax + D = Cz + F
By + E = Cz + F

将它们转换为最小二乘形式

[A, -B,  0]                 [E - D]
[A,  0, -C] @ [[x, y, z]] = [F - D]
[0,  B, -C]                 [F - E]

问题在于最小二乘问题的大小随着行数的增加呈二次增长。我想知道是否有更有效的方法来解决n元等最小二乘线性问题？

我在考虑上面提供的By + E = Cz + F的必要性，但由于这个问题没有精确的解（即它们实际上不相交），我认为这样做会对某些变量产生更多的“权重”？

谢谢你的帮助！

编辑

我刚刚使用以下代码测试了将第一个术语与所有其他术语配对的n元等式（而没有其他配对）

def lineIntersect(k, b):
    "k, b: N-by-D matrices describing N D-dimensional lines: k[i] * x + b[i]"

    # Convert the problem to least-square form `Ax = B`
    # A is temporarily defined 3-dimensional for convenience
    A = np.zeros((len(k)-1, k.shape[1], len(k)), k.dtype)
    A[:,:,0] = k[0]
    A[range(len(k)-1),:,range(1,len(k))] = -k[1:]

    # Convert to 2-dimensional matrix by flattening first two dimensions
    A = A.reshape(-1, len(k))

    # B should be 1-dimensional vector
    B = (b[1:] - b[0]).ravel()

    x = np.linalg.lstsq(A, B, None)[0]

    return (x[:,None] * k + b).mean(0)

下面的结果表明这样做是 不正确 的，因为n元等式中的第一项“权重不同”。

第一个输出是正常结果和不同输入顺序（行顺序无关紧要）的结果之间的差异，其中第一项没有改变。

第二个输出与第一项已更改相同。

k = np.random.rand(10, 100)
b = np.random.rand(10, 100)
print(np.linalg.norm(lineIntersect(k, b) - lineIntersect(np.r_[k[:1],k[:0:-1]], np.r_[b[:1],b[:0:-1]])))
print(np.linalg.norm(lineIntersect(k, b) - lineIntersect(k[::-1], b[::-1])))

导致

7.889616961715915e-16
0.10702479853076755

- ZisIsNotZis

---

那么，究竟是什么阻止你使用最小二乘法呢？为什么要进行所有这些无关的计算？ - n. m.

@n.m. 我相信这些东西并非无关紧要，因为它们是最小二乘法的输入。使用最小二乘法的问题在于它似乎具有Θ(n^2)的复杂度，其中n是线条数，这非常慢。我想知道是否有更有效的方法。 线性代数最小二乘需要两个输入，一个是上面最左边的矩阵，另一个是上面最右边的向量。 - ZisIsNotZis

经过编辑后现在更清晰了。复杂度应该是与维度数量的平方成正比，与行数成线性关系。 - n. m.

@n.m. 我认为它应该是与行数的平方成正比，因为在一个n元等式问题中，有n*(n-1)/2对相等的元素。如果我们还考虑LHS矩阵的行长度，那么矩阵实际上的形状是(n*(n-1)/2 * d, n)，这似乎是行数的立方和维度的线性关系？ - ZisIsNotZis

你需要最小化一个d变量函数（你正在寻找的点的坐标）。雅可比矩阵是一个d×d矩阵，矩阵的每个元素都有n个项。写下一些极端情况会有所帮助。1000000“行”，一维。2行，1000000维（1行在这里不起作用）。 - n. m.

Answer 1

另一个“几乎交点”的标准是一个点x，使得该点到这些直线的距离的平方和最小。与你的标准一样，如果这些直线实际上相交，则几乎交点将是实际交点。不过，我认为距离平方和的标准可以很容易地计算出所需的点：

假设我们用一个点和沿着该线的单位向量来表示一条线。因此，如果一条线由p，t表示，则该线上的点形式为

p + l*t for scalar l

点x到线p,t的距离平方为：

(x-p)'*(x-p) - square( t'*(x-p))

如果我们有N条线段p[i]，t[i]，那么从点x到这些线段的距离平方和为：

   Sum { (x-p[i])'*(x-p[i]) - square( t[i]'*(x[i]-p[i]))}

将其扩展，我得到了上面的结果

x'*S*x - 2*x'*V + K

在哪里

S = N*I - Sum{ t[i]*t[i]'}
V = Sum{ p[i] - (t[i]'*p[i])*t[i] }

而且K不依赖于x

除非所有的线都是平行的，否则S将是（严格）正定的，因此可逆，在这种情况下，我们的距离平方和为

(x-inv(S)*V)'*S*(x-inv(S)*V) + K - V'*inv(S)*V

因此，最小化的 x 是：

inv(S)*V

所以，步骤如下：将您的“方向向量”标准化（并使用与缩放方向相同的因子缩放每个点），按上述方法形成 S 和 V，解决。

S*x = V for x

- dmuir

---

太完美了！我在尝试将 x'*S*x - 2*x'*V + K 转换回 Ax=b 形式时卡住了，而你回答了这一部分！但我认为有一个打字错误：V 应该是 *t[i] } 而不是 *p[i] }。 - ZisIsNotZis

Answer 2

这个问题可能更适合数学堆栈交换。此外，有人有没有好的方法在这里格式化数学公式？很抱歉它很难阅读，我已经尽力使用Unicode。

编辑：我误解了@ZisIsNotZis在Ax+C行中的意思，所以请忽略下一段。

我不确定你的方法是否正确陈述。你介意发布你的代码和一个小的输出示例（也许是3或4条线的2D图形，这样我们就可以绘制它）吗？当你试图找到两条线的交点时，你不应该做Ax+C = Bx+D吗？如果你做Ax+C=By+D，你可以在第一条线上选一些x，在第二条线上选一些y，并完全满足两个方程。因为这里的x和y应该与空间的维度一样大，而不是标量。

寻找一点，使其与所有直线的距离尽可能接近的问题有许多理解方式。我认为最自然的方式是将每条直线到该点欧几里得距离的平方和最小化。

假设我们有一条线在R^n中: c^Tz + d = 0 (其中c是单位长度)，还有另一个点x。那么从x到该线的最短向量是：(I-cc^T)(x-d)，因此从x到该线的距离的平方是║(I-cc^T)(x-d)║^2。我们可以通过最小化这个距离来找到最靠近该线的点。请注意，这是一个标准的最小二乘问题，形式为min_x ║b-Ax║_2。

现在，假设我们有由 c_iz+d_i 给出的线，其中 i=1,...,m。点 x 到第 i 条线的平方距离为 d_i^2 = ║(I-cc^T)(x-d)║_2^2。我们现在想要解决 min_x \sum_{i=1}^{m} d_i^2 问题。

以矩阵形式表示为：

      ║ ⎡ (I-c_1 c_1^T)(x-d_1) ⎤ ║
      ║ | (I-c_2 c_2^T)(x-d_2) | ║
min_x ║ |      ...             | ║
      ║ ⎣ (I-c_n c_n^T)(x-d_n) ⎦ ║_2

这句话的意思是：这个系统的形式仍然是 min_x ║b - Ax║_2，因此有很好的求解器可用。每个块的大小为n（空间的维数），共有m个块（行数）。因此，该系统的大小为mn乘以n。特别地，它在行数上是线性的，在空间维数上是二次的。此外，它还有一个优点，即如果添加一条线路，只需将另一个块添加到最小二乘系统中即可。这也提供了在添加线路时迭代更新解决方案的可能性。我不确定是否有专门针对这种最小二乘系统的求解器。请注意，每个块都是单位矩阵减去秩一矩阵，因此可能会给出一些额外的结构，可以用来加速计算。尽管如此，我认为使用现有的求解器几乎总是比编写自己的求解器更好，除非您具有相当多的数值分析背景或需要解决非常专业的系统类别。

- tch

---

实际上，我试图在 min_x = mean(K[i] * x[i] + B[i]) 中找到 x，然后计算 min_x。我想直接计算 min_x 更有意义。 - ZisIsNotZis

关于 Ax+C = By+D，其中 A 是方向，C 是起点。因此，使用标量 x 的 Ax+C 简单地表示第一条线上的任意点。右侧遵循相同的规则。由于方向向量的大小和起点值没有限制（线上的任何点都可以），因此 x 和 y 不一定相等。 - ZisIsNotZis

我明白了，我误解你所说的行的含义。但是，一旦您为每条线找到标量，您不会以每条线一个点而不是单个点结束吗？ - tch

是的，然后我对它们取平均值，在两行的情况下完全有意义，但在 n 行的情况下就不确定了。 - ZisIsNotZis

对于再次在形式上的min_x ║b - Ax║_2部分，我认为它现在试图最小化一个矩阵范数，在该范数中，x的相同元素出现在不同的位置。由于最小二乘法不能像(ax+b)^2+(cx+d)^2=(ex+f)^2那样合并，我不确定该怎么做。 - ZisIsNotZis

Answer 3

不是解决方案，只是一些想法：

如果n维空间中的线具有参数方程（使用单位Dir向量）

 L(t) = Base + Dir * t

点 P 到此直线的平方距离为

W = P - Base
Dist^2 = (W - (W.dot.Dir) * Dir)^2

如果可以将Min(Sum(Dist[i]^2))以适合LSQ方法的形式编写（通过每个点坐标进行偏导数），则得到的系统可以解决(x1..xn)坐标向量。

（情况类似于通常LSQ的许多点和单条线的反转）

- MBo

Answer 4

你说你有两条“高维”线。这意味着表示线的矩阵具有比行更多的列。

如果是这种情况，并且您可以有效地找到一个低秩分解，使得A=LRᵀ，那么您可以将最小二乘问题的解重写为x=(Rᵀ RLᵀ L)⁻¹ Lᵀ y。

如果m是线的数量，n是线的维数，则将最小二乘时间复杂度从O(mn²+nʷ)减少到O(nr²+mr²)，其中r=min(m,n)。

问题在于如何找到这样的分解。

- Richard

---

这看起来很有趣。我需要一点时间来考虑一下。 - ZisIsNotZis

嗨，根据 @dmuir 的回答，现在我有一个矩阵 A，它是一个方阵（长度等于维数）。你知道任何更快的方法吗？ - ZisIsNotZis

@ZisIsNotZis：如果那个答案有一些合成数据可以玩耍就太好了。这个矩阵有什么特殊的结构吗？如果有，可能有加速的方法。如果没有，那么你最好找一个好的密集矩阵库，也许还要配备一张GPU。 - Richard