Python迭代Hessian函数优化：预生成随机矩阵与计时功能改进

Hey there! Let's work through your two optimization goals for the iterative_hessian function, while keeping the out-of-core data constraint (single pass only) front and center.

1. 预生成独立随机矩阵 + 单次数据遍历 #

You’re absolutely right that re-generating sketch matrices in each iteration and re-scanning data is inefficient—especially for large datasets that can’t fit in memory. Here’s a practical fix that aligns with your requirements:

• Pre-generate independent random states: Instead of creating full sketch matrices upfront (impossible when we don’t know total sample count n for out-of-core data), we initialize separate random number generators for each iteration. This guarantees independence across sketches without wasting memory on full-sized matrices.
• Incremental preprocessing: Traverse your data chunks once, accumulating all values needed for iterations:
  • ATy (Aᵀb)
  • Covariance matrix AᵀA
  • For each iteration’s sketch, accumulate S@A and S@y (where S is the chunk-specific sketch matrix generated from the pre-seeded random state)

This way, during iterations, you never touch the original data again—all computations use preprocessed values stored in memory.

2. 基于类的可选计时功能 #

Using a class is the perfect approach here: it encapsulates timing state, lets you toggle timing on/off easily, and avoids overhead when timing isn’t needed. Here’s how it works:

• Add an enable_timing flag in the class constructor to activate/deactivate timing.
• Store individual solve step durations in a list when timing is enabled.
• Provide a method to retrieve the average solve time after running iterations.

This keeps your code clean, reusable, and efficient—no unnecessary timing calls when you don’t need them.

完整实现代码 #

import numpy as np
import time

class IterativeHessianSolver:
    def __init__(self, sketch_method, sketch_size, num_iters, enable_timing=False):
        self.sketch_method = sketch_method
        self.sketch_size = sketch_size
        self.num_iters = num_iters
        self.enable_timing = enable_timing
        self.solve_times = []
        
        # 预生成独立随机状态，保证每个迭代的 sketch 独立
        self.sketch_rngs = [np.random.RandomState() for _ in range(num_iters)]
        
    def _generate_chunk_sketch(self, n_chunk, iter_idx):
        """为单个数据块和迭代生成 sketch 矩阵"""
        m = int(self.sketch_size)
        rng = self.sketch_rngs[iter_idx]
        
        if self.sketch_method == "gaussian":
            return m**(-0.5) * rng.normal(size=(m, n_chunk))
        elif self.sketch_method == "srht":
            # 替换为你实际的 SRHT 实现（比高斯变换更快）
            # 这里是占位示例
            return m**(-0.5) * rng.randn(m, n_chunk)
        elif self.sketch_method == "sparse":
            # 稀疏随机 sketch（例如 1/m 的非零元素）
            S = rng.choice([-1, 1], size=(m, n_chunk)) * np.sqrt(m)
            S[rng.rand(m, n_chunk) > 1/m] = 0
            return S
        else:
            raise ValueError(f"不支持的 sketch 方法: {self.sketch_method}")
    
    def process_data_chunk(self, data_chunk, target_chunk):
        """增量处理单个数据块（支持核外数据）"""
        n_chunk, d = data_chunk.shape
        y_chunk = target_chunk
        
        # 第一次处理数据块时初始化累加器
        if not hasattr(self, "_initialized"):
            self.ATy = np.zeros(d)
            self.covariance_mat = np.zeros((d, d))
            self.S_dot_A_list = [np.zeros((int(self.sketch_size), d)) for _ in range(self.num_iters)]
            self.S_dot_y_list = [np.zeros(int(self.sketch_size)) for _ in range(self.num_iters)]
            self._initialized = True
        
        # 累加 Aᵀy
        self.ATy += data_chunk.T @ y_chunk
        
        # 累加协方差矩阵 AᵀA
        self.covariance_mat += data_chunk.T @ data_chunk
        
        # 累加每个迭代对应的 S@A 和 S@y
        for iter_idx in range(self.num_iters):
            S_chunk = self._generate_chunk_sketch(n_chunk, iter_idx)
            self.S_dot_A_list[iter_idx] += S_chunk @ data_chunk
            self.S_dot_y_list[iter_idx] += S_chunk @ y_chunk
    
    def solve(self):
        """执行迭代优化并返回解"""
        if not hasattr(self, "_initialized"):
            raise ValueError("还未处理任何数据！请先调用 process_data_chunk 方法。")
            
        d = self.ATy.shape[0]
        x0 = np.zeros(d)
        
        for iter_idx in range(self.num_iters):
            S_dot_A = self.S_dot_A_list[iter_idx]
            
            # 计算 B 和 z（修正了原代码的形状不匹配问题）
            B = S_dot_A.T @ S_dot_A
            z = self.ATy - self.covariance_mat @ x0 + S_dot_A.T @ (S_dot_A @ x0)
            
            # 计时（如果启用）
            if self.enable_timing:
                start = time.perf_counter()
                x_new = np.linalg.solve(B, z)
                end = time.perf_counter()
                self.solve_times.append(end - start)
            else:
                x_new = np.linalg.solve(B, z)
            
            x0 = x_new
        
        return np.ravel(x0)
    
    def get_average_solve_time(self):
        """返回每次 solve 步骤的平均耗时（仅在计时启用时有效）"""
        if not self.enable_timing:
            raise ValueError("初始化时未启用计时功能。")
        return np.mean(self.solve_times) if self.solve_times else 0.0

使用示例 #

# 初始化求解器：使用 SRHT sketch，sketch 大小 200，10 次迭代，启用计时
solver = IterativeHessianSolver(sketch_method="srht", sketch_size=200, num_iters=10, enable_timing=True)

# 处理数据块（示例：2 个数据块）
chunk1 = np.random.randn(1000, 50)  # 1000 个样本，50 个特征
target1 = np.random.randn(1000)
solver.process_data_chunk(chunk1, target1)

chunk2 = np.random.randn(1500, 50)
target2 = np.random.randn(1500)
solver.process_data_chunk(chunk2, target2)

# 运行优化
solution = solver.solve()

# 获取平均求解耗时
avg_time = solver.get_average_solve_time()
print(f"每次迭代的平均求解耗时：{avg_time:.6f} 秒")

关键细节说明 #

• 核外数据兼容: process_data_chunk 方法支持逐个处理数据块，无需将整个数据集载入内存。
• 高效 sketch 生成: 预生成随机状态保证了迭代间 sketch 的独立性，同时避免了预分配大内存的问题。
• 无 overhead 计时: 仅当 enable_timing=True 时才会进行计时，不会影响正常运行的性能。

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