Projection Matrix and Linear Regression

本文详细介绍了向量投影的数学原理,并将其应用于线性回归的最小二乘法求解过程中,解释了为何选择最小二乘作为评估标准。

假如想把一个向量投影到一个空间,我们想得到投影后的向量,那么我只需要用一个投影矩阵乘以原来的向量就能得到投影后的矩阵,那么如何得到这个投影矩阵呢,下面开始推导,首先从一维空间开始:

对于上面这幅图来说,a,b都是向量,x是一个数,可以得到:


解出:

           

因为aTb是一个数字,所以我们能将数字和向量a交换,得到投影后的向量xa。对于投影到一维空间来说,

就是投影矩阵。

投影到多维空间时,有:


x1*a1+x2*a2+...+xn*an则是投影后的向量,写成矩阵形式就是:


可以解出:


其中


就是投影矩阵,A=(a1,a2,...,an),a1,a2,an为该空间的基。X=(x1,x2,..xn)T,x1,x2..xn分别为基的线性组合的系数,可以发现多维的投影矩阵形式和一维有点相似,一维也可以写成多维的这种形式。

投影有个非常重要的性质就是:假如我想在A的列空间找一个离b最近的向量(最近指的是,向量每一维差值平方和最小),那么最近的向量就是b投影到A空间的投影向量。

我之所以说重要,因为他关系到线性回归最小二乘法的解,假如现在我们要做直线(多维情况下应该说超平面)拟合,方程为:AX=y,一般情况下,这个方程是无解的,不然也不需要拟合了,拟合效果好坏的评价就是最小二乘,就是求解Ax=y',其中,y'在A的列空间里(即AX=y'有解),并且要保证,y与y'的距离最小。根据前面投影的性质,这个y’就是y在A的列空间的投影。那么我们我们求解Ax=y',解出的x向量就是最佳直线(超平面)。有:

            

解出这个有解的方程,就得到最佳拟合直线(超平面)了。

对于最小二乘的来源,我觉得还是有必要说明下,为什么我们要用最小二乘来衡量优劣,而不用四次方,八次方,或者绝对值的和呢,据我所知,这是根据假设yi满足均值为y'i,方差为σ的高斯分布所得出来的结果,即y服从

分布(其中,y=(y1,y2,...ym)T,y'=(y'1,y'2,...y'm)T,i=1..m,m为点的个数,然后可以写出似然函数:



发现中括号里面的前面部分不依赖y'i,不影响求最大值,可以省去,最后我们其实求


这个式子的最小值,这与最小二乘是同一个目的。

最后我想说的是,有时候虽然叫线性回归的拟合,其实可以是多项式拟合或者其他拟合,比如现在有m个点,每个点数据都是n维的,那么最简单的拟合方法就是,将每一维直接映射过来(注意可以多加一维,每个点都设为1,作为常数项),但是我也可以把第二维的平方,第三维的立方....映射过来,这样的结果就是多项式拟合了。

另外,这样线性回归的方法会导致过拟合的问题,关于过拟合的意思,prml里面有详细的讲解,我的理解大概就是,我为了尽量拟合给定数据集,而丧失了数据本身的规律,而利用贝叶斯思想,能一定程度上克服这个问题。


import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as sns sns.set(context="notebook", style="whitegrid", palette="deep") # 读入数据集 data = pd.read_csv('data/multiple_linear_regression_dataset.csv') fig = plt.figure(figsize=(10, 8)) ax = fig.add_subplot(111, projection='3d') # 提取数据 age = data['age'] experience = data['experience'] income = data['income'] # 绘制三维散点图 ax.scatter(age, experience, income, c='blue', marker='o', label='Data Points') # 添加标签和标题 ax.set_xlabel('Age') ax.set_ylabel('Experience') ax.set_zlabel('Income') ax.set_title('Three-Dimensional Scatter Plot') # 显示图形 plt.show() # 提取 x 向量和 y 向量 X = data[['age', 'experience']] X['bias'] = 1 # 添加额外的 1 列 Y = data['income'] # 判断X.T * X是否是满秩矩阵 print(np.linalg.matrix_rank(X.T @ X) == min(X.shape[0], X.shape[1])) import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from mpl_toolkits.mplot3d.art3d import Poly3DCollection import numpy as np # 定义直线方程 def linear_equation(age, experience): return k[0] * age + k[1] * experience + k[2] # 选择一部分数据进行绘制 sample_size = min(50, len(age)) # 选择样本数量为50或总体数量的最小值 sample_indices = np.random.choice(len(age), size=sample_size, replace=False) # 创建三维图形 fig = plt.figure(figsize=(10, 8)) ax = fig.add_subplot(111, projection='3d') # 绘制散点图,调整透明度和颜色映射 scatter = ax.scatter(age[sample_indices], experience[sample_indices], income[sample_indices], c=income[sample_indices], cmap='viridis', marker='o', s=80, alpha=0.6, label='Data Points', zorder=2) # 生成一些数据点用于绘制,调整曲面高度 age_points = np.linspace(20, 60, 200) experience_points = np.linspace(0, 20, 200) age_grid, experience_grid = np.meshgrid(age_points, experience_points) income_pred = linear_equation(age_grid, experience_grid) # 绘制曲面,调整透明度 ax.plot_trisurf(age_grid.flatten(), experience_grid.flatten(), income_pred.flatten(), cmap='viridis', alpha=0.2, linewidth=0, zorder=1) # 设置轴标签 ax.set_xlabel('Age') ax.set_ylabel('Experience') ax.set_zlabel('Income') # 调整视角 ax.view_init(elev=40, azim=30) # 显示图例 ax.legend(handles=[scatter], labels=['Data Points'], loc='upper left') # 显示图形 plt.show() 将上面这段代码改为在pytorch环境下的代码
最新发布
03-19
<think>嗯,用户给了一段Python代码,主要是用numpy、pandas、matplotlib和seaborn进行数据处理和可视化,特别是3D散点图和曲面的绘制。然后用户希望将这段代码转换为PyTorch环境下的代码。首先,我需要理解原代码的功能,然后找出需要替换的部分。 首先,原代码导入了必要的库,读取数据,绘制三维散点图,添加了偏置列,检查矩阵的秩,然后定义了一个线性方程,并绘制了拟合的平面。用户想用PyTorch来实现,可能需要用PyTorch的张量操作代替numpy的部分,但其他部分如pandas和matplotlib可能不需要改动,因为PyTorch主要处理的是张量计算,而绘图还是依赖matplotlib。 首先,原代码中的X和Y是pandas的DataFrame和Series,这里需要转换为PyTorch张量。比如,X = data[['age', 'experience']],添加bias列后,转换为张量。Y同理。然后判断X.T @ X的秩,这里PyTorch有torch.matrix_rank函数可以使用。不过原代码中的判断是否满秩的逻辑可能需要调整,因为PyTorch的张量操作和numpy略有不同。 接下来,定义线性方程的部分,原代码中的k应该是系数,可能在PyTorch中需要用模型参数来表示,比如使用torch.nn.Linear层,或者手动定义参数。但原代码中的k可能是通过正规方程计算得到的,所以需要检查原代码是否使用了正规方程。原代码中没有显示计算k的部分,可能在用户提供的代码中被省略了?或者可能在后面的部分,比如用np.linalg求解的? 用户提供的代码中有判断X.T @ X是否是满秩矩阵,但没有显示计算系数的步骤。可能用户的实际代码中,k是通过正规方程得到的,比如k = np.linalg.inv(X.T @ X) @ X.T @ Y。那么在转换为PyTorch时,需要将矩阵运算转换为PyTorch的张量运算,比如使用torch.inverse()函数。 所以,步骤大概是: 1. 将X和Y转换为PyTorch张量。 2. 计算X的转置与X的乘积,判断其秩。 3. 使用正规方程计算系数k。 4. 绘制结果时,可能需要将张量转换回numpy数组,因为matplotlib需要numpy的数据。 需要注意的点: - PyTorch张量需要设置dtype,比如torch.float32。 - 矩阵相乘在PyTorch中使用@运算符或者torch.matmul()。 - 计算逆矩阵时,使用torch.inverse(),但需要确保矩阵是可逆的。 - 当需要将张量转换为numpy数组时,使用.detach().numpy(),尤其是当张量有梯度时。 另外,绘图部分大部分不需要改动,只需要确保输入到matplotlib的数据是numpy数组。因此,在计算完k之后,需要将其转换为numpy数组,或者在定义线性方程时使用PyTorch张量,然后在生成预测值时转换为numpy。 可能的错误点: - 忘记将张量从GPU移到CPU(如果使用了GPU),但用户可能只是在CPU上运行,所以暂时不需要考虑。 - 矩阵运算中的维度不匹配问题,需要确保转置和乘法的正确性。 - 正规方程在PyTorch中的实现是否正确,比如是否正确处理了浮点精度问题。 总结转换步骤: 1. 导入PyTorch库,替换numpy相关操作为PyTorch张量操作。 2. 将数据转换为张量,并添加偏置列。 3. 计算X的转置和乘积,判断秩。 4. 计算系数k,使用PyTorch的逆矩阵函数。 5. 在绘图时,将张量转换回numpy数组,保持绘图代码不变。 需要检查用户代码中是否有计算k的部分,可能在用户提供的代码中没有显示,但根据后面的定义,原代码中k应该是一个系数向量,所以需要补充这部分计算。 例如,原代码可能在使用正规方程求解线性回归参数,所以转换到PyTorch时需要做类似的操作: X_tensor = torch.tensor(X.values, dtype=torch.float32) Y_tensor = torch.tensor(Y.values, dtype=torch.float32).unsqueeze(1) XTX = X_tensor.T @ X_tensor k = torch.inverse(XTX) @ X_tensor.T @ Y_tensor 然后k就是参数向量,可能需要转置或者调整形状,具体取决于维度是否正确。 绘图时,使用k的numpy数组形式来生成预测值。 最后,确保所有张量计算正确,并且在绘图部分转换回numpy。</think>我将逐步指导你将这段代码转换为PyTorch版本,以下是修改后的代码及详细解释: ```python import torch import pandas as pd import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # 读取数据(无需修改) data = pd.read_csv('data/multiple_linear_regression_dataset.csv') # ========== 修改点1:数据转为PyTorch张量 ========== # 提取特征并添加偏置项 X = torch.tensor(data[['age', 'experience']].values, dtype=torch.float32) X = torch.cat([X, torch.ones(X.size(0), 1)], dim=1) # 添加全1列作为偏置项 Y = torch.tensor(data['income'].values, dtype=torch.float32).unsqueeze(1) # ========== 修改点2:使用PyTorch进行矩阵运算 ========== # 检查矩阵秩 XTX = X.T @ X print("是否满秩:", torch.matrix_rank(XTX) == min(X.shape)) # 计算正规方程解(当矩阵满秩时) k = torch.linalg.inv(XTX) @ X.T @ Y # 参数向量 [w1, w2, b] # ========== 绘图部分(保持matplotlib逻辑,添加张量转换) ========== def linear_equation(age, experience): # 将输入转换为张量进行计算 inputs = torch.tensor([age, experience, 1.0], dtype=torch.float32) return (inputs * k.squeeze()).sum().item() # 随机采样(保持原始逻辑) sample_size = min(50, len(data)) sample_indices = torch.randperm(len(data))[:sample_size] age_sample = data['age'].iloc[sample_indices] exp_sample = data['experience'].iloc[sample_indices] inc_sample = data['income'].iloc[sample_indices] # 创建三维图形 fig = plt.figure(figsize=(10, 8)) ax = fig.add_subplot(111, projection='3d') # 绘制散点图 scatter = ax.scatter(age_sample, exp_sample, inc_sample, c=inc_sample, cmap='viridis', marker='o', s=80, alpha=0.6) # 生成预测曲面 age_pts = torch.linspace(20, 60, 200) exp_pts = torch.linspace(0, 20, 200) age_grid, exp_grid = torch.meshgrid(age_pts, exp_pts) inc_pred = k[0]*age_grid + k[1]*exp_grid + k[2] # 转换为numpy数组供matplotlib使用 ax.plot_surface(age_grid.numpy(), exp_grid.numpy(), inc_pred.numpy(), cmap='viridis', alpha=0.2, antialiased=True) # 设置标签和视角 ax.set_xlabel('Age') ax.set_ylabel('Experience') ax.set_zlabel('Income') ax.view_init(40, 30) plt.show() ``` 主要修改点说明: 1. **数据转换**: - 使用`torch.tensor`替代numpy数组 - 通过`torch.cat`添加偏置项列 - 使用`unsqueeze(1)`将y转换为列向量 2. **矩阵运算**: - 使用`@`运算符进行矩阵乘法 - `torch.matrix_rank`替代`np.linalg.matrix_rank` - `torch.linalg.inv`替代`np.linalg.inv` 3. **绘图适配**: - 保持原始matplotlib接口 - 将最终计算结果转换为numpy数组 - 使用`.item()`提取标量值 4. **性能优化**: - 使用`torch.randperm`进行高效随机采样 - 利用广播机制进行向量化计算 需要注意的特殊情况处理: ```python # 当XTX不满秩时的处理(添加正则化项) if torch.matrix_rank(XTX) < X.shape[1]: XTX += torch.eye(X.shape[1]) * 1e-6 # L2正则化 k = torch.linalg.inv(XTX) @ X.T @ Y ``` 这个修改后的版本完整保留了原始代码的视觉效果,同时: 1. 利用PyTorch的自动微分功能(后续可扩展为神经网络) 2. 支持GPU加速(将数据移到GPU即可) 3. 与深度学习流程更好兼容 如果需要进一步与PyTorch生态系统集成,可以考虑: 1. 使用`DataLoader`进行数据加载 2. 使用`nn.Module`定义模型 3. 使用优化器进行梯度下降训练
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值