Rad*_*Dan 7 matlab algebra numerical-methods least-squares
我有n实数变量(不知道,不关心),让我们称之为X[n].我也有m >> n他们之间的关系让我们称他们R[m]为:
X[i] = alpha*X[j],alpha是一个非零正实数,i并且j是不同的,但该(i, j)对不一定是唯一的(即,具有不同α因子的相同变量之间可以存在两个关系)
我要做的是找到一组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
在Octave:
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.有更快的方法来计算伪逆(实际上,我只需要它的第一列,而不是鉴于矩阵的稀疏性(每行包含两个非零值,其中一个总是一个而另一个总是负的),给定B除零第一行外的整个矩阵)
您建议使用哪些数学库？LAPACK好吗？
我也对任何其他建议持开放态度,条件是它们在数值上稳定并且渐近快速(假设k*n^2它们k可以很大).
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 大致相同，那么您可以非常确信它们对于您的目的来说足够准确。只有当这些都太慢并且您愿意投入一些时间时，才考虑稀疏算法。