7

我有n实数变量（不知道，不在乎），我们称它们为X[n]. 我也有m >> n他们之间的关系，我们称他们R[m]为 ，形式为：

X[i] = alpha*X[j],alpha是一个非零正实数，i并且j是不同的，但该(i, j)对不一定是唯一的（即具有不同 alpha 因子的相同变量之间可能存在两种关系）

我正在尝试做的是找到一组alpha参数，以某种最小二乘的方式解决超定系统。理想的解决方案是最小化每个方程参数与其选择值之间的差异平方和，但我对以下近似值感到满意：

如果我将 m 个方程转换为 n 个未知数的超定系统，任何基于伪逆的数值求解器都会给我一个明显的解（全为零）。所以我目前所做的是在混合中添加另一个方程x[0] = 1（实际上任何常数都可以）并使用 Moore-Penrose 伪逆以最小二乘法求解生成的系统。虽然这试图最小化 的总和(x[0] - 1)^2和 的平方和x[i] - alpha*x[j]，但我发现它是我的问题的一个很好且数值稳定的近似值。这是一个例子：

a = 1
a = 2*b
b = 3*c
a = 5*c

在八度：

A = [
  1  0  0;
  1 -2  0;
  0  1 -3;
  1  0 -5;
]

B = [1; 0; 0; 0]

C = pinv(A) * B or better yet:
C = pinv(A)(:,1)

a这会产生, b, c:的值，[0.99383; 0.51235; 0.19136] 这给了我以下（合理的）关系：

a = 1.9398*b
b = 2.6774*c
a = 5.1935*c

所以现在我需要在 C/C++/Java 中实现这个，我有以下问题：

有没有更快的方法来解决我的问题，或者我在生成超定系统和计算伪逆方面是否走在正确的轨道上？

m我目前的解决方案需要一个奇异值分解和三个矩阵乘法，考虑到可能是 5000 甚至 10000有点多。是否有更快的方法来计算伪逆（实际上，我只需要它的第一列，而不是给定矩阵的稀疏性（每行恰好包含两个非零值，其中一个始终为一个，另一个始终为负）

您建议为此使用哪些数学库？LAPACK好吗？

我也对任何其他建议持开放态度，只要它们在数值上稳定且渐近快速（比如说k*n^2，哪里k可以很大）。

4

2 回答 2

4

你的问题不合适。如果您将问题视为 n 个变量的函数（差的最小二乘），则该函数恰好具有一个全局最小超平面。

该全局最小值将始终包含零，除非您将其中一个变量固定为非零，或以其他方式减少函数域。

如果您想要的是解决方案超平面的参数化，您可以从 Moore-Penrose Pseudo-Inverse http://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverse获得，并查看获取所有解决方案。

（请注意，我以技术上不正确的方式使用了“超平面”一词。我的意思是参数空间的“平面”无界子集......一条线，一个平面，可以由一个或多个向量参数化的东西。出于某种原因，我找不到此类对象的通用名词）

于 2011-08-19T09:32:36.813 回答
3

SVD 方法在数值上非常稳定，但速度不是很快。如果您使用 SVD，那么 LAPACK 是一个很好的库。如果它只是一次性计算，那么它可能已经足够快了。

如果你需要一个更快的算法，你可能不得不牺牲稳定性。一种可能性是使用 QR 分解。您必须阅读此内容才能查看详细信息，但部分原因如下。如果 AP = QR（其中 P 是置换矩阵，Q 是正交矩阵，R 是三角矩阵）是 A 的经济 QR 分解，则方程 AX = B 变为 QRP^{-1} X = B并且解决方案是 X = PR^{-1} Q^T B。以下 Octave 代码使用与您的代码中相同的 A 和 B 来说明这一点。

[Q,R,P] = qr(A,0)
C(P) = R \ (Q' * B)

这样做的好处是您可以通过进行稀疏 QR 分解来利用 A 的稀疏性。Octave 帮助中有一些关于 qr 功能的解释，但它并没有立即对我有用。

更快（但更不稳定）是使用正规方程：如果 AX = B 则 A^TAX = A^T B。矩阵 A^TA 是（希望）满秩的方阵，所以你可以使用任何线性方程的求解器。八度代码：

C = (A' * A) \ (A' * B)

同样，可以在这种方法中利用稀疏性。求解稀疏线性系统的方法和库有很多；一个流行的似乎是UMFPACK。

后来补充：我对这个领域了解的不够多，无法量化。整本书都写在这上面。也许 QR 的速度大约是 SVD 的 3 或 5 倍，而正规方程的速度又快了两倍。对数值稳定性的影响取决于您的矩阵 A。稀疏算法可以快得多（比如 m 的一个因子），但它们的计算成本和数值稳定性在很大程度上取决于问题，其方式有时并不十分清楚。

在您的用例中，我的建议是尝试使用 SVD 计算解决方案，看看需要多长时间，如果可以接受，那么就使用它（我想 n=1000 和 m=10000 大约需要一分钟）。如果你想进一步研究它，也可以试试 QR 和正规方程，看看它们的速度有多快，有多准确；如果它们提供与 SVD 大致相同的解决方案，那么您可以非常确信它们对于您的目的来说足够准确。仅当这些都太慢并且您愿意花一些时间在其中时，请查看稀疏算法。

于 2011-08-20T02:56:21.377 回答