请求协助计算百万级观测下(O_Temp,O_Salin)与(D_Temp,D_Salin)的Mahalanobis距离

核心原理 #

你需要计算每行中两个二维向量（(O_Temp, O_Salin) 和 (D_Temp, D_Salin)）的马氏距离，本质是对每行的差值向量计算其马氏长度。公式如下：
对于第i行的向量 $\boldsymbol{x}i = (O{Temp,i}, O_{Salin,i})$ 和 $\boldsymbol{y}i = (D{Temp,i}, D_{Salin,i})$，差值向量 $\boldsymbol{d}_i = \boldsymbol{x}_i - \boldsymbol{y}_i$，马氏距离为：
$$MD_i = \sqrt{\boldsymbol{d}_i^T \Sigma^{-1} \boldsymbol{d}_i}$$
其中 $\Sigma$ 是所有差值向量的协方差矩阵（若需用样本协方差可调整自由度，需根据业务场景明确）。

Python 高效实现（适配300万条数据） #

全程用矢量化操作避免循环，保证计算效率：

import numpy as np
import pandas as pd

# 读取数据（替换为你的数据路径）
df = pd.read_csv("your_data.csv")

# 生成差值矩阵：每行对应一组向量的差值
diff_matrix = df[["O_Temp", "O_Salin"]].values - df[["D_Temp", "D_Salin"]].values

# 计算协方差矩阵的逆矩阵（ddof=0用总体协方差，样本协方差改ddof=1）
cov_mat = np.cov(diff_matrix, rowvar=False, ddof=0)
# 处理奇异矩阵：若协方差不可逆，添加微小正则项
reg = 1e-6
cov_inv = np.linalg.inv(cov_mat + reg * np.eye(cov_mat.shape[0]))

# 矢量化计算马氏距离
md_squared = np.einsum('ij,jk,ik->i', diff_matrix, cov_inv, diff_matrix)
md_values = np.sqrt(md_squared)

# 将结果存入原数据框
df["Mahalanobis_Distance"] = md_values

R 语言高效实现 #

同样用矩阵运算替代逐行循环：

# 读取数据（替换为你的数据路径）
df <- read.csv("your_data.csv")

# 生成差值矩阵
diff_matrix <- df[, c("O_Temp", "O_Salin")] - df[, c("D_Temp", "D_Salin")]

# 计算协方差矩阵的逆（ddof=0用总体协方差，样本协方差改ddof=1）
cov_mat <- cov(diff_matrix, ddof = 0)
# 处理奇异矩阵
reg <- 1e-6
cov_inv <- solve(cov_mat + reg * diag(ncol(cov_mat)))

# 计算马氏距离
md_squared <- rowSums((diff_matrix %*% cov_inv) * diff_matrix)
md_values <- sqrt(md_squared)

# 将结果存入原数据框
df$Mahalanobis_Distance <- md_values

关键注意事项 #

• 协方差矩阵选择：若比较的是同一分布下的向量差异，用差值向量的协方差；若O、D来自不同分布，可改用两组数据的联合协方差，需根据业务逻辑调整。
• 大数据优化：避免用逐行循环（如Python的apply、R的apply），矢量化/矩阵运算能将计算速度提升几个数量级。
• 奇异矩阵处理：当两个变量完全线性相关时，协方差矩阵不可逆，必须添加正则项或对变量降维。

内容的提问来源于stack exchange，提问作者Adel KACIMI