多变量线性回归的求解方法

多变量线性回归是一种统计方法，用于预测因变量 $y$ 与多个自变量 $x_1,x_2,\dots ,x_k$ 之间关系。目标是找到一组参数 $w_1,w_2,\dots ,w_k$和 $b$，使得模型 $z=x_1{w}_1+x_2{w}_2+\cdots +x_k{w}_k+b$ 能够最好地拟合数据。例如，在房价预测中，影响因素可能包括居住面积、地理位置、朝向、学区房、周边设施、建筑年份等。

在建立多元线性回归模型时，为了确保模型具有良好的解释能力和预测效果，自变量的选择应遵循以下准则：

• 自变量必须对因变量有显著影响，并呈现紧密的线性相关性；
• 自变量与因变量之间的相关性应为真实的，而非形式上的；
• 自变量之间应具有适度的互斥性，其相关性不应超过自变量与因变量之间的相关性；
• 自变量应具备完整的统计数据，以便预测值容易确定；

接下来，将讨论两种求解多变量线性回归参数的方法：正规方程解法和梯度下降解法。

---

正规方程解法

正规方程解法是通过解一个线性方程组来找到参数的精确解。这种方法基于最小二乘法，即最小化预测值和实际值之间的平方误差之和。

给定数据集 $x_1,x_2,\dots ,x_k,y$，假设我们有 $n$ 个样本，我们可以将数据写成矩阵形式：

$$X=\left[\begin{array}{cccc}{x}_{11}& x_{12}& \cdots & x_{1k}\\ x_{21}& x_{22}& \cdots & x_{2k}\\ ⋮& ⋮& \ddots & ⋮\\ x_{n1}& x_{n2}& \cdots & x_{nk}\end{array}\right],\mathbf{y}=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_n\end{array}\right]$$

其中 $X$ 是设计矩阵，$\mathbf{y}$ 是目标向量。

正规方程为：

$$(X^{T}X)\mathbf{w}=X^T\mathbf{y}$$

其中 $\mathbf{w} $ 是参数向量 $[w_1,w_2,\dots ,w_k{]}^T$，$b$ 通常通过在模型中添加截距项来处理。

解正规方程得到：

$$\mathbf{w}=(X^{T}X{)}^{-1}{X}^T\mathbf{y}$$

Python代码实现：

import numpy as np

class LinearRegression:
    def __init__(self):
        self.w = None
        self.b = None

    def fit(self, X, y):
        # 添加截距项
        X = np.concatenate([np.ones((X.shape[0], 1)), X], axis=1)
        # 计算参数
        self.w = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
        self.b = self.w[0]

    def predict(self, X):
        # 添加截距项
        X = np.concatenate([np.ones((X.shape[0], 1)), X], axis=1)
        return X.dot(self.w)


# 示例数据
X = np.array([[1, 2], [2, 3], [3, 4], [4, 5]])
y = np.array([1, 2.5, 3.5, 5])
model = LinearRegression()
model.fit(X, y)
print("Weights:", model.w)
print("Bias:", model.b)

---

梯度下降解法

梯度下降解法是一种迭代优化算法，用于找到最小化成本函数的参数。在多变量线性回归中，成本函数通常是平方误差之和：

$$J(\mathbf{w},b)=\frac{1}{2n}\sum _{i=1}^n(y_i-(x_{i1}{w}_1+x_{i2}{w}_2+\cdots +x_{ik}{w}_k+b){)}^2$$

梯度下降通过迭代更新参数来最小化这个成本函数：

$$\mathbf{w}:=\mathbf{w}-\alpha \frac{\partial J}{\partial \mathbf{w}},b:=b-\alpha \frac{\partial J}{\partial b}$$

其中 $\alpha$ 是学习率，$\frac{\partial J}{\partial \mathbf{w}}$ 和 $\frac{\partial J}{\partial b}$ 分别是成本函数相对于 $\mathbf{w}$ 和 $b$ 的偏导数。

Python代码实现：

class LinearRegressionSGD:
    def __init__(self, learning_rate=0.01, max_iter=1000):
        self.learning_rate = learning_rate
        self.max_iter = max_iter
        self.w = None
        self.b = None

    def fit(self, X, y):
        # 添加截距项
        X = np.concatenate([np.ones((X.shape[0], 1)), X], axis=1)
        # 初始化参数
        self.w = np.zeros(X.shape[1])
        
        for _ in range(self.max_iter):
            gradients = 2/X.shape[0] * X.T.dot(X.dot(self.w) - y)
            self.w -= self.learning_rate * gradients
            if _ % 100 == 0:
                loss = np.mean((y - X.dot(self.w)) ** 2)
                print(f"Iteration {_}, Loss: {loss}")

    def predict(self, X):
        # 添加截距项
        X = np.concatenate([np.ones((X.shape[0], 1)), X], axis=1)
        return X.dot(self.w)


# 示例数据
X = np.array([[1, 2], [2, 3], [3, 4], [4, 5]])
y = np.array([1, 2.5, 3.5, 5])
model = LinearRegressionSGD(learning_rate=0.01, max_iter=1000)
model.fit(X, y)
print("Weights:", model.w)
print("Bias:", model.w[0])

---

总结

• 正规方程解法：提供了参数的精确解，计算复杂度为 $ O(k^3) $，适用于较小的数据集。
• 梯度下降解法：是一种迭代方法，可以处理更大的数据集，但需要选择合适的学习率和迭代次数。

在实际应用中，选择哪种方法取决于数据集的大小、特征的数量以及计算资源。对于大规模数据集，通常更倾向于使用梯度下降或其变体（如随机梯度下降）。

方法	正规方程	梯度下降
原理	几次矩阵运算	多次迭代
特殊要求	$X^{\top }X$ 的逆矩阵存在	需要确定学习率
复杂度	$O(n^3)$	$O(n^2)$
适用样本数	$m<10000$	$m\ge 10000$