我想编写一个程序，给定一个由浮点型x、y、z坐标数组表示的3D空间中的点列表，输出该空间中最佳拟合线。这条线可以/应该以单位向量和线上的一点形式表示。
问题是我不知道如何实现。我找到的最接近的方法是这个链接，但说实话我并不理解他是如何从方程推导出其他方程的，而当我们讲到矩阵时，我已经有些迷失了。
是否有简单的2D线性回归的概括方法可用于此，或者是否有人能够解释（数学上）上述链接中的方法如何工作（以及计算使用它来计算最佳拟合线需要做什么）？
3D线性回归
16
- Jimmy
3
这个类似的问题/答案可能会有所帮助 https://dev59.com/LnrZa4cB1Zd3GeqP46OU#20888844 只有2D示例，但算法仍适用于您的需求。使用方向检查方法（将3D编码的单位切向量转换为标量）。 - Spektre
那种方法不一定会给出最佳拟合线，是吗？ - Jimmy
是的，但它会过滤掉错误的点，并为您提供一个不错的起点，以加快线拟合的速度。 - Spektre
5个回答
42
线性回归
N维线性回归有一个标准公式，如下所示：
其中结果
是大小为n + 1的向量，给出最适合数据的函数的系数。
在你的情况下n = 3。X是一个mx(n + 1)的矩阵，称为设计矩阵，在你的情况下是mx4。构建设计矩阵时，只需将每个数据点坐标值(x1,x2,...)复制到X的一行中，并在每行的第1列上放置数字1。向量y具有与这些坐标相关联的值。项
和
是“X的转置”和“X的转置和积的逆”。计算后一项可能需要大量的计算量，因为求矩阵的逆是O(n^3)，但对于n = 4的情况而言，只要n小于5000就没有问题。
一个例子
假设你有数据点(6,4,11) = 20, (8,5,15) = 30, (12,9,25) = 50和(2,1,3) = 7。 在这种情况下，
然后你只需要将它们相乘，就可以直接得到。矩阵乘法是简单的，虽然更复杂，但求一个矩阵的逆是相当简单的(可参考此处)。 然而，对于Matlab、Octave和Julia等科学计算语言，只需一行代码即可完成。
julia> X = [1 6 4 11; 1 8 5 15; 1 12 9 25; 1 2 1 3]
4x4 Array{Int64,2}:
1 6 4 11
1 8 5 15
1 12 9 25
1 2 1 3
julia> y = [20;30;50;7]
4-element Array{Int64,1}:
20
30
50
7
julia> T = pinv(X'*X)*X'*y
4-element Array{Float64,1}:
4.0
-5.5
-7.0
7.0
正在验证...
julia> 12*(-5.5) + 9*(-7.0) + 25*(7) + 4
50.0
在Julia、Matlab和Octave中，可以使用*来简单地对矩阵进行乘法运算，而转置操作符是'。请注意，在这里我使用了pinv(伪逆)，当数据过于冗余并导致一个非可逆的X-Xtranspose时，这是必需的(不是这次)，如果你选择自己实现矩阵求逆，请记住这一点。
使用PCA代替
主成分分析(PCA)是一种降维技术，其目标是从n维空间中找到一个k维空间，使投影误差最小化。在一般情况下，n和k是任意的，但在这个例子中，n=3，k=1。有四个主要步骤。第1步：数据预处理
为了使标准方法有效，首先必须执行均值归一化，并可能对数据进行缩放，以避免由于浮点误差导致算法失败。在后一种情况下，这意味着如果某个维度的值范围与另一个维度相比巨大，可能会出现问题(例如，在一个维度中-1000到1000与-0.1到0.2相比)。通常它们足够接近。均值归一化只是表示对于每个维度，从每个数据点中减去平均值，以便得到的数据集围绕原点居中。将结果取出，并将每个数据点(x1、x2、...xn)作为一个大矩阵X中的一行储存。X = [ 6 4 11; 8 5 15; 12 9 25; 2 1 3]
4x3 Array{Int64,2}:
6 4 11
8 5 15
12 9 25
2 1 3
找出平均值
y = convert(Array{Float64,1},([sum(X[1:4,x]) for x = 1:3])/4')
3-element Array{Float64,1}:
7.0
4.75
13.5
标准化...
julia> Xm = X .- y'
4x3 Array{Float64,2}:
-1.0 -0.75 -2.5
1.0 0.25 1.5
5.0 4.25 11.5
-5.0 -3.75 -10.5
步骤2：计算协方差矩阵
协方差矩阵sigma可以简单地表示为：
其中m是数据点的数量。
步骤3：执行奇异值分解
在这一步，最好找到一个库，将协方差矩阵作为输入并输出答案。有很多这样的库，以下是其中几个：Python、R、Java，当然还有Octave、Julia以及像R和Matlab一样的另外一个一行代码的函数svd。
对协方差矩阵执行SVD。
(U,S,V) = svd((1/4)*Xm'*Xm);
步骤四：找到这一行
取第一个组件（对于k维数据，您需要取前k个组件）。
Ureduce = U[:,1]
3-element Array{Float64,1}:
-0.393041
-0.311878
-0.865015
这是将投影误差最小化的线。
额外加分：回溯
你甚至可以恢复原始值的近似，但它们将全部排列在同一条线上并进行投影。连接这些点以获得一条线段。
获取X中每个数据点的降维（因为1-D将只有1个值）：
z= Ureduce' * Xm'
1x4 Array{Float64,2}:
2.78949 -1.76853 -13.2384 12.2174
回到另一条路上；原始值，但全部位于同一（最佳）线上。
julia> (Ureduce .* z .+ y)'
4x3 Array{Float64,2}:
5.90362 3.88002 11.0871 6 4 11
7.69511 5.30157 15.0298 versus 8 5 15
12.2032 8.87875 24.9514 12 9 25
2.19806 0.939664 2.93176 2 1 3
- waTeim
11
@waTeim 我突然有一种想给你买啤酒的冲动。 - easymoden00b
2@easymoden00b 哇，我就知道这个stackoverflow的东西迟早会有回报！ - waTeim
2@waTeim 尽管我对第4步中的三个值如何代表三维空间中一条线的值有些困惑。但我该如何使用这三个值来建立一个 (x,y) -> (z) 函数呢？ - easymoden00b
3有人赐予这个人一块饼干，因为他/她值得获得！ - Dennis Braga
@DennisBraga 你怎么从最后一部分获取一行？就像 easymoden00b 说的那样？ - hownowbrowncow
显示剩余6条评论
7
@WaTeim提供了很好的答案。
这是我的Python代码，适用于需要它的人。与提供的数值示例一起使用。
def regr(X):
y= np.average(X, axis=0)
Xm = X-y
u, s, v = np.linalg.svd((1./X.shape[0])*np.matmul(Xm.T,Xm))
# Extra Credit: Going back
z= np.matmul(u[:,0].T, Xm.T)
c = np.array([z*n for n in u[:,0]])
d = np.array(y.tolist()*c.shape[1]).reshape(c.shape[1],-1).T
e = (c+d).T
return u,s,v
regr(np.array([[6, 4, 11],[8,5,15],[12,9,25],[2,1,3]]))
顺便问一下，有人能告诉我为什么numpy的np.cov()与1./X.shape[0])*np.matmul(Xm.T,Xm)会给出不同的结果吗？
- Chris
2
在三维空间中找到适合给定点列表的最佳拟合直线是一项相当困难的任务。可以使用两个向量定义三维空间中的一条直线：点a，它位于直线上，以及线（归一化）方向n。可以使用以下方程描述它，其中t是实数。
以下是使用Armadillo库（C ++）对点集{(1,1,1), (2,2,2), (3,3,3)}进行演示的示例：
以下是使用Armadillo库（C ++）对点集{(1,1,1), (2,2,2), (3,3,3)}进行演示的示例：
#include<armadillo>
#include<vector>
using namespace arma;
int main()
{
std::vector<vec3> points {{
{1, 1, 1},
{2, 2, 2},
{3, 3, 3}
}};
int N = points.size();
vec3 mean = {0, 0, 0};
mat33 corr(fill::zeros);
for(auto p : points)
{
mean += p;
for(int i = 0; i < 3; i++)
for(int j = i; j < 3; j++)
corr(i, j) += p(i) * p(j);
}
corr /= N;
mean /= N;
mat33 cov {{
{corr(0, 0) - mean(0) * mean(0), corr(0, 1) - mean(0) * mean(1), corr(0, 2) - mean(0) * mean(2)},
{corr(0, 1) - mean(0) * mean(1), corr(1, 1) - mean(1) * mean(1), corr(1, 2) - mean(1) * mean(2)},
{corr(0, 2) - mean(0) * mean(2), corr(1, 2) - mean(2) * mean(1), corr(2, 2) - mean(2) * mean(2)}
}};
vec3 eigval; mat33 eigvec;
eig_sym(eigval, eigvec, cov);
mean.print("\nPoint: ");
eigvec.col(eigval.index_max())
.print("\nDirection:");
}
- andywiecko
1
这可以通过使用NPM的ml-matrix一行代码来实现：
const { Matrix, solve } = require('ml-matrix');
solve(this.DataX, Matrix.columnVector(this.DataY[0]));
- Jason Mullings
2
这个问题与编程语言无关，所以至少你应该说一下你指的是哪种语言... - m93a
1我的猜测会是JavaScript，因为有提到npm，如果这有帮助的话 @m93a - CPSuperstore
-1
我在 QT 代码中使用这个简单的方法：
QPair<QVector3D, QVector3D> getLineByLeastSquares(const QVector<QVector3D>& points)
{
if (points.size() <= 1)
return QPair<QVector3D, QVector3D>();
QVector3D avg;
for (const QVector3D& p : points)
avg += p;
avg /= static_cast<float>(points.size());
float nX = 0.0F, nY = 0.0F, nZ = 0.0F;
for (const QVector3D& p : points)
{
const QVector3D tmp = p - avg;
nX += tmp.x() * tmp.x();
nY += tmp.x() * tmp.y();
nZ += tmp.x() * tmp.z();
}
return QPair<QVector3D, QVector3D>(avg, QVector3D(nX, nY, nZ).normalized());
}
结果的第一个组件是QPair<QVector3D，QVector3D>，它是线的点，第二个组件是线的法向量。
- Michal Steller
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接