algorithm - 3D 线性回归

Question

我想编写一个程序，给定 3D 空间中的点列表，表示为浮点中的 x、y、z 坐标数组，在此空间中输出最佳拟合线。线可以/应该是单位向量和线上的一个点的形式。

问题是我不知道如何做到这一点。我发现的最接近的是 this链接，虽然老实说我不明白他是如何从一个方程到另一个方程的，当我们到达矩阵时，我已经很迷茫了。

是否有我可以使用的简单 2D 线性回归的概括/有人可以(数学上)解释上述链接方法是否/如何工作(以及使用它计算最佳拟合线必须做什么)？

Answer 1

线性回归

N维线性回归的标准公式为



结果在哪里，是一个大小为 n + 1 的向量，给出了最适合数据的函数系数。

在您的情况下 n = 3。而 X 是称为设计矩阵的 mx(n+1) 矩阵——在您的情况下为 mx4。要构 build 计矩阵，您只需将每个数据点坐标值 (x1,x2,...) 复制到 X 的一行中，此外，将数字 1 放在每行的第 1 列中。矢量 y 具有与这些坐标相关联的值。条款 和 是“X 的转置”和“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。
在这种情况下，



然后你只需要把东西相乘，你就可以得到 直接地。矩阵相乘很简单，虽然更复杂，但求矩阵的逆相当简单( see here for example )。然而，对于像 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-X 转置时，这是必要的(不是这次)，如果您选择自己实现矩阵求逆，请记住这一点。

取而代之的是 PCA

主成分分析(PCA)是一种降维技术，目标是从一个n维空间中找到一个k维空间，使得投影误差最小。在一般情况下，n 和 k 是任意的，但在这种情况下，n = 3 和 k = 1。有 4 个主要步骤。

第 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:执行奇异值分解

在这里最好只找到一个库，它采用协方差矩阵并吐出答案。有很多，这里有一些； in python in R , in Java ，当然，在 Octave、Julia、Matlab(如 R)中，它是另一个衬垫 svd .

对协方差矩阵执行 SVD

(U,S,V) = svd((1/4)*Xm'*Xm);

第 4 步:找到线

取第一个组件(对于 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