使用 `drange` 拟合分布式逻辑回归

本示例演示了如何使用分布式数组和 for-drange 循环对大规模数据进行逻辑回归分析。

当数据集过大而无法容纳在客户端内存中，或者需要在一个并行池中的多个工作单元上执行数据并行计算时，可以使用分布式数组。这种方法有助于：

持久化数据访问：您可以在两次计算之间将数据保存在工作单元上，以避免来回传输大型数组。
集合运算：每个工作单元存储数组的一部分，但可以与其他工作单元通信，对整个数据集执行 sum、mean 或 norm 等归约操作。
数据并行循环：您可以在 spmd 模块内使用 for-drange-loop 高效地遍历分布式范围。在此循环中，每个工作单元独立执行其本地迭代，工作单元之间不进行通信。

本示例使用一种基于共识的交替方向乘子法 (ADMM) 算法，对合成数据拟合一个 $ℓ_{1}$ 正则化的逻辑回归模型。该实现将中间结果存储在工作单元上，并并行执行本地计算，在每次迭代之间进行集合归约以更新回归变量。分布式共识 ADMM 算法基于 Boyd 等人 [1][2] 论文中描述的串行实现方案。

在本地计算机上使用 Processes 集群配置文件启动一个并行池。如果您可以访问一个集群，则可以在该集群上启动一个并行池，方法是在 parpool 命令中将 Processes 替换为您集群配置文件的名称，并通过第二个参量指定工作单元的数量。

parpool("Processes",6);

Starting parallel pool (parpool) using the 'Processes' profile ...
Connected to parallel pool with 6 workers.

生成用于训练的合成数据

创建一个模拟大型稀疏特征矩阵的、采用模块划分方式的合成数据集。目标是仿真一个具有已知底层模型的大规模数据集，以便 ADMM 算法能够恢复用于生成该数据的权重和截距。

该数据集包含 100,000 个观测值。每个观测数据包括：

一个 1×1000 的特征向量 A
集合 {-1,1} 中的对应标签 b

定义问题规模

设置每个观测值的特征数量、数据集中的观测值总数，以及用于划分数据集的模块数量。将观测值均匀分配到各模块中，以符合分配方案。

nFeatures = 1E3;
nObs = 100000;
nBlocks = 100;
nObsPerBlock = nObs/nBlocks;

创建基准模型

定义一个稀疏权重向量和截距，以生成合成标签。该模型是 ADMM 求解器试图恢复的目标。权重向量模拟了一个具有少量有效特征的现实稀疏模型。将权重向量存储为稀疏数组，既避免了为零值分配不必要的空间，又能加快运算速度。

使用默认算法和种子设置随机数生成器。将权重向量定义为一个随机 nFeatures 行 1 列的稀疏矩阵，其中约有 100 个非零值，这些值是从标准正态分布中抽样的。将截距定义为来自标准正态分布的一个标量。

rng("default");
modelWeight = sprandn(nFeatures,1,100/nFeatures);  
modelIntercept = randn(1);

在工作单元上生成数据模块

为每个模块并行创建稀疏特征矩阵和标签，并将数据存储在工作单元上，以避免客户端内存限制。

为特征模块和标签预分配数组时，应确保其布局与 drange 及分布式计算相兼容。稀疏数组不支持超过两维的数组，因此请为特征模块使用列分布的元胞数组。此外，请使用列分布的数值数组作为标签。对两个数组采用相同的分布方案，可以简化未来的 for-drange-loop 索引操作以及 sum 和 mean 等全局归约操作。

A = distributed.cell(1,nBlocks);
b = zeros(nObsPerBlock,nBlocks,"distributed");

使用一个带有 for-drange- 循环的 spmd 模块。对于每个模块，生成稀疏观测值，添加标签噪声，计算带符号特征，并将它们存储在分布式数组中。

spmd
    for idx = drange(1:nBlocks)
        observations = sprandn(nObsPerBlock,nFeatures,10/nFeatures);
        noise = sqrt(0.1)*randn(nObsPerBlock,1);
        bi = sign(observations*modelWeight + modelIntercept + noise);
        Ai = spdiags(bi,0,nObsPerBlock,nObsPerBlock)*observations;

        A(1,idx) = {Ai};
        b(:,idx) =  bi;
    end
end

计算正则化参数

根据全局类统计信息，估计 $ℓ_{1}$ 正则化强度 mu。按照博伊德 [1] 的方法，计算 $λ_{\max}$ ，并将其按 0.1 进行缩放，以适应 ADMM 算法。

预分配分布式数组，用于存储每个模块的计数和求和结果。为确保分布式数组采用相同的分布方案，请将 nBlocks 定义为列数。

positiveCount = zeros(1,nBlocks,"distributed");
positiveSum = zeros(nFeatures,nBlocks,"distributed");
negativeSum = zeros(nFeatures,nBlocks,"distributed");

在 spmd 模块内部，使用 for-drange-loop 在各工作单元上本地处理每个模块。提取特征矩阵和标签向量，统计正标签的数量，并按标签分组计算各行的总和。

spmd
    for idx = drange(1:nBlocks)
        Ai = cell2mat(A(idx));
        bi = b(:,idx);

        % Count positive labels in this block
        positiveCount(1,idx) = sum(bi == 1);

         % Feature sums by label
        positiveSum(:,idx) = sum(Ai(bi == 1,:),1)';
        negativeSum(:,idx) = sum(Ai(bi == -1,:),1)';
    end
end

当您执行需要数组所有分区数据的全局归约操作时，分布式数组会自动管理工作单元之间的通信。计算所有模块的总计数和总和。计算数据集中正样本的比例。

totalPositiveCount = sum(positiveCount);
totalPositiveSum = sum(positiveSum,2);
totalNegativeSum = sum(negativeSum,2);
ratio = totalPositiveCount/(nObsPerBlock*nBlocks);

计算正则化参数 mu。正则化参数是正和与负和的加权组合 $λ$ 的无穷大范数的一个分数。

$λ_{\max} \approx | | (1 - r a t i o) \cdot p o s i t i v e S u m + r a t i o \cdot n e g a t i v e S u m | |_{\infty}$

计算 $λ$ ，并将其乘以 0.1/(nObsPerBlock*nBlocks) 以得到 mu.

lambda = norm((1-ratio)*totalPositiveSum + ...
ratio*totalNegativeSum,"inf");
mu = gather((0.1/(nObsPerBlock*nBlocks))*lambda)

mu = 
3.6668e-04

运行共识 ADMM 循环

使用共识 ADMM 循环，对数据集拟合 $ℓ_{1}$ 正则化的逻辑回归模型。当满足可行性容差时，停止循环。更多详情，请参阅 ADMM 算法详情。

设置 ADMM 超参数

ADMM 循环采用了 Boyd[2] 中给出的常数和容差。 $ρ$ 是惩罚（增广拉格朗日）参数， $α$ 是超松弛参数。指定最大迭代次数、绝对容差和相对容差，以控制迭代循环。

rho = 1.0;
alpha = 1.0;
maxIter = 1000;
abstol = 1e-4;
reltol = 1e-2;

执行逻辑回归

创建分布式数组，用于初始化本地截距和权重向量 x、缩放对偶向量 u 以及共识变量 z。为确保分布式数组采用相同的分布方案，请将 nBlocks 定义为列数。

x = zeros(nFeatures+1,nBlocks,"distributed");
u = zeros(nFeatures+1,nBlocks,"distributed");
z = zeros(nFeatures+1,nBlocks,"distributed");

预分配一个结构体来跟踪收敛情况。

history = struct;

使用共识 ADMM 拟合 $ℓ_{1}$ 正则化的逻辑回归模型。该过程在每次迭代中执行以下操作：

本地更新：每个工作单元都使用牛顿求解器对对数正态损失函数进行求解，从而更新其所属模块的参数。
过度放松：将本地更新与之前的共识相结合，以提高收敛稳定性。
共识更新：对各模块取平均值，并对权重应用软阈值处理，以实现 $ℓ_{1}$ 正则化。
双重更新：调整缩放后的双变量以达成共识。
诊断：跟踪残差，并在满足可行性容差时停止。

为了减少数据传输开销，分布式数组 A、b、u 和 x 在每次迭代之间仍保留在工作单元上。只有共识向量 z 会在迭代之间进行缩减，然后重新分配。updateX、saveHistory 和 plotHistory 辅助函数作为辅助文件附于本示例之后。

for step = 1:maxIter
    if step == 1
        disp("Starting ADMM loop...")
    end
    % 1) Local x-update in parallel over blocks
    spmd
        for idx = drange(1:nBlocks)
            Ai = cell2mat(A(idx));
            % solve the damped Newton subproblem for logistic loss
            x(:,idx) = updateX(Ai,b(:,idx),u(:,idx), ...
                z(:,idx),rho,x(:,idx));
        end
    end

    %  2) Over-relaxation
    zOld = z;    
    xHat = alpha*x + (1-alpha)*zOld;

    % 3) Global z-update 
    zTilde = mean(xHat+u,2); % average across blocks
    kappa = (nObsPerBlock*nBlocks)*mu/(rho*nBlocks); % threshold scalar
    a = zTilde(2:end); % L1 soft-thresholding on weights only
    zTilde(2:end) = max(0,a-kappa) - max(0,-a-kappa);
    
    % Redistribute new consensus z to all workers
    z = zTilde*ones(1,nBlocks,"distributed");

    % 4) Dual update
    u = u + (xHat-z); 

    % 5) Calculate diagnostics and stopping criteria
    zDiff = z - zOld;
    history = saveHistory(step,A,b,mu,x,z,rho,zDiff, ...
        u,nBlocks,abstol,reltol,history);
    history = plotHistory(history);

    % Stop when primal and dual residuals meet tolerances
    if history.primalRes(step) < history.EpsPrimal(step) && ...
            history.dualRes(step) < history.EpsDual(step)
        break;
    end
end

Starting ADMM loop...

提取和验证模型

收敛后，从共识变量 z 中获取拟合截距和权重向量。

interceptEst = z(1,1);
weightEst = z(2:end,1);

创建一个新的数据集来评估拟合模型，并采用相同的真实标签生成流程来生成标签。

observations = sprandn(nObsPerBlock,nFeatures,10/nFeatures);
noise = sqrt(0.1)*randn(nObsPerBlock,1);
bValid = sign(observations*modelWeight + modelIntercept + noise);
AValid = spdiags(bValid,0,nObsPerBlock,nObsPerBlock)*observations;

构建线性得分，应用逻辑映射，并预测其符号。对于验证集，您可以在 [0,1] 中获取概率，在 {-1,1} 中获取预测标签。

scores = AValid*weightEst + interceptEst;
probabilities = 1./(1+exp(-scores));
bPred = sign(scores);

计算预测的准确度。该模型的准确率约为 88%，对于带有噪声的合成数据而言，这一结果是合理的。博伊德 [1] 的原始示例在 $λ = 0.1 λ_{\max}$ 处使用更多观测数据时，准确率达到了约 90%。

accuracy = mean(bPred == bValid);
fprintf("Prediction accuracy: %.2f%%n",accuracy*100);

Prediction accuracy: 88.11%

将真实标签转换为逻辑值，计算模型残差并绘制其直方图。残差集中在零附近，表明概率校准效果更好。

yBinary = (bValid == 1);
residuals = yBinary - probabilities;

figure;
histogram(residuals,BinWidth=0.05,FaceColor=[0 0.4 0.8],FaceAlpha=0.8);
xlabel("Residual");
ylabel("Frequency");
title("Histogram of Residuals");
grid on;

ADMM 算法详情

本示例采用交替方向乘子法 (ADMM)[1] 来求解带 $ℓ_{1}$ 正则化的分布式逻辑回归问题。符号标记遵循参考文献中的规定。ADMM 用于求解如下形式的问题：

$\min I (x) + g (z) subject to x - z = 0 .$

问题表述

对于采用 $ℓ_{1}$ 正则化的逻辑回归，其目标函数为：

$\min I (x) + g (z) = \sum_{i = 1}^{m} \log (1 + \exp (- b_{i} (a_{i}^{⊤} w + v))) + μ ‖ w ‖_{1},$

受 $x - z = 0$ 的影响，权重向量为 $w \in R^{n}$ ，截距项为 $v \in R$ 。

$x_{i} = (v_{i}, w_{i})$ 是每个工作单元上的本地副本， $z = (v, w)$ 是全局共识变量， $μ$ 控制 $ℓ_{1}$ 的惩罚值。训练数据集由一组组观测数据对（ $a_{i}$ , $b_{i}$ ）组成，其中 $a_{i} \in R^{n}$ 是 $i^{th}$ 向量， $b_{i} \in {- 1, 1}$ 是相应的标签。

迭代更新

逻辑斯蒂损失函数是光滑且凸的，因此迭代更新相当于求解一个本地优化问题（牛顿法或准牛顿法），并对 $ℓ_{1}$ 项应用近端算子。每次迭代中执行的更新包括：

$\begin{array}{l} x_{i}^{k + 1} = \arg \min_{x_{i}} (\sum_{j \in D_{i}} \log (1 + \exp (- b_{j} (a_{j}^{⊤} w_{i} + v_{i}))) + \frac{ρ}{2} ‖ x_{i} - z^{k} + u_{i}^{k} ‖_{2}^{2}) \\ z^{k + 1} = (v_{}^{‾}, S_{τ} (w_{}^{‾})), x_{}^{‾} = \frac{1}{K} \sum_{i = 1}^{K} (x_{i}^{k + 1} + u_{i}^{k}), \\ u_{i}^{k + 1} = u_{i}^{k} + x_{i}^{k + 1} - z^{k + 1}, \end{array}$

其中 $ρ$ 是惩罚参数， $u_{i}$ 是归一化的对偶变量， $S_{τ}$ 是软阈值运算符：

$S_{τ} (a) = {\begin{cases} a - τ, & a > τ, \\ 0, & | a | \leq τ, \\ a + τ, & a < - τ, \end{cases}$ 与 $τ = \frac{m μ}{ρ K}$

$x_{i}$ 的本地更新是独立的，并在各个工作单元上并行计算。只有共识变量 $z$ 会在每次迭代中进行聚合和重新分配。每个本地求解器都以前一次迭代结果作为起点。随着收敛性的提高，所需的牛顿步数会减少，从而使后续的 ADMM 迭代速度更快。

残差、容差和停止准则

ADMM 监控两个残差：

原始残差： $r^{k} = {x_{i}^{k} - z^{k}}_{i = 1}^{K}, ‖ r^{k} ‖_{2} = \sqrt{\sum_{i = 1}^{K} ‖ x_{i}^{k} - z^{k} ‖_{2}^{2}} .$
双残差： $s^{k} = ρ K (z^{k} - z^{k - 1}), ‖ s^{k} ‖_{2} = ρ K ‖ z^{k} - z^{k - 1} ‖_{2} .$

在以下情况下，迭代将停止：

$‖ r^{k} ‖_{2} \leq ε_{pri}, and ‖ s^{k} ‖_{2} \leq ε_{dual}$ ,

其中

$\begin{array}{l} ε_{pri} = \sqrt{K} abstol + reltol \max {‖ x^{k} ‖_{2}, ‖ z^{k} ‖_{2}}, \\ ε_{dual} = \sqrt{K} abstol + reltol ‖ ρ u^{k} ‖_{2} \\ a b s t o l = 1 0^{- 4}, \\ r e l t o l = 1 0^{- 2} \end{array}$

参考资料

[1] 博伊德，S.“基于交替方向乘子法的分布式优化与统计学习。”《计算机学习的基础与趋势®》 第 3 卷第 1 期（2010 年）：1–122. https://doi.org/10.1561/2200000016.

[2] “基于交替方向乘子法的分布式优化与统计学习代码。”访问日期：2025 年 11 月 10 日。https://web.stanford.edu/~boyd/papers/admm_distr_stats.html。

另请参阅

for (drange) | distributed

使用 drange 拟合分布式逻辑回归