矩阵运算、线性代数与性能

一、矩阵乘法的含义和维度

星号 A*B 表示逐元素乘法，遵守广播规则；A@B 表示矩阵乘法，二维情况下要求 A 的列数等于 B 的行数。若 X 是 n 个样本 p 个特征的 (n,p) 矩阵，系数 beta 是 (p,)，预测 X@beta 为 (n,)。每个预测是该行各特征乘系数后加总，并不是把整张矩阵相同位置两两相乘。

小例子 X=[[1,2],[1,3]]、beta=[10,5]，两行预测分别二十、二十五；第一列全一表达截距，第二列是特征。多个目标的系数矩阵 (p,k) 可得到 (n,k)。高维 @ 还涉及批次广播，初学时先限定二维并检查形状。不要通过转置试到“能运行”为止，因为维度合法仍可能把样本与特征的业务角色颠倒。

二、方程求解优先 solve

对方阵 A 和向量 b，目标是找到 x 满足 Ax=b。推荐 np.linalg.solve(A,b)，而不是先求 A 的逆再乘 b；求逆做了不必要工作，也可能放大数值误差。二维方阵必须可解，若奇异会抛 LinAlgError。求得结果后检查 A@x 与 b 的残差，不应只因为没有异常就相信每个系数都准确。

例如方程两倍 x 加 y 等于五，x 加三倍 y 等于十，解为 x 一、y 三。A 是 [[2,1],[1,3]]，b 是 [5,10]，手算消元与代码结果应一致。若第二行变成第一行的倍数，方程失去独立信息；即使 b 也成比例，有无穷多个解而非唯一解。业务特征重复会产生类似问题，不能简单删异常提示继续解释系数。

三、最小二乘使用 lstsq

大多数回归是样本多于特征的矩形系统，不存在精确满足每行的解。lstsq 找到最小化平方残差总和的系数，返回系数、残差信息、秩和奇异值。示例 x 为零、一、二，y 为一、三、五，加入截距列后拟合得到截距一、斜率二，预测完全吻合。使用 rcond=None 明确采用当前库默认阈值，不手工求 (X.T@X) 的逆。

显式形成正规方程会使条件数平方，可能降低稳定性；lstsq 的数值分解通常更合适。返回 residuals 有时为空，例如秩不足或样本数不大于特征数，空并不等价于误差为零；总可以自行计算 y-X@beta 再求平方和。解释模型时同时看秩、特征尺度和共线性，不能把一个数值解当成唯一可靠的因果结论。

四、条件数与数据尺度

条件数衡量问题对输入扰动的敏感程度，较大的值提醒系数可能对小变化很敏感。矩阵数值都有限并不等于问题稳定，接近重复的特征列就可能有很大条件数。标准化、合理特征选择或正则化可以改善问题，但这些措施也改变建模解释，需要说明。单凭任意固定阈值删除数据不是可靠统计决策。

大额和微小量相加可能丢失低位信息，float64 也有有限精度；提高 dtype 不是所有问题的万能解。计算金额总额、几何运算、线性求解分别有不同精度需求。检查残差时使用与单位相匹配的绝对与相对容差；如果数据尺度变化一百万倍，固定一微单位容差不再具有相同意义。遇到异常先检查矩阵形状、有限值、秩和尺度。

五、性能要测量完整成本

向量化矩阵乘法通常比 Python 双重循环更快，但创建大数组、复制、转置和 BLAS 线程调度都有成本。小矩阵一次运行的微秒差异不能证明某写法永久更优；使用 perf_counter 多次测量，明确是否包含数据创建和暖机。对角权重不必构造完整 n×n 对角矩阵，可以按行广播权重，避免平方内存消耗。

配套代码验证矩阵预测、solve 方程与 lstsq 回归，并展示奇异矩阵被拒绝。再用重复列构建秩不足反例，观察系数可不唯一而预测仍相同。掌握本章的标准是能够把形状、公式和检查连接起来：为什么是 @，为什么不用逆矩阵，残差如何解释，哪里可能共线或溢出。性能优化应在语义和正确性稳定之后进行，不要为了少一行代码取消必要检查。

完整示例与实际输出

import numpy as np
X=np.array([[1.,2.],[1.,3.]]); beta=np.array([10.,5.])
pred=X@beta
assert pred.tolist()==[20.,25.]
print('prediction',pred.tolist())
A=np.array([[2.,1.],[1.,3.]]); b=np.array([5.,10.])
solution=np.linalg.solve(A,b)
assert np.allclose(solution,[1,3]) and np.allclose(A@solution,b)
print('solve',solution.tolist(),'condition',float(np.linalg.cond(A)))
x=np.array([0.,1.,2.]); y=np.array([1.,3.,5.])
design=np.column_stack([np.ones_like(x),x])
fitted,residuals,rank,singular=np.linalg.lstsq(design,y,rcond=None)
error=y-design@fitted
assert np.allclose(fitted,[1,2]) and np.allclose(error,0) and rank==2
print('lstsq',np.round(fitted,6).tolist(),'rank',rank,'squared error',float(error@error))
try:
    np.linalg.solve(np.array([[1.,2.],[2.,4.]]),np.array([3.,6.]))
except np.linalg.LinAlgError:
    print('singular system correctly rejected')
else:
    raise AssertionError('singular solve should fail')

实际运行通过CPU Python execution.

prediction [20.0, 25.0]
solve [1.0, 3.0] condition 2.6180339887498953
lstsq [1.0, 2.0] rank 2 squared error 1.0970096963229695e-30
singular system correctly rejected

五道自测：先回答，再核对

两项动手任务

import numpy as np
x=np.arange(4,dtype=float); y=np.array([2.,5.,8.,11.])
X=np.column_stack([np.ones(4),x])
beta,_,rank,_=np.linalg.lstsq(X,y,rcond=None)
assert rank==2 and np.allclose(beta,[2,3]) and np.allclose(X@beta,y)
print(np.round(beta,6).tolist())

import numpy as np
X=np.array([[1.,1.],[2.,2.],[3.,3.]])
y=np.array([1.,2.,4.])
beta,residuals,rank,_=np.linalg.lstsq(X,y,rcond=None)
sse=float(np.sum((y-X@beta)**2))
assert rank==1 and residuals.size==0 and sse>0
print(beta.tolist(),rank,sse)

官方资料与掌握状态

• NumPy 线性代数
• NumPy lstsq

完成两项任务后，按独立完成程度记录。状态会尽量保存在当前浏览器，可导出备份。