内存受限下，能否像IncrementalPCA一样批量执行PowerTransformer幂变换？

问题背景 #

• 数据集规模：(15000, 850000)，拆分存储为多个.npz（可转HDF/NetCDF）文件，每个文件包含200-1000行（时间步）与全部850000列（网格单元）
• 核心需求：先通过幂变换让各特征列接近高斯分布（依据Lim(2012)、Zhao(2022)的研究，非高斯分布会导致PCA丢失大量信息），再用IncrementalPCA做增量PCA分析
• 当前瓶颈：sklearn.preprocessing.PowerTransformer默认逐特征处理，全量加载数据需50-100GB内存，本地仅16GB内存，无法一次性完成变换

可行的批量处理方案 #

1. 用partial_fit增量拟合（sklearn 1.2+版本支持） #

这是最直接的方案，sklearn 1.2及以上版本为PowerTransformer新增了partial_fit方法，支持像IncrementalPCA那样分批拟合特征的统计量和变换参数：

from sklearn.preprocessing import PowerTransformer
import h5py
import numpy as np

# 初始化PowerTransformer，优先用Yeo-Johnson（支持正负数据，适配气候场景）
pt = PowerTransformer(method="yeo-johnson")

# 第一步：遍历所有数据文件，增量拟合变换参数
for file_path in your_data_file_list:
    # 加载单批次数据，shape=(n_samples, 850000)
    data_batch = np.load(file_path)["arr_0"].astype(np.float32)
    # 增量更新每个特征的lambda参数和统计量
    pt.partial_fit(data_batch)

# 第二步：遍历所有文件，批量执行变换并保存
for file_path in your_data_file_list:
    data_batch = np.load(file_path)["arr_0"].astype(np.float32)
    transformed_batch = pt.transform(data_batch)
    # 用HDF5存储变换后的数据，节省空间且支持分块读写
    with h5py.File(f"transformed_{file_path.split('/')[-1].replace('.npz', '.h5')}", "w") as f:
        f.create_dataset("data", data=transformed_batch, dtype=np.float32)

2. 手动分批统计+批量变换（兼容低版本sklearn） #

如果你的sklearn版本低于1.2，可以手动分两步完成：先遍历所有批次收集每个特征的必要统计量，再基于这些统计量批量变换数据：

步骤1：收集特征统计量（以Yeo-Johnson为例）

Yeo-Johnson需要每个特征的样本分位数（默认10%和90%）以及用于计算lambda的统计信息。我们可以分批累积样本总和、平方和，再合并分位数：

import numpy as np

n_features = 850000
sample_counts = np.zeros(n_features, dtype=np.int64)
sum_vals = np.zeros(n_features, dtype=np.float32)
sum_sq_vals = np.zeros(n_features, dtype=np.float32)
temp_quantile_files = []

# 遍历所有批次，累积统计量
for idx, file_path in enumerate(your_data_file_list):
    data_batch = np.load(file_path)["arr_0"].astype(np.float32)
    # 更新样本数、总和、平方和
    sample_counts += data_batch.shape[0]
    sum_vals += np.nansum(data_batch, axis=0)
    sum_sq_vals += np.nansum(data_batch ** 2, axis=0)
    # 保存当前批次的分位数到临时文件
    batch_quantiles = np.nanpercentile(data_batch, [10, 90], axis=0)
    temp_quantile_path = f"temp_quantiles_{idx}.npy"
    np.save(temp_quantile_path, batch_quantiles)
    temp_quantile_files.append(temp_quantile_path)

# 计算整体均值、标准差
mean_vals = sum_vals / sample_counts
std_vals = np.sqrt((sum_sq_vals / sample_counts) - (mean_vals ** 2))

# 合并所有批次的分位数，计算整体10%/90%分位数
all_quantiles = []
for f_path in temp_quantile_files:
    all_quantiles.append(np.load(f_path))
all_quantiles = np.concatenate(all_quantiles, axis=1)
global_quantiles = np.array([
    np.nanpercentile(all_quantiles[0], 10),
    np.nanpercentile(all_quantiles[1], 90)
])

步骤2：初始化PowerTransformer并批量变换

手动将统计量传入PowerTransformer，再逐个批次处理数据：

from sklearn.preprocessing import PowerTransformer

pt = PowerTransformer(method="yeo-johnson")
# 手动设置统计量（注意sklearn内部属性可能随版本变化，需对应调整）
pt._quantiles = global_quantiles
pt.mean_ = mean_vals
pt.scale_ = std_vals

# 拟合lambda参数（可以用一个小批次数据快速拟合）
data_sample = np.load(your_data_file_list[0])["arr_0"].astype(np.float32)
pt.fit(data_sample)

# 批量变换并保存
for file_path in your_data_file_list:
    data_batch = np.load(file_path)["arr_0"].astype(np.float32)
    transformed_batch = pt.transform(data_batch)
    # 保存到HDF5文件
    with h5py.File(f"transformed_{file_path.split('/')[-1].replace('.npz', '.h5')}", "w") as f:
        f.create_dataset("data", data=transformed_batch, dtype=np.float32)

3. 特征分块处理（极端内存紧张时用） #

如果上述方案仍有内存压力，可以将850000列特征拆分为多个小块（比如每10000列一块），对每个块单独执行批量变换：

from sklearn.preprocessing import PowerTransformer
import h5py
import numpy as np

n_features = 850000
block_size = 10000
n_blocks = (n_features + block_size - 1) // block_size

for block_idx in range(n_blocks):
    start_col = block_idx * block_size
    end_col = min((block_idx + 1) * block_size, n_features)
    
    # 增量拟合当前特征块的PowerTransformer
    pt_block = PowerTransformer(method="yeo-johnson")
    for file_path in your_data_file_list:
        data_block = np.load(file_path)["arr_0"].astype(np.float32)[:, start_col:end_col]
        pt_block.partial_fit(data_block)
    
    # 批量变换当前特征块并保存
    for file_path in your_data_file_list:
        data_block = np.load(file_path)["arr_0"].astype(np.float32)[:, start_col:end_col]
        transformed_block = pt_block.transform(data_block)
        # 追加到对应HDF5文件的特征块
        hdf_path = f"transformed_{file_path.split('/')[-1].replace('.npz', '.h5')}"
        with h5py.File(hdf_path, "a") as f:
            f.create_dataset(f"block_{block_idx}", data=transformed_block, dtype=np.float32)

关键注意事项 #

• 优先选择Yeo-Johnson变换：气候数据常包含负值（如气温距平），Box-Cox仅支持非负数据，Yeo-Johnson适配性更强
• 内存优化：全程用float32格式存储数据，可将内存占用减半；HDF5/NetCDF比NPZ更适合大规模数据的分块读写和压缩存储
• 缺失值处理：如果数据存在缺失值，需提前用np.nanmean等方式填充，或确保使用支持缺失值的统计函数（如np.nansum、np.nanpercentile）

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